/* file ephem, solar ephemeris routines */
#include<stdio.h>
#include<string.h>
#include <math.h>
 static	float	rq = 57.2957795, rqi = 0.0174532925;
 static	float	pi = 3.141592654, b, r, d, p, cl, crotation;
 /*--------------------------------------------------------------------------*/
int julian(iy,im,id)
 int	iy,im,id;
 /*  (j.meeus: astr. formulae for calculators, willmann bell inc. 1982)
 c   parameters in :
 c         iy  -  year no. after 1900 
 c         im  -  month no. of the year    
 c         id  -  day no. of the month  
 c   parameters out :
 c         jd  -  julian day
 */
 /* this returns an integer Julian day for the start of the UTC day, you
 should add 0.5 to the result to get the actual Julian time for the
 start (00:00 UTC) of the day entered */
 {
 double	b, did, a;
 int	ii, iyy, ij;
 iyy  =  1900. + iy;
 if (im <= 2) { iyy = iyy - 1;	im = im+12; }
 a = ((float) iyy)/100.;	b = 2.-a+a/4.;	did = id;
 ii = 365.25 * (float) iyy;	ij = 30.6001* (float) (im+1);
 return	(int) (ii + ij + b + 1720994.5 + did);
 }
 /*--------------------------------------------------------------------------*/
int sephem_year_day(ny,day)
 int	ny;
 float	day;
 {
 /*	computes solar b angle (in radians) and solar radius (in arcsec)
	 input is year and day of year (including fraction of day)*/
 /*
   time variable for newcomb's formulae:
   fraction of julian centuries elapsed since 1900.05 (noon 1st of january)
   =  j.d.2415020.0) 
 */
 double	sday, jd, h0, h, hh, ehel, eks, sml, anm, cc, el, sl, san, av, om, ba;
 double	year, eincl, t, theta;
 int	id, im, idoy, iy;
 /* the day is done as a fp value such that the first day is 0 to 0.9999999,
 this means that to use the julian date function, we have to add 1,
 we finally get the f.p. day since 1900.05 */
 idoy = (int) (day+1.0);	iy = ny;	sday = (day- (int) day);
 admo(idoy, iy, &id, &im);
 /* printf("iy,im,id = %d %d %d\n",iy,im,id); */
 jd = (double) julian(iy,im,id)+0.5+sday;
 /* need years since 1850 to calculate long. of ascending node */
 year = 1900.+ny+(day/365.);
 return sephem(jd, year);
 }
 /*--------------------------------------------------------------------------*/
int sephem(jd, year)
 double	jd, year;
 {
 double h0, h, hh, ehel, eks, sml, anm, cc, el, sl, san, av, om, ba;
 double	eincl, t, theta;
 //printf("jd = %g, year = %g\n", jd, year);
 h0 = (jd - 2415020.0)/36525.0;
 h = (jd - 2415020.0)/36525.0;	hh = h * h;
 //printf("sday, jd, h0, h, hh = %f, %10.1f %g %g %g\n",sday, jd, h0, h, hh);
	 /* newcomb's formulae. (page 98 explanatory suppl. to the ephemeris)
	    mean obliquity of the ecliptic*/
 ehel = 0.4093198 - 2.2703e-4 * h - 2.86e-8 * hh;
					 /* eccentricity of earth's orbit */
 eks = 0.01675104 - 0.0000418*h;
						 /* mean longitude of sun*/
 sml = 279.6967 + 36000.769*h;
 sml = sml - 360.* floor( sml/360.);
 sml = sml * rqi;
 /*printf("sml = %f\n",sml);*/
						 /* mean anomaly*/
 anm = 358.4758+35999.0498*h - 0.00015*hh;
 anm = anm*rqi;
					 /* true longitude of sun (sl)*/
 cc = (1.91946-0.00479*h) * sin(anm) + 0.020 * sin(2*anm);
 cc = cc * rqi;
 sl = sml + cc;
							 /* true anomaly*/
 /*printf("anm, cc = %f %f\n", anm, cc);*/
 san = anm + cc;
 el = sl;
					 /* no correction for abberation */
					 /* distance to sun*/
 /*printf("san, eks = %f %f\n",san,eks);*/
 av = (1-eks*eks)/(1+eks*cos(san));
							 /* radius of sun*/
 r = atan(0.0046555/av)*(rq)*3600.;
							 /* in arcseconds */
 /* need years since 1850 to calculate long. of ascending node*/
 /* longitude of ascending node (page 171 smart)*/
 om = (73.66666+0.01395833*(year-1850.0)) * rqi;
 /* inclination of solar equator on the ecliptic (7.25 degrees ('constant'))*/
 eincl = 0.12653637;
 
 /* B angle, heliographic latitude of centre of disc*/
 ba = sin(el-om) * sin(eincl); b = asin(ba) * rq;
 /* D angle */
 ba = sin(el-om+0.5*pi) * sin(eincl); d = asin(ba) * rq;
 /* P angle, position angle of northern rotation pole of sun*/
 t = atan(-1.0*cos(el)*tan(ehel))+atan(-1.0*cos(el-om)*tan(eincl));
 p = t * rq;					/* p angle in degrees*/
 /* Carrington longitude (parts taken from LS's code) */
 /* rotation of sun since reference back in 1853, note sideral rate here */
 theta = rqi*(jd - 2398220.0)*360.0/25.38;
 cl = fmod(atan2(-sin(el-om)*cos(eincl), -cos(el-om)) - theta, 2.*pi)*rq;
 if (cl < 0.0) cl += 360.;
 return 1;
 }