/***
 *  $Header: /home/tom/C/astro/RCS/tropical_year.c,v 1.3 2001/11/21 23:24:33 tom Exp tom $
 *
 *	tropical_year
 *
 * Compute the length of the tropical year, and its rate,
 * as function of the mean anomaly of the reference point on the ecliptic
 * and of the time since the epoch.
 *
 * A.R. (Tom) Peters
 * Created 20011112..
 * Version 20011121
 *
 * Source code license: GPL
 */

/*
Rather than doing this purely numerically, I do a semi-analytical derivation of the year length and rate of change for instructive purposes.

Formulae:

time:
	t

Excentricity:
	 E = E0 + E1*t +   E2*t^2 +   E3*t^3 +   E4*t^4 +   E5*t^5
    change:
	dE =      E1   + 2*E2*t   + 3*E3*t^2 + 4*E4*t^3 + 5*E5*t^4

Mean Ecliptic Longitude of Perihelion:
	 P = P0 + P1*t +   P2*t^2 +   P3*t^3 +   P4*t^4 +   P5*t^5
    change:
	dP =      P1   + 2*P2*t   + 3*P3*t^2 + 4*P4*t^3 + 5*P5*t^4

Mean Anomaly:
	 M = M0 + M1*t +   M2*t^2 +   M3*t^3 +   M4*t^4 +   M5*t^5
    mean orbital angular velocity:
	dM =      M1   + 2*M2*t   + 3*M3*t^2 + 4*M4*t^3 + 5*M5*t^4

Equation of Center:
	 C = C1*sin(M) + C2*sin(2*M) + C3*sin(3*M) + C4*sin(4*M) + ...
		C1 = 2*E - (1/4)*E^3 + ...
		C2 = (5/4)*E^2 - (11/24)*E^4 + ...
		C3 = (13/12)*E^3 - ...
		C4 = (103/96)*E^4 - ...
    derivative:
	dC = dM*X
	 X = C1*cos(M) + 2*C2*cos(2*M) + 3*C3*cos(3*M) + 4*C4*cos(4*M) + ...

Mean Ecliptic Longitude:
	 L = M + P
	   = L0 + L1*t +   L2*t^2 +   L3*t^3 +   L4*t^4 +   L5*t^5
    mean angular velocity:
	dL = dM + dP
	   =      L1   + 2*L2*t   + 3*L3*t^2 + 4*L4*t^3 + 5*L5*t^4

True Anomaly:
	 A = M + C
    true orbital angular velocity:
	dA = dM + dC = dM*(1+X)

Mean Tropical Year:
	J = 2*pi/dL = 2*pi/(dM+dP)

Anomalistic Year:
	K = 2*pi/dM

Mean difference:
	DJK = J-K = -2*pi*dP/(dM*(dM+dP)) = -J*dP/dM

Difference in orbit angle after 1 mean tropical year:
	DH = dM*J - dL*J = -dP*J = -2*pi*dP/(dM+dP)

True time to run the difference:
	DT = -dP*J/dA = -2*pi*dP*(1/(dM+dP))*(1/dA)

1/dA = 1/(dM*(1+X)) = (1/dM)*(1 - X + X^2 - X^3 + X^4 - ...)
     = (1+Q)/dM
 Q = -C1*cos(M) -2*C2*cos(2*M) -3*C3*cos(3*M) -4*C4*cos(4*M)
     +C1^2*cos(M)^2 +4*C1*C2*cos(M)*cos(2*M) +6*C1*C3*cos(M)*cos(3*M) +4*C2^2*cos(2*M)^2
     -C1^3*cos(M)^3 -6*C1^2*C2^2*cos(M)^2*cos(2*M)
     +C1^4*cos(M)^4
dQ = dM*Z
 Z = +C1*sin(M) +4*C2*sin(2*M) + ...
     -2*C1^2*cos(M)*sin(M) - ...

So:
	DT = -J*(dP/dM)*(1+Q)

Therefore the difference between the mean tropical year, and a tropical year for an ecliptic reference point with ecliptic longitude L and mean anomaly M, is:

	DU = DT - DJK = -J*(dP/dM)*Q = -2*pi*dP*(1/dM)*(1/dL)*Q

Time derivative:
*! These appear to be in error !*
	d(J)/dt = dJ = 2*pi*[ -(1/(dM+dP))*(1/(dM+dP)) * (d2M+d2P) ]
	             = -J*(d2M+d2P)/(dM+dP) = -J*d2L/dL
	d(DU)/dt = dU = -2*pi*[
		+d2P*(1/dM)*(1/dL)*Q
		-dP*(1/dM)^2*d2M*(1/dL)*Q
		-dP*(1/dM)*(1/dL)^2*d2L*Q
		+dP*(1/dM)*(1/dL)*dQ
		             ]
		      = J*[
		-d2P*(1/dM)*Q
		+dP*(1/dM)^2*d2M*Q
		+dP*(1/dM)*(1/dL)*d2L*Q
		-dP*(1/dM)*Q
		          ]
		The first 4 components are negligeable.
	approximations:
		d(E)/dt = dE = 0
		d(P)/dt = dP = P1
		d^2(P)/dt^2 = d2P = P2
		d(M)/dt = dM = M1
		d^2(M)/dt^2 = d2M = M2
		d(L)/dt = dL = L1
		d^2(L)/dt^2 = d2L = L2
		d(Q)/dt = dQ = dM*Z (to 2nd order)
	then:
	           dU = -2*pi*[
		+P2*(1/M1)*(1/L1)*Q
		-P1*(1/M1)^2*M2*(1/L1)*Q
		-P1*(1/M1)*(1/L1)^2*L2*Q
		+P1*(1/L1)*Z
		              ]
		      = J*[
		-P2*(1/M1)*Q
		+P1*(1/M1)^2*M2*Q
		+P1*(1/M1)*(1/L1)*L2*Q
		-P1*Z
		          ]
*/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <ctype.h>
#include <string.h>

