/*
/*
* Program to calculate the analemma as a function of the sun's declination.
*
* Extreme values for the difference in right ascension or hour angle of the sun
* and mean sun are computed.
*
* Two methods are used for computation of the true anomaly.
*
* Approximate method gives extremes of -16 min 20 sec and +14 min 11 sec
* Exact method gives extremes of -16 min 24 sec and +14 min 16 sec
*
* Influence of other planets, parallax, refraction, precession
* and aberration are not taken into account.
*
* The notion of a mean sun involves a value for its RA during the year that
* on the average is as close to that of the sun as possible. One method of determining
* the location of the mean sun is to estimate its RA for perigee, compute the EoT for the
* year and sum the values of the EoT. Ideally the sum should be zero. If it is not
* one may alter the value of the RA and perform the calculation again. Repeated such
* calculations allows one to converge on a value for the RA of the mean sun at
* perigee that gives a value of zero for the some of the EoT.
* That process leads to a value that is numerically very close
* to the numerical value of the ecliptic longitude for the sun at perigee.
* The actual value for the RA of the sun at perigee that gives a zero sum of the EoT differs
* from W by 0.00039735 degrees.
*
* Reference for values and more information:
* "Practical Astronomy with your calculator", Third Edition, Peter Duffett-Smith
*/
#include <stdio.h>
#include <stdlib.h>
double W=282.768422; /* ecliptic angle of vernal equinox with respect to perihelion (in degrees) (epoch 1990) (ecliptic longitude) */
double e=.016713; /* eccentricity of earth/sun elipse (epoch 1990) */
double obl=23.441884; /* obliquity of ecliptic (in degrees) (epoch 1990) */
double RA0=282.768422; /* value of the RA of the mean sun at perigee that zeros the sum of EoT over the year */
double R=360/365.242191; /* angle of movement in siderial day (in degrees) (epoch 1990) */
/* this value is not needed in the calculation but is used to give some relationship to the movement of the sun in a day */
double eps=.00000001; /* precision used for calculation of anomaly */
double F=.1; /* FR is the interval between output values */
int method=1; /* method=0 approximation for v, method!=0 use more accurate method */
double PI=3.14159265358979323846; /* value of pi */
FILE *fp;
void doCalculation(),getAngles();
double getV(),floor(),sin(),cos(),tan(),sqrt(),atan(),fabs();
double maxDiff,minDiff;
int main(void)
{
int min,sec;
printf("Starting calculation with the following values --\n\n");
printf("W=%lf\n",W);
printf("e=%lf\n",e);
printf("R=%.8lf\n",R);
printf("obl=%lf\n",obl);
printf("\n");
printf("eps=%.8lf\n",eps);
printf("F=%lf\n",F);
if(method==0)printf("Using approximate calculation of E\n\n");
else printf("Using more accurate calculation of E\n\n");
if((fp=fopen("output","w"))==NULL)
{fprintf(stderr,"Cannot open file (output)!\n");exit(1);}
doCalculation();
min=floor(maxDiff*60/15); /* calculate diff in minutes and seconds */
sec=(maxDiff*60/15-min)*60;
printf("Extreme differences are plus %d minutes %d seconds and -",(int)min,(int)sec);
minDiff=-minDiff;
min=floor(minDiff*60/15); /* calculate diff in minutes and seconds */
sec=(minDiff*60/15-min)*60;
printf("%d minutes %d seconds\n",min,sec);
fclose(fp);
return(0);
}
void doCalculation()
{
double M; /* mean anomaly in degrees */
double Mr; /* mean anomaly in radians */
double Vr; /* true anomaly in radians */
double Lr; /* ecliptic longitude */
double Wr; /* ecliptic longitude of perigee */
double alpha; /* right ascension of sun with respect to vernal equinox */
double delta; /* declination of sun with respect to plane of equator */
double alphaM; /* right ascension of mean sun */
double Rx,Ry,Rz; /* vector pointing to sun in equatorial coordinates */
double diff; /* difference in right ascension of sun and mean sun */
double delM; /* to adjust W to give a zero sum of EoT over the year */
double sinObl;
double cosObl;
double del,piover180;
int min,sec;
piover180=PI/180;
Wr=W*piover180;
sinObl=sin(obl*piover180);
cosObl=cos(obl*piover180);
maxDiff=minDiff=0;
del=F*R; /* increment in mean anomaly */
for(M=0;M<=360;M+=del)
{
Mr=piover180*M;
if(method==0) Vr=Mr+2*e*sin(Mr);
else Vr=getV(Mr);
Lr=Vr+Wr; /* ecliptic longitude of sun */
/* convert from eliptic coordinates to equatorial coordinates */
Rx=cos(Lr);
Ry=cosObl*sin(Lr);
Rz=sinObl*sin(Lr);
getAngles(Rx,Ry,Rz,&alpha,&delta);
alphaM=M+RA0; /* right ascension of mean sun in equatorial coordinates */
diff=alpha-alphaM; /* difference in right ascension between sun and mean sun */
if(diff<0)diff=diff+360;
/* compute extreme values of difference */
if(diff>maxDiff)maxDiff=diff;
if(diff<minDiff)minDiff=diff;
/* output values to file */
fprintf(fp,"%lf %lf ",diff,delta);
if(diff<0)
{
diff=-diff;
min=floor(diff*60/15); /* calculate diff in minutes and seconds */
sec=(diff*60/15-min)*60;
fprintf(fp,"-%d:%d\n",min,sec);
}
else
{
min=floor(diff*60/15); /* calculate diff in minutes and seconds */
sec=(diff*60/15-min)*60;
fprintf(fp,"%d:%d\n",min,sec);
}
}
return;
}
double getV(Mr)
double Mr;
{
double Er,delEr,del,Vr,tanVover2;
/* solve the equation Er-e*sin(Er)=Mr (Kepler's equation) */
/* Er is the eccentric anomaly in radians */
Er=Mr; /* set initial approximation for Er */
while(1)
{
del=Er-e*sin(Er)-Mr;
if(fabs(del)<eps)break;
delEr=del/(1-e*cos(Er));
Er=Er-delEr;
}
/* compute Vr from Er */
tanVover2= tan(Er/2)*sqrt((1+e)/(1-e));
Vr= atan(tanVover2)*2;
return(Vr);
}
void getAngles(Rx,Ry,Rz,pAlpha,pDelta)
double Rx,Ry,Rz;
double *pAlpha,*pDelta;
{
double RAr,PDr,cosPD,tanPD;
if(Rx != 0)
{
RAr= atan(Ry/Rx);
if(RAr > 0 && Ry < 0 && Rx < 0)
RAr= PI+RAr;
if(RAr < 0)
{
if(Ry < 0 && Rx > 0)
RAr= 2*PI+RAr;
if(Ry > 0 && Rx < 0)
RAr= PI+RAr;
}
if(RAr == 0 && Rx < 0)
RAr= PI;
}
else
{
if(Ry != 0)
{
if(Ry > 0)
RAr=PI/2;
else
RAr= 3*PI/2;
}
else
RAr= 0;
}
cosPD=Rz;
if(cosPD != 0)
{
tanPD= sqrt(1-cosPD*cosPD)/cosPD;
PDr= atan(tanPD);
if(cosPD < 0)
PDr= PI+PDr;
}
else
PDr= PI/2;
*pDelta= 90-PDr*180/PI;
*pAlpha= RAr*180/PI;
return;
}