/*
################################################  
 SOLAR EPHEMERIS CALCULATION TO GET
               P, B0 AND L0, SEMI DIAMETER
               RA and decl (Epoch is given date)

     Algorithm from:  Duffett-Smith, 1988 
 
            by: Anders Johannesson 1997-10-08

     Mod. Oct. 25 AJO  Equatorial Coordinates  
################################################ 
*/ 

#include "soleph.h"
  
void soleph(TimeType Etime, EphType *Ephemeris_p)
{

/* DECLARATIONS */ 
short year, month, day, hour, minute, second;
short ye, mo;
long A, B, C, D;
double d_in, jd, DD, N, M, eccentricity, delta, E, x, y;
double de, tanvo2, vrad, v, longitude, latitude, S, T, bdelta, omega;
double AA, Mprim, MM, L0, B0, P1, P2, P, dE, epsilon;
double decl, RA, T0, jd_ut0, ut, gst, lst, H, Si;

/* TIME */
year=Etime.iy;
month=Etime.im;
day=Etime.id;
hour=Etime.ih;
minute=Etime.imi;
second=Etime.is;

/* JULIAN DATE*/ 
d_in=day+hour/24.0+minute/(24.0*60.0)+second/(24.0*3600.0); 
if ((month == 1) | (month == 2)) 
{ 
   ye=year-1; 
   mo=month+12; 
}
else
{ 
   ye=year; 
   mo=month; 
}
if (ye > 1582)
{ 
   A=(long)(ye/100.0); 
   B=2-A+(long)(A/4.0); 
}
else
{ 
   B=0; 
} 
C=(long)(365.25*ye); 
D=(long)(30.6001*(mo+1)); 
jd=B+C+D+d_in+1720994.5; 
 
/* THE NUMBER OF DAYS SINCE THE SELECTED REFERENCE EPOCH */
DD=jd-2447891.5; 
 
/* CALCULATE MEAN ANOMALY */ 
N=360.0/365.242191*DD; 
N=fmod(N,360.0); 
if (N < 0.0)
   N=N+360.0; 
M=N+279.403303 -282.768422; 
if (M < 0) 
   M=M+360.0; 
M=M/180.0*DPI; 
 
/* SOLVE KEPLERS EQUATION */
eccentricity=0.016713; 
delta=1.0; 
E=M; 
while (fabs(delta) > 1.0E-6)
{ 
   delta=E-eccentricity*sin(E)-M; 
   dE=delta/(1.0 -eccentricity*cos(E)); 
   E=E-dE; 
} 
 
/* FIND TRUE ANOMALY */
tanvo2=sqrt((1.0 +eccentricity)/(1.0 -eccentricity))*tan(E/2.0); 
vrad=atan(tanvo2)*2.0; 
v=vrad*180.0/DPI; 
 
/* FIND THE SUN'S ECLIPTIC LONGITUDE */
longitude=v+282.768422; 
longitude= fmod(longitude, 360.0); 
if (longitude < 0.0) 
   longitude=longitude+360.0; 

/* FIND SEMI DIAMETER */
S=(1.0 +eccentricity*cos(vrad))/(1-eccentricity*eccentricity); 
S=S*0.533128*3600.0 /2.0; 
 
/* HELIOGRAPHIC COORDINATES */
T=(jd-2415020.0)/36525.0;  /* centuries since reference Epoch */ 
bdelta=84.0*T/60.0; 
omega=74.0 +22.0/60.0 +bdelta; 
y=sin((omega-longitude)/180.0*DPI)*cos(7.25/180.0*DPI); 
x=-cos((omega-longitude)/180.0*DPI); 
AA=atan2(y,x)*180.0/DPI; 
Mprim=360.0/25.38*(jd-2398220.0); 
Mprim= fmod(Mprim, 360.0);
if (Mprim < 0.0) 
   Mprim=Mprim+360.0; 
MM=360.0 -Mprim; 
L0=MM+AA; 
L0=fmod(L0,360.0); 
if (L0 < 0.0) 
   L0=L0+360.0;  
B0=asin(sin((longitude-omega)/180.0*DPI)*sin(7.25/180.0*DPI)); 
B0=B0*180.0/DPI; 
P1=atan(-cos(longitude/180.0*DPI)*tan(23.441884 /180.0*DPI)); 
P2=atan(-cos((omega-longitude)/180.0*DPI)*tan(7.25 /180.0*DPI)); 
P=(P1+P2)*180.0/DPI;

/* EQUATORIAL COORDINATES (Epoch is todays date) */
latitude=0.0; /* radians */
T=(jd-2451545.0)/36525.0;  /* centuries since reference Epoch */ 
epsilon=(23.439292-(46.815*T+0.0006*T*T-0.00181*T*T*T)/3600.0)/180.0*DPI;
decl=asin(sin(latitude)*cos(epsilon)+cos(latitude)*sin(epsilon)*sin(longitude/180.0*DPI));
decl=decl*180.0/DPI;
y=sin(longitude/180.0*DPI)*cos(epsilon)-tan(latitude)*sin(epsilon);
x=cos(longitude/180.0*DPI);
RA=atan2(y,x)*180.0/DPI/15.0;
if (RA < 0) {
   RA = RA + 24.0;
} 

/* CALCULATE HOUR ANGLE (USES SOME OF THE CALCULATIONS ABOVE) */
 
/* JULIAN DATE FOR 0h UT */ 
d_in=day+0.0; 
jd_ut0=B+C+D+d_in+1720994.5;

/* GREENWICH MEAN SIDEREAL TIME */
Si=jd_ut0-2451545.0;
T=Si/36525.0;
T0=6.697374558+(2400.051336*T)+(0.000025862*T*T);
T0=fmod(T0,24.0); 
if (T0 < 0.0) 
   T0=T0+24.0;
ut=hour+minute/60.0+second/3600.0;
gst=1.002737909*ut+T0;
gst=fmod(gst,24.0); 
if (gst < 0.0) 
   gst=gst+24.0;

/* LOCAL SIDEREAL TIME */
lst=gst-7.79433335;
if (lst > 24.0) 
   lst=lst-24.0;
if (lst < 0)
   lst=lst+24.0;

/* HOUR ANGLE */
H=lst-RA;
if (H < 0)
   H=H+24.0;

/* FILL EPHEMERIS ARRAY */
Ephemeris_p->rad=S;
Ephemeris_p->p=P;
Ephemeris_p->b=B0;
Ephemeris_p->l0=L0;
Ephemeris_p->decl=decl;
Ephemeris_p->RA=RA;
Ephemeris_p->H=H;

}