#undef ANALYTICAL

/* some constants */
#define DAY     86400.0   	/* day in (SI) seconds */
#define GYR       365.2425	/* average Gregorian year in days */
#define JYR       365.25  	/* average Julian year in days */
#define JCY     36525.0   	/* 100 Julian years in days */
#define JMIL   365250.0   	/* 1000 Julian years in days */
#define J1900 2415020     	/* Julian day of epoch J1900 */
#define J2000 2451545     	/* Julian day of epoch J2000 */

const double vlight_d = 5.775518331E-3;	/* 1/c (light speed), in days per AU */
const double a_earth = 1.0000010178;	/* mean orbit radius of Earth, A.U. */

char *DoDebug;

/* Orbit parameters from:
 * J.L.Simon et al., Astron.Astroph. 282, 663..683 (1994)
 */
/* Excentricity */
double E[] = { 0.0167086342, -420.3654E-6, -12.6734E-6, 144.4E-9, -0.2E-9, 0.3E-9 };
/* Ecliptic Longitude of the Perihelion to Mean Equinox of Date */
double P[] = { 1.7965956473, 0.3001023491, 797.41170E-6, -308.58E-9, -586.14E-9, 14.45E-9 };
/* Mean Anomaly */
double M[] = { 6.2400601194, 6283.01955171586, -268.19834E-6, 657.98E-9, -554.63E-9, -23.17E-9 };
 
/*
The excentricity, longitude of the perihelion (to a fixed frame of reference), and the precession, actually have few secular terms; rather, they vary with long-period terms around some mean value.  For a list, see for instance:
	A.L. Berger (1977): Celest.Mech. 15, 53..74
	A.L. Berger (1979): J.Atmos.Sc. 35, 2362..2367
I prefer using an approximation by a polynomial in time rather than a Fourier series because:
- it is more convenient in the analytical expressions, especially differentiation;
- the expressions of Simon e.a. are more accurate for several thousands of years around present.
However, for a full precession period, periodic terms like from Berger should be used.
*/


/***
 * Mathematical utility functions.
 * T.P. created 199?
 * Modified 950804
 */

#define SGN(x)	((x)<0?-1:((x)==0?0:+1))

typedef struct {
                unsigned h;
                unsigned m;
                float    s;
                int      t;
               } SEXAGES ;

long DIV(long num, long den)
  {
  /*
   *  q=num/den  is defined only for num>0, den>0 ;
   * DIV does not truncate, but rounds down.
   * Consistent with MOD.
   */
  long q;
 
  /* avoid ambiguities, use positive operands only */
  if (den<0) {
    if (num<0) {
      q = (-num)/(-den) ;
    } else {
      q = -(num/(-den)) ;
      if (num%(-den) != 0) q -= 1 ;
      }
  } else if (den>0) {
    if (num<0) {
      q = -((-num)/den) ;
      if ((-num)%den != 0) q -= 1 ;
    } else {
      q = num/den ;
      }
  } else { /* den==0 */
    q = SGN(num)*LONG_MAX ;
    }
 
  return(q);
  }

