/*
* Program calculates the analemma as a function of the sun's declination
* for each hour and projects the results onto a user defined plane surface
*
* Two methods are used for computation of the true anomaly.

* Influence of other planets, parallax, refraction, precession
* and aberration are not taken into account.
*
* There is an implicit assumption that the sun and the mean sun have the same right ascension at perigee, ie. the
* sun and mean sun are in synch at perigee. Given that the right ascension of the mean sun is a calculated
* number via Newcomb's tables this assumption seems to need some confirmation.
*
* 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; /* angle of vernal equinox with respect to perihelion (in degrees) (epoch 1990) */
double e=.016713; /* eccentricity of earth/sun elipse (epoch 1990) */
double R=360/365.242191; /* angle of movement in siderial day (in degrees) (epoch 1990) */
double obl=23.441884; /* obliquity of ecliptic (in degrees) (epoch 1990) */

double eps=.00000001; /* precision used for calculation of anomaly */
double F=1; /* fraction of day 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(),project(),rotateZ(),rotateY();
double getV(),floor(),sin(),cos(),tan(),sqrt(),atan(),fabs();
double piover180;

int main(void)
{
 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();

 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 sinObl;
 double cosObl;
 double del;
 int j;

 piover180=PI/180;

 Wr=W*piover180;

 sinObl=sin(obl*piover180);
 cosObl=cos(obl*piover180);

 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 longitute of sun */

/* convert from eliptic coordinates to equatorial coordinates */

  Rx=cos(Lr);
  Ry=cosObl*sin(Lr);
  Rz=sinObl*sin(Lr);

/* Rx, Ry, Rz are the position coordinates of the sun in the equatorial system */

  getAngles(Rx,Ry,Rz,&alpha,&delta);

  alphaM=M+W; /* 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 for each hour and output values to file */
  for(j=-12;j<=12;j++)
   {
   project(diff+j*15,delta);
   }
 }
 return;
}

double del=0; /* displacement of local meridian from local noon time zone meridian in degrees */
double dec=0; /* local declination in degrees */

/*
unit vector defining the projection plane
horizontal plane
Rxp,Ryp,Rzp = 0,0,1

equatorial plane
Rxp,Ryp,Rzp = sin(dec*piover180),0,cos(dec*piover180)

vertical plane facing south
Rxp,Ryp,Rzp = 1,0,0

vertical facing east of south by 20 degrees
Rxp,Ryp,Rzp = cos(20*piover180),sin(20*piover180),0
*/

double Rxp=0; /* unit vector defining the projection plane */
double Ryp=0;
double Rzp=1;

void project(diff,delta)
double diff,delta;
{
 double diffr,pdSunr;
 double Rx,Ry,Rz;
 double delr,pdr;
 double alphap,deltap;
 double x,y;
 double scale;

/*
consider a right handed earth centered coordinate system with the z axis toward the north pole
the x axis toward the local noon time zone meridian and the y axis perpendicular to the other two

the positions of the sun during the year with respect to this reference meridian and the equator
are given by diff and delta (pd is polar distance)
the position vector of the sun with repect to this system is
*/

 diffr=diff*piover180;
 pdSunr=(90-delta)*piover180;

 Rx=cos(diffr)*sin(pdSunr);
 Ry=sin(diffr)*sin(pdSunr);
 Rz=cos(pdSunr);

/*
Let the local meridian be measured with repect to this noon meridian with a displacement toward the east
being a positive number. Let the local polar distance pd be 90-dec where dec is the local declination.

If del is that displacement then a rotation about the z axis by -del and a rotation about the y axis by
-pd gives the position vector of the sun with respect to the local position on the earth.
*/

 delr=del*piover180;
 pdr=(90-dec)*piover180;

 rotateZ(&Rx,&Ry,&Rz,-pdr);
 rotateY(&Rx,&Ry,&Rz,-delr);

/*
now consider the local position vector of the sun reversed in direction
and follow that vector from the sun to the origin of the coordinate system to some arbitrary
surface. The surface might be a horizontal plane, a vertical plane with some orientation,
a cube, a cylinder, a hemi-sphere, etc.

the resulting projections of the sun's positions at each hour onto the surface of choice can represent
analemmic time lines for a sundial

a small sphere at the origin of the coordinate system can cast a shadow to indicate the time

as an example, we calculate the projection of the suns positions at noon onto a plane defined
by a unit vector perpendicular to the plane

this problem has meaning to someone designing a sundial to be placed on a slanting wall
oriented in an abitrary direction

we can think of the system as a unit sphere centered on the origin of the coordinate system
with the plane tangential to the unit sphere in a direction defined by the negative of the vector

Rxp,Ryp,Rzp is given in the local coordinate system

for Rxp=Ryp=0, the surface is a horizonal plane
if Rzp= 1 the plane is below the sphere, if Rzp= -1 the plane is above the sphere

for Ryp=Rzp=0, the surface is a vertical plane oriented toward the south or north
depending upon the sign of Rxp

for Rxp=Rzp=0, the surface is a vertical plane oriented toward the the east or west
depending upon the sign of Ryp

we want to rotate the local sun vector to the coordinate system defined by the point of intersection
of the plane's unit vector with the unit sphere

first get the angles involved in this rotation
*/

 getAngles(Rxp,Ryp,Rzp,&alphap,&deltap);

/* rotate point */

 rotateZ(&Rx,&Ry,&Rz,-alphap*piover180);
 rotateY(&Rx,&Ry,&Rz,-deltap*piover180);

/* project the point onto the plane */

/*
the origin of the x,y system on the plane is the point of intersection of the position
vector that defines the plane

Rz is the cos of the projection angle

tan of that angle is the length of the projection vector
*/

 if(Rz<=0)
  {
  return;
/* ignore points that project to infinity or project onto the negative side of the plane */
  }

 x= -Rx/Rz; /* take the negative of the vector during the projection */
 y= -Ry/Rz;

 scale=600; /* set scale as needed */
 x*=scale;
 y*=scale;

 fprintf(fp,"%lf,%lf\n",x,y); /* output the point for use on dial */

return;
}

void rotateY(pRx,pRy,pRz,del) /* rotate vector about the Y axis by del in radians */
double *pRx,*pRy,*pRz;
double del;
{
 double x,z,cosDel,sinDel;

 cosDel= cos(del);
 sinDel= sin(del);

 x= cosDel*(*pRx)+sinDel*(*pRz);
 z= -sinDel*(*pRx)+cosDel*(*pRz);
 *pRx= x;
 *pRy= *pRy;
 *pRz= z;

 return;
}

void rotateZ(pRx,pRy,pRz,del) /* rotate vector about the Z axis by del in radians */
double *pRx,*pRy,*pRz;
double del;
{
 double x,y,cosDel,sinDel;

 cosDel= cos(del);
 sinDel= sin(del);

 x= cosDel*(*pRx)-sinDel*(*pRy);
 y= sinDel*(*pRx)+cosDel*(*pRy);
 *pRx= x;
 *pRy= y;
 *pRz= *pRz;
 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;
}