/*
* 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;
}