long MOD(long num, long den)
  {
  /*
   *  m=num%den  defined only for num>=0, den>0 ;
   * MOD returns 0<=m<den if den>0 ;
   *             0>=m>den if den<0 .
   * Consistent with DIV .
   */
  long m;
 
  /* avoid ambiguities, use positive operands only */
  if (den<0) {
    if (num<0) {
      m = -( (-num)%(-den) ) ;
    } else {
      m = -(   ( -( den + (num%(-den)) )  )   %   (-den)   ) ;
      }
  } else if (den>0) {
    if (num<0) {
      m = ( den - ((-num)%den) )  %  den ;
    } else {
      m = num % den ;
      }
  } else if (den==0) {
    m = num ;
    }
 
  return(m);
  }

long TRUNC(double f)
  /* identical to modf(), integer part only */
  {
  double i;
  (void)modf(f,&i);
  return((long)i);
  }

long INT(double f)
  {
  double i,r;
  r=modf(f,&i);
  if (r!=0.0 && f<0.0) i -= 1.0 ;
  return((long)i);
  }

double RINT(double f)
  {
  double i,r;
  r=modf(f,&i);
  if (r!=0.0 && f<0.0) i -= 1.0 ;
  return((double)i);
  }

double RTRUNC(double f)
  /* identical to modf(), integer part only */
  {
  double i;
  (void)modf(f,&i);
  return((double)i);
  }

double ROUNDR(double f)
  {
  double i;
  i=RINT(f+0.5);
  return(i);
  }

double RMED(double num, double den)
  {
  /*
   *  RMED  returns -den <= m <= den
   */
  double m,dum;

  m=num;
  if (fabs(num)>fabs(den)) {
    dum=num/(den*2);
    m=num-ROUNDR(dum)*den*2;
  }
  return(m);
  }

double RMOD(double num, double den)
  {
  /*
   * standard function  fmod()  is not always appropriate.
   *  RMOD  returns 0<=m<den if den>0 ;
   *                0>=m>den if den<0 .
   */
  double m,dum;

  dum=num/den;
  m=num-RINT(dum)*den;
  return(m);
  }

SEXAGES *R_SEXA(double f, SEXAGES *hms)
  {
  double  dh,dm;
 
  hms->t=SGN(f);
  f=fabs(f);
  f=(60.0*modf(f,&dh));
  hms->h=(unsigned)dh;
  hms->s=(float)(60.0*modf(f,&dm));
  hms->m=(unsigned)dm;
 
  return(hms);
  }


/***
 *  Julian Day conversion Algorithms from :
 *
 *  D.A.Hatcher : Generalized Equations for Julian Day Numbers
 *                and Calendar Dates
 *  Quarterly Journal of the Royal astronomical Society
 *    (1985) 26, 151-155
 *
 * T.P. created 199?
 * Modified 950804
 * Adapted for calendar type 20011120
 */

void datjul(double jd, char *Cp, int *y, int *m, double *d)
{
  long J1,g;
  int  Y1,M1,M2,T1,D1,K,L;
  double f;

  /* first correction for Gregorian Calendar */
  if ( (*Cp=='G') || (*Cp=='g') ||
       ( ((*Cp=='C')||(*Cp=='c')) && (jd >= 2299160.5) )
     )
  {
    g = ( (long)floor((jd+68569.5)/36524.25) )*3L/4L - 38L ;
    *Cp='G';
  } else {
    g=0 ;
    *Cp='J';
  }

  J1 = (long)floor(jd+1401.5) + g ;
  f  = (jd+0.5)-floor(jd+0.5);
  Y1 = (int)DIV(J1*4L+3L,1461L);   /* integer divisions also */ 
  K  = (int)MOD(J1*4L+3L,1461L);   /* for negative numbers ! */
  T1 = K/4;
  M1 = (T1*5+2)/153;
  L  = (T1*5+2)%153;
  D1 = L/5;
  M2 = (M1+2)%12 + 1 ;

  if (d!=NULL) *d = (double)D1+1.0+f;
  if (m!=NULL) *m = M2 ;
  if (y!=NULL) *y = Y1-4716+(14-M2)/12;
}


/***
 * Application functions
 */

double analytical(double Linp, int O_A, double yinp, double *tya, double *dtya)
{
  int i,j;
  double y,t,r,tab,jda;
  double corr;
  double Ph,Ma,Mb,Ea,A,Lo,Ex;
  double C,C1,C2,C3,C4,X,Q,Z;
  double dP,dM,dL,dC,dQ;
  double d2P,d2M,d2L;
  double J,K,DH,DT,DKJ,DU,dU,dJ;
  double P1,P2,M1,M2,L1,L2;
  double F1,F2,F3,F4;

  if (DoDebug) fprintf(stderr,"D: begin analytical\n");

  t = (yinp-2000.0)*0.001; /* time since epoch J2000 in 1000's of Julian years */

  /* Initial estimates */

  for (i=5, Ex=0.0; i>0; i--) Ex = (Ex+E[i])*t ; Ex+=E[0];
  for (i=5, Ph=0.0; i>0; i--) Ph = (Ph+P[i])*t ; Ph+=P[0];
  Ph += M_PI;   /* Want longitude of Sun, not Earth */
  /* Mean Anomaly for estimated time t : */
  for (i=5, Mb=0.0; i>0; i--) Mb = (Mb+M[i])*t ; Mb+=M[0];
  Mb = RMOD(Mb,M_PI*2);
  for (i=5, dM=0.0; i>1; i--) dM = (dM+i*M[i])*t ; dM+=M[1];

  /* Use Keplers equations forward */
   A = RMOD(Linp-Ph,M_PI*2);   /* estimated True Anomaly */
  Ea = 2*atan2( sin(A*0.5), cos(A*0.5)*sqrt((1+Ex)/(1-Ex)) );
  Ma = Ea - Ex*sin(Ea);   /* target Mean Anomaly */
  if (O_A)
  { /* Linp is apparent longitude, so correct for aberration (light time) */
    r = a_earth*(1.0-Ex*Ex)/(1.0+Ex*cos(A));
    tab = r*vlight_d/JMIL;
  } else tab=0.0;
  corr = RMED( Mb - (Ma+tab*dM) , M_PI)/dM;
  t-=corr;
  /* alternatively (but less accurate) use Ma=Linp */

  /* Now use these values */
  for (i=5, Ex=0.0; i>0; i--) Ex = (Ex+E[i])*t ; Ex+=E[0];
  for (i=5, Ph=0.0; i>0; i--) Ph = (Ph+P[i])*t ; Ph+=P[0];
  for (i=5, dP=0.0; i>1; i--) dP = (dP+i*P[i])*t ; dP+=P[1];
  for (i=5, d2P=0.0; i>2; i--) d2P = (d2P+i*(i-1)*P[i])*t ; d2P+=2*P[2];
  for (i=5, Mb=0.0; i>0; i--) Mb = (Mb+M[i])*t ; Mb+=M[0];
  Mb = RMOD(Mb,M_PI*2);
  for (i=5, dM=0.0; i>1; i--) dM = (dM+i*M[i])*t ; dM+=M[1];
  for (i=5, d2M=0.0; i>2; i--) d2M = (d2M+i*(i-1)*M[i])*t ; d2M+=2*M[2];
  dL = dM + dP;
  d2L = d2M + d2P;
  C1 = (2 - (1.0/4.0)*Ex*Ex)*Ex;
  C2 = ((5.0/4.0) - (11.0/24.0)*Ex*Ex)*Ex*Ex;
  C3 = (13.0/12.0)*Ex*Ex*Ex;
  C4 = (103.0/96.0)*Ex*Ex*Ex*Ex;
  J = 2*M_PI/dL;   /* mean tropical year expressed in Julian milennia */
  Q = -C1*cos(Mb) -2*C2*cos(2*Mb) -3*C3*cos(3*Mb) -4*C4*cos(4*Mb)
      +C1*C1*cos(Mb)*cos(Mb) +4*C1*C2*cos(Mb)*cos(2*Mb) +6*C1*C3*cos(Mb)*cos(3*Mb) +4*C2*C2*cos(2*Mb)*cos(2*Mb)
     -C1*C1*C1*cos(Mb)*cos(Mb)*cos(Mb) -6*C1*C1*C2*C2*cos(Mb)*cos(Mb)*cos(2*Mb)
     +C1*C1*C1*C1*cos(Mb)*cos(Mb)*cos(Mb)*cos(Mb) ;
  DU = -J*dP*Q/dM ;

  Z = +C1*sin(Mb) +4*C2*sin(2*Mb)
      -2*C1*C1*cos(Mb)*sin(Mb) ;
  P1=P[1]; P2=P[2]; M1=M[1]; M2=M[2]; L1=M1+P1; L2=M2+P2;
  F1 = +P2*(1/M1)*Q;
  F2 = -P1*(1/M1)*(1/M1)*M2*Q;
  F3 = -P1*(1/M1)*(1/L1)*L2*Q;
  F4 = +P1*Z;

  if (DoDebug) fprintf(stderr,"D: dU components: %lg %lg %lg %lg\n",F1,F2,F3,F4);
  /* F1,F2,F3 are negligeable compared to F4 */
  dU = -J*(F1+F2+F3+F4); /* scaling? */
  dJ = -J*d2L/dL; /* scaling ? */
  if (DoDebug) fprintf(stderr,"D: dJ=%lg , dU=%lg\n",dJ,dU);

  jda = t*JMIL+J2000; /* not particularly good */
   *tya = (J + DU)*JMIL; /* tropical year expressed in days */
  *dtya = dJ + dU;

  if (DoDebug) fprintf(stderr,"D: end analytical\n");
  return(jda);
}


double numerical(double Linp, int O_A, double yinp, double *tyn, double *dtyn)
{
  int i,j;
  double y,t,r,tab;
  double corr,acc=0.001/DAY/JMIL;   /* 1ms/86400/365250 = 3E-14 */
  double jdarr[6];
  double Ph,Ma,Mb,Ea,A,Lo,Ex,dM,dL;

  if (DoDebug) fprintf(stderr,"D: begin numerical\n");

  for (j=0,y=yinp-1.0; j<3; y+=1.0, j++)
  { /* for target year, and years before and after */
    if (DoDebug) fprintf(stderr,"D: for j: j=%i, y=%lf\n",j,y);
    t = (y-2000.0)*0.001;   /* time since epoch J2000 in 1000's of Julian years */
    do
    {
      if (DoDebug) fprintf(stderr,"D: while\n");
      for (i=5, Ex=0.0; i>0; i--) Ex = (Ex+E[i])*t ; Ex+=E[0];
      for (i=5, Ph=0.0; i>0; i--) Ph = (Ph+P[i])*t ; Ph+=P[0];
      Ph += M_PI;   /* Want longitude of Sun, not Earth */
      /* Mean Anomaly for estimated time t : */
      for (i=5, Mb=0.0; i>0; i--) Mb = (Mb+M[i])*t ; Mb+=M[0];
      Mb = RMOD(Mb,M_PI*2);
      for (i=5, dM=0.0; i>1; i--) dM = (dM+i*M[i])*t ; dM+=M[1];

      /* Use Keplers equations forward */
       A = RMOD(Linp-Ph,M_PI*2);   /* estimated True Anomaly */
      Ea = 2*atan2( sin(A*0.5), cos(A*0.5)*sqrt((1+Ex)/(1-Ex)) );
      Ma = Ea - Ex*sin(Ea);   /* target Mean Anomaly */
      if (O_A)
      { /* Linp is apparent longitude, so correct for aberration (light time) */
        r = a_earth*(1.0-Ex*Ex)/(1.0+Ex*cos(A));
	tab = r*vlight_d/JMIL;
      } else tab=0.0;
      corr = RMED( Mb - (Ma+tab*dM) , M_PI)/dM;
      t-=corr;
      if (DoDebug) fprintf(stderr,"D: bottom of while: A=%lg, Ma=%lg, Mb=%lg, corr=%lg\n",A,Ma,Mb,corr);
    } while (fabs(corr)>acc);
    jdarr[j] = t*JMIL+J2000;
  }

  /* tabular differences */
  jdarr[3]=jdarr[1]-jdarr[0];
  jdarr[4]=jdarr[2]-jdarr[1];
  jdarr[5]=jdarr[4]-jdarr[3];
   *tyn = (jdarr[3]+jdarr[4])*0.5;
  *dtyn = jdarr[5];

  if (DoDebug) fprintf(stderr,"D: end numerical\n");
  return(jdarr[1]);
}


double mean(double Linp, int O_A, double yinp, double *tym, double *dtym)
{
  int i,j;
  double y,t,tab,jdm;
  double corr,acc=0.001/DAY/JMIL;   /* 1ms/86400/365250 = 3E-14 */
  double Ph,Ma,Mb,dP,dM,dL,d2P,d2M,d2L;

  if (DoDebug) fprintf(stderr,"D: begin mean\n");

    t = (yinp-2000.0)*0.001;   /* time since epoch J2000 in 1000's of Julian years */
    do
    {
      if (DoDebug) fprintf(stderr,"D: while\n");
      for (i=5, Ph=0.0; i>0; i--) Ph = (Ph+P[i])*t ; Ph+=P[0];
      Ph += M_PI;   /* Want longitude of Sun, not Earth */
      for (i=5, dP=0.0; i>1; i--) dP = (dP+i*P[i])*t ; dP+=P[1];
      for (i=5, d2P=0.0; i>2; i--) d2P = (d2P+i*(i-1)*P[i])*t ; d2P+=2*P[2];
      /* Mean Anomaly for estimated time t : */
      for (i=5, Mb=0.0; i>0; i--) Mb = (Mb+M[i])*t ; Mb+=M[0];
      Mb = RMOD(Mb,M_PI*2);
      for (i=5, dM=0.0; i>1; i--) dM = (dM+i*M[i])*t ; dM+=M[1];
      for (i=5, d2M=0.0; i>2; i--) d2M = (d2M+i*(i-1)*M[i])*t ; d2M+=2*M[2];
      dL = dM + dP;
      d2L = d2M + d2P;

      Ma = RMOD(Linp-Ph,M_PI*2);   /* estimated Mean Anomaly */
      if (O_A)
      { /* Linp is apparent mean longitude, so correct for aberration (light time) */
	tab = a_earth*vlight_d/JMIL;
      } else tab=0.0;
      corr = RMED( Mb - (Ma+tab*dM) , M_PI)/dM;
      t-=corr;
      if (DoDebug) fprintf(stderr,"D: bottom of while: Ma=%lg, Mb=%lg, corr=%lg\n",Ma,Mb,corr);
    } while (fabs(corr)>acc);
  jdm = t*JMIL+J2000;

   *tym = 2*M_PI*JMIL/dL;
  *dtym = -2*M_PI*JYR*d2L/(dL*dL);

  if (DoDebug) fprintf(stderr,"D: end mean\n");
  return(jdm);
}


void usage()
{
  if (DoDebug) fprintf(stderr,"D: begin usage\n");
  printf(
"Usage:\n"
"tropical_year [-AMNFEUCJGS] <longitude> <startyear> [numyear [stepsize [calyear]]]\n"
"  where:\n"
"\t-\tany combination of 1-letter options:\n"
"\t  A\tApparent longitude (incl. aberration = light time correction)\n"
"\t  M\tMean longitude, mean tropical year\n"
"\t  N\tNumerical approximation for elliptic orbit (default)\n"
"\t  F\tAnalytical approximation with equation of center formula\n"
"\t  E\tresults in Ephemeris Time scale (default)\n" 
"\t  U\tresults in Universal Time scale\n"
"\t  C|J|G\tresults with calendar date: conventional, Julian, Gregorian\n"
"\t  S\ttime results in sexagesimal notation (hh:mm:ss)\n"
"\tlongitude\tlongitude of the ecliptic reference point in degrees\n"
"\tstartyear\tfractional julian year (we use J2000 as epoch)\n"
"\tnumyear \tnumber of years for output (optional; default=1)\n"
"\tstepsize\tstep in years for each increment (opt.; def.=1.0)\n"
"\tcalyear \tcalendar year as reference length (opt.; def.=365.2425)\n"
  );
  printf(
"Output:\n"
"mean ET: JD of mean Sun at specified longitude; ;\n"
"\tmean TY; ; annual change of TY; mean-calyr\n"
"num. ET: JD of true Sun at specified longitude; true-mean JD;\n"
"\tTY for spec. longitude; true-mean TY; annual change of TY; true-calyr\n"
"ana. ET: JD from analytical expression; ana-num JD;\n"
"\tanalytical TY; ana-num TY; annual change of TY; ana-calyr\n"
"mean UT: JD of mean Sun at specified longitude; UT-ET;\n"
"\tmean TY; TY(UT)-TY(ET); annual change of TY; mean-cal\n"
"num. UT: JD of true Sun at specified longitude; true-mean JD;\n"
"\tTY for spec. longitude; true-mean TY; annual change of TY; true-calyr\n"
"ana. UT: JD from analytical expression; ana-mean JD;\n"
"\tanalytical TY; ana-num TY; annual change of TY; ana-calyr\n"
  );
  if (DoDebug) fprintf(stderr,"D: end usage\n");
}


void output(char outS[], int O_S, char Ctype, char *yrtyp, char *tmtyp, double jd, double jdif, double ty, double tydif, double dty, double caldif)
{
  int yi,mi,di;
  double d;
  SEXAGES hms;
  char C,sgnC;
  char LeadS[30],jdS[17],dateS[21],timS[13],
       jdifS[12],tyS[13],tydifS[12],dtyS[16],caldifS[10];

  C=Ctype;
  if ( (C=='C')||(C=='c') ||
       (C=='J')||(C=='j') ||
       (C=='G')||(C=='g') )
  { /* print calendar date */
    datjul(jd,&C,&yi,&mi,&d);
    if (O_S)
    {
      di=INT(d); d-=di;
      (void)R_SEXA(d*24.0,&hms);
      sprintf(timS, "%02i:%02i:%02i", hms.h, hms.m, INT(ROUNDR(hms.s)) );
      sprintf(dateS,"%c%5i/%02i/%02i;%s",C,yi,mi,di,timS);
      sprintf(LeadS,"%.4s: %s%.2s",yrtyp,dateS,tmtyp);
    } else {
      sprintf(dateS,"%c%5i/%02i/%08.5lf",C,yi,mi,d);
      sprintf(LeadS,"%.4s %.2s: %s",yrtyp,tmtyp,dateS);
    }
  } else { /* print Julian day */
    if (O_S)
    {
      di=INT(jd+0.5); /* day number */
       d=jd+0.5-di;   /* fraction of day since midnight */
      (void)R_SEXA(d*24.0,&hms);
      sprintf(timS, "%02i:%02i:%02i", hms.h, hms.m, INT(ROUNDR(hms.s)) );
      sprintf(LeadS,"%.4s: %7i;%s%.2s",yrtyp,di,timS,tmtyp);
    } else {
      sprintf(LeadS,"%.4s %.2s: %13.5lf",yrtyp,tmtyp,jd);
    }
  }

  if (O_S)
  {
    di=TRUNC(jdif);
    d=jdif-di;
    (void)R_SEXA(fabs(d*24.0),&hms);
    sprintf(timS, "%02i:%02i:%02i", hms.h, hms.m, INT(ROUNDR(hms.s)) );
    sprintf(jdifS,"%+2id%s",di,timS);
  } else {
    sprintf(jdifS,"%+8.5lf",jdif);
  }

  if (O_S)
  {
    di=TRUNC(ty);
    d=ty-di;
    (void)R_SEXA(d*24.0,&hms);
    sprintf(timS, "%02i:%02i:%02i", hms.h, hms.m, INT(ROUNDR(hms.s)) );
    sprintf(tyS,"%3id%s",di,timS);
  } else {
    sprintf(tyS,"%11.7lf",ty);
  }

  if (O_S)
  {
    (void)R_SEXA(tydif*24.0,&hms);
    sgnC=(hms.t<0?'-':' ');
    sprintf(tydifS, "%c%02i:%02i:%02i", sgnC, hms.h, hms.m, INT(ROUNDR(hms.s)) );
  } else {
    sprintf(tydifS,"%+10.7lf",tydif);
  }

  if (O_S)
  {
    (void)R_SEXA(dty*24.0,&hms);
    sgnC=(hms.t<0?'-':' ');
    sprintf(dtyS, "%c%1i:%02i:%06.3f", sgnC, hms.h, hms.m, hms.s );
  } else {
    sprintf(dtyS,"%+12.9lf",dty);
  }

  if (O_S)
  {
    (void)R_SEXA(caldif*24.0,&hms);
    sgnC=(hms.t<0?'-':' ');
    sprintf(caldifS, "%c%1i:%02i:%02i", sgnC, hms.h, hms.m, INT(ROUNDR(hms.s)) );
  } else {
    sprintf(caldifS,"%+10.7lf",caldif);
  }

  sprintf(outS,"%s %s %s %s %s %s",LeadS,jdifS,tyS,tydifS,dtyS,caldifS);
}


int main(int argc, char *argv[])
{
  int O_A=0,O_M=0,O_N=0,O_F=0,O_E=0,O_U=0,O_C=0,O_J=0,O_G=0,O_S=0,Off=0;
  char *pc,Ctype='\0',optS[9]="",outS[99]="";
  double linp,ystart,yinp,ystep,calyr;
  double jdm,tym,dtym,jda,tya,dtya,jdn,tyn,dtyn;
  double tj19;
  double delT,delD;
  int Nyr,i;
  int yi,mi;
  double d,di;
  SEXAGES hms;

  DoDebug = getenv("DEBUG");
  if (DoDebug) fprintf(stderr,"D: begin main\n");
  
  if ( (argc<3) || (argc>7) )
  { 
    usage();
    exit(argc);
  }

  
  pc = argv[1];
  if (*pc=='-'||isupper((int)*pc))
  { /* output options */
    Off=1;
    while (*pc!=0)
    {
      if (*pc=='A') { O_A=1; strcat(optS,"A"); }
      if (*pc=='M') { O_M=1; strcat(optS,"M"); }
      if (*pc=='N') { O_N=1; strcat(optS,"N"); }
      if (*pc=='F') { O_F=1; strcat(optS,"F"); }
      if (*pc=='E') { O_E=1; strcat(optS,"E"); }
      if (*pc=='U') { O_U=1; strcat(optS,"U"); }
      if (*pc=='C') { O_C=1; strcat(optS,"C"); Ctype=*pc;}
      if (*pc=='J') { O_J=1; strcat(optS,"J"); Ctype=*pc;}
      if (*pc=='G') { O_G=1; strcat(optS,"G"); Ctype=*pc;}
      if (*pc=='S') { O_S=1; strcat(optS,"S"); }
      pc++;
    }
  } else { /* defaults */
    Off=0;
    O_A=0;
    O_M=0;
    O_N=1;
    O_F=0;
    O_E=1;
    O_U=0;
    O_C=0;
    O_J=0;
    O_G=0;
    O_S=0;
  }
  /* Sanity check */
  if (O_C+O_J+O_G>1)
  {
    fprintf(stderr,"E: specify only 1 type of output date ('C' or 'J' or 'G'\n");
    exit (-1);
  }
  if ( (O_E==0) && (O_U==0) ) O_E=1;
  if ( (O_M==0) && (O_N==0) && (O_F==0) ) O_N=1;


  linp = atof(argv[1+Off])*M_PI/180.0;
  ystart = atof(argv[2+Off]);
  if (argc>=(4+Off)) Nyr = atoi(argv[3+Off]);
  else Nyr = 1;
  if (argc>=(5+Off)) ystep = atof(argv[4+Off]);
  else ystep = 1.0;
  if (argc==(6+Off)) calyr = atof(argv[5+Off]);
  else calyr = GYR;
  if (DoDebug) 
    fprintf(stderr,"D: input pars: opts=\"%s\", linp=%lg, ystart=%lg, Nyr=%i, ystep=%lg, calyr=%lg\n",optS,linp,ystart,Nyr,ystep,calyr);


  /* Print output headers */

  printf(
"Input: Options=\"%s\"  Longitude=%lg  Start year=%lg\n"
"\tNumber=%i  Step=%lg  Reference Calendar year=%11.7lf\n",
         optS,linp,ystart,Nyr,ystep,calyr
        );
  printf("\n");

  if (O_C||O_J||O_G)
    if (O_S)
      printf(
"Output: date of event          diff       year length   diff     annual change  -calyr\n"
            );
    else
      printf(
"Output:  date of event       diff    year length  diff      annual change  -calyr  \n"
            );
  else
    if (O_S)
      printf(
"Output:    JD of event    diff       year length   diff     annual change  -calyr\n"
            );
    else
      printf(
"Output:  JD of event    diff    year length  diff      annual change  -calyr  \n"
            );


  for (i=0, yinp=ystart; i<Nyr; yinp+=ystep, i++)
  {
    jdm = mean(linp,O_A,yinp,&tym,&dtym);
    jdn = numerical(linp,O_A,yinp,&tyn,&dtyn);
  /* Not as accurate, and annual change in error */
    jda = analytical(linp,O_A,yinp,&tya,&dtya);

    /* Assume Delta-T=0 and mean solar day equal to 86400 SI seconds at J1900 */
    tj19 = (jdn-J1900)/JYR;
    delT = 36E-9*tj19*tj19; /* estimated [ephemeris time - universal time] */
    delD = 1.0/(1.0+2.0E-10*tj19);
  /* estimated length of ephemeris day = SI day, expressed in mean solar days */

    /* Output of results */
    printf("\n");
    if (O_E)
    {
      if (O_M) {
        output(outS,O_S,Ctype,"mean","ET",jdm,0.0,tym,0.0,dtym,tym-calyr);
        puts(outS);
      }
      if (O_N) {
        output(outS,O_S,Ctype,"num.","ET",jdn,jdn-jdm,tyn,tyn-tym,dtyn,tyn-calyr);
        puts(outS);
      }
      if (O_F) {
        output(outS,O_S,Ctype,"ana.","ET",jda,jda-jdn,tya,tya-tyn,dtya,tya-calyr);
        puts(outS);
      }
    }
    if (O_U)
    {
      if (O_M) {
        output(outS,O_S,Ctype,"mean","UT",jdm-delT,-delT,tym*delD,tym*(delD-1.0),dtym*delD,tym*delD-calyr);
        puts(outS);
      }
      if (O_N) {
output(outS,O_S,Ctype,"num.","UT",jdn-delT,(jdn-jdm)*delD,tyn*delD,(tyn-tym)*delD,dtyn*delD,tyn*delD-calyr);
        puts(outS);
      }
      if (O_F) {
        output(outS,O_S,Ctype,"ana.","UT",jda-delT,(jda-jdn)*delD,tya*delD,(tya-tyn)*delD,dtya*delD,tya*delD-calyr);
        puts(outS);
      }
      /* notes:
         - the crossterm for the transformation t(ET)->t(UT) 
           is incorporated into the year length
         - deltaT is in ET time scale, difference with UT time scale 
           is negligeable
         - the differences (jdn-jdm) and (jda-jdn) for ET and UT are 
           expressed in SI seconds and current mean solar day seconds 
           respectively; the difference is of course negligeable.
      */
    }

  } /* next i */


  if (DoDebug) fprintf(stderr,"D: end main\n");
  return(0);
}