wmmoonclock/src/CalcEphem.c

   1 #include "CalcEphem.h"
   2 #include <math.h>
   3 #include <string.h>
   4
   5 void CalcEphem(date, UT, c)
   6 long int         date;  /* integer containing the date (e.g. 960829) */
   7 double           UT;    /* Universal Time */
   8 CTrans          *c;     /* structure containing all the relevent coord trans info */
   9 {
  10
  11     int    year, month, day;
  12     double TU, TU2, TU3, T0;
  13     double varep, varpi;
  14     double eccen, epsilon;
  15     double days, M, E, nu, lambnew;
  16     double r0, earth_sun_distance;
  17     double RA, DEC, RA_Moon, DEC_Moon;
  18     double TDT, AGE, LambdaMoon, BetaMoon, R;
  19     double jd(), hour24(), angle2pi(), angle360(), kepler(), Moon(), NewMoon();
  20     double Ta, Tb, Tc, frac();
  21     double SinGlat, CosGlat, Tau, lmst, x, y, z;
  22     double SinTau, CosTau, SinDec, CosDec;
  23
  24
  25
  26
  27
  28     c->UT = UT;
  29     year = (int)(date/10000);
  30     month = (int)( (date - year*10000)/100 );
  31     day = (int)( date - year*10000 - month*100 );
  32     c->year = year;
  33     c->month = month;
  34     c->day = day;
  35
  36     c->doy = DayofYear(year, month, day);
  37     c->dow = DayofWeek(year, month, day, c->dowstr);
  38
  39
  40
  41     /*
  42      *  Compute Greenwich Mean Sidereal Time (gmst)
  43      *  The TU here is number of Julian centuries
  44      *  since 2000 January 1.5
  45      *  From the 1996 astronomical almanac
  46      */
  47     TU = (jd(year, month, day, 0.0) - 2451545.0)/36525.0;
  48     TU2 = TU*TU;
  49     TU3 = TU2*TU;
  50     T0 = (6.0 + 41.0/60.0 + 50.54841/3600.0) + 8640184.812866/3600.0*TU
  51             + 0.093104/3600.0*TU2 - 6.2e-6/3600.0*TU3;
  52     T0 = hour24(T0);
  53     c->gmst = hour24(T0 + UT*1.002737909);
  54
  55     lmst = 24.0*frac( (c->gmst - c->Glon/15.0) / 24.0 );
  56
  57
  58
  59
  60
  61
  62     /*
  63      *
  64      *   Construct Transformation Matrix from GEI to GSE  systems
  65      *
  66      *
  67      *   First compute:
  68      *          mean ecliptic longitude of sun at epoch TU (varep)
  69      *          elciptic longitude of perigee at epoch TU (varpi)
  70      *          eccentricity of orbit at epoch TU (eccen)
  71      *
  72      *   The TU here is the number of Julian centuries since
  73      *   1900 January 0.0 (= 2415020.0)
  74      */
  75     TDT = UT + 59.0/3600.0;
  76     TU = (jd(year, month, day, TDT) - 2415020.0)/36525.0;
  77     varep = (279.6966778 + 36000.76892*TU + 0.0003025*TU*TU)*RadPerDeg;
  78     varpi = (281.2208444 + 1.719175*TU + 0.000452778*TU*TU)*RadPerDeg;
  79     eccen = 0.01675104 - 0.0000418*TU - 0.000000126*TU*TU;
  80     c->eccentricity = eccen;
  81
  82
  83
  84     /*
  85      * Compute the Obliquity of the Ecliptic at epoch TU
  86      * The TU in this formula is the number of Julian
  87      * centuries since epoch 2000 January 1.5
  88      */
  89     TU  = (jd(year, month, day, TDT) - jd(2000, 1, 1, 12.0))/36525.0;
  90     epsilon = (23.43929167 - 0.013004166*TU - 1.6666667e-7*TU*TU
  91                 - 5.0277777778e-7*TU*TU*TU)*RadPerDeg;
  92     c->epsilon = epsilon;
  93
  94
  95     /*
  96      * Compute:
  97      *          Number of Days since epoch 1990.0 (days)
  98      *          The Mean Anomaly (M)
  99      *          The True Anomaly (nu)
 100      *      The Eccentric Anomaly via Keplers equation (E)
 101      *
 102      *
 103      */
 104     days  = jd(year, month, day, TDT) - jd(year, month, day, TDT);
 105     M = angle2pi(2.0*M_PI/365.242191*days + varep - varpi);
 106     E = kepler(M, eccen);
 107     nu = 2.0*atan( sqrt((1.0+eccen)/(1.0-eccen))*tan(E/2.0) );
 108     lambnew = angle2pi(nu + varpi);
 109     c->lambda_sun = lambnew;
 110
 111
 112     /*
 113      *  Compute distance from earth to the sun
 114      */
 115     r0 = 1.495985e8;  /* in km */
 116     earth_sun_distance  = r0*(1-eccen*eccen)/(1.0 + eccen*cos(nu))/6371.2;
 117     c->earth_sun_dist = earth_sun_distance;
 118
 119
 120
 121
 122
 123     /*
 124      * Compute Right Ascension and Declination of the Sun
 125      */
 126     RA = angle360(atan2(sin(lambnew)*cos(epsilon), cos(lambnew))*180.0/M_PI);
 127     DEC = asin(sin(epsilon)*sin(lambnew))*180.0/M_PI;
 128     c->RA_sun = RA;
 129     c->DEC_sun = DEC;
 130
 131
 132
 133
 134
 135
 136     /*
 137      * Compute Moon Phase and AGE Stuff. The AGE that comes out of Moon()
 138      * is actually the Phase converted to days. Since AGE is actually defined
 139      * to be time since last NewMoon, we need to figure out what the JD of the
 140      * last new moon was. Thats done below....
 141      */
 142     TU = (jd(year, month, day, TDT) - 2451545.0)/36525.0;
 143     c->MoonPhase = Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE);
 144     LambdaMoon *= RadPerDeg;
 145     BetaMoon *= RadPerDeg;
 146
 147
 148     RA_Moon  = angle360(atan2(sin(LambdaMoon)*cos(epsilon)-tan(BetaMoon)*sin(epsilon), cos(LambdaMoon))*DegPerRad);
 149     DEC_Moon = asin( sin(BetaMoon)*cos(epsilon) + cos(BetaMoon)*sin(epsilon)*sin(LambdaMoon))*DegPerRad;
 150     c->RA_moon = RA_Moon;
 151     c->DEC_moon = DEC_Moon;
 152
 153
 154     /*
 155      *  Compute Alt/Az coords
 156      */
 157     Tau = (15.0*lmst - RA_Moon)*RadPerDeg;
 158     CosGlat = cos(c->Glat*RadPerDeg); SinGlat = sin(c->Glat*RadPerDeg);
 159     CosTau = cos(Tau); SinTau = sin(Tau);
 160     SinDec = sin(DEC_Moon*RadPerDeg); CosDec = cos(DEC_Moon*RadPerDeg);
 161     x = CosDec*CosTau*SinGlat - SinDec*CosGlat;
 162     y = CosDec*SinTau;
 163     z = CosDec*CosTau*CosGlat + SinDec*SinGlat;
 164     c->A_moon = DegPerRad*atan2(y, x);
 165     c->h_moon = DegPerRad*asin(z);
 166     c->Visible = (c->h_moon < 0.0) ? 0 : 1;
 167
 168
 169
 170     /*
 171      * Compute accurate AGE of the Moon
 172      */
 173     Tb = TU - AGE/36525.0; /* should be very close to minimum */
 174     Ta = Tb - 0.4/36525.0;
 175     Tc = Tb + 0.4/36525.0;
 176     c->MoonAge = (TU - NewMoon(Ta, Tb, Tc))*36525.0;
 177
 178
 179
 180     /*
 181      * Compute Earth-Moon distance
 182      */
 183     c->EarthMoonDistance = R;
 184
 185
 186
 187 }
 188
 189
 190
 191
 192
 193 double kepler(M, e)
 194 double M, e;
 195 {
 196         int n=0;
 197         double E, Eold, eps = 1.0e-8;
 198
 199
 200
 201         E = M + e*sin(M);
 202         do{
 203                 Eold = E;
 204                 E = Eold + (M-Eold+e*sin(Eold))
 205                         /(1.0-e*cos(Eold));
 206                 ++n;
 207         }while((fabs(E-Eold) > eps) && (n < 100));
 208         return(E);
 209 }
 210
 211
 212
 213
 214 int DayofYear(year, month, day)
 215 int year, month, day;
 216 {
 217         double jd();
 218         return((int)(jd(year, month, day, 0.0) - jd(year, 1, 0, 0.0)));
 219 }
 220
 221
 222
 223
 224 int DayofWeek(year, month, day, dowstr)
 225 int year, month, day;
 226 char dowstr[];
 227 {
 228         double JD, A, Afrac, jd();
 229         int n, iA;
 230
 231         JD = jd(year, month, day, 0.0);
 232         A = (JD + 1.5)/7.0;
 233         iA = (int)A;
 234         Afrac = A - (double)iA;
 235         n = (int)(Afrac*7.0 + 0.5);
 236         switch(n){
 237                 case 0:
 238                         strcpy(dowstr, "Sunday");
 239                         break;
 240                 case 1:
 241                         strcpy(dowstr, "Monday");
 242                         break;
 243                 case 2:
 244                         strcpy(dowstr, "Tuesday");
 245                         break;
 246                 case 3:
 247                         strcpy(dowstr, "Wednesday");
 248                         break;
 249                 case 4:
 250                         strcpy(dowstr, "Thursday");
 251                         break;
 252                 case 5:
 253                         strcpy(dowstr, "Friday");
 254                         break;
 255                 case 6:
 256                         strcpy(dowstr, "Saturday");
 257                         break;
 258         }
 259         return(n);
 260 }
 261
 262
 263
 264
 265
 266 /*
 267  *  Compute the Julian Day number for the given date.
 268  *  Julian Date is the number of days since noon of Jan 1 4713 B.C.
 269  */
 270 double jd(ny, nm, nd, UT)
 271 int ny, nm, nd;
 272 double UT;
 273 {
 274         double A, B, C, D, JD, day;
 275
 276         day = nd + UT/24.0;
 277
 278
 279         if ((nm == 1) || (nm == 2)){
 280                 ny = ny - 1;
 281                 nm = nm + 12;
 282         }
 283
 284         if (((double)ny+nm/12.0+day/365.25)>=(1582.0+10.0/12.0+15.0/365.25)){
 285                         A = ((int)(ny / 100.0));
 286                         B = 2.0 - A + (int)(A/4.0);
 287         }
 288         else{
 289                         B = 0.0;
 290         }
 291
 292         if (ny < 0.0){
 293                 C = (int)((365.25*(double)ny) - 0.75);
 294         }
 295         else{
 296                 C = (int)(365.25*(double)ny);
 297         }
 298
 299         D = (int)(30.6001*(double)(nm+1));
 300
 301
 302         JD = B + C + D + day + 1720994.5;
 303         return(JD);
 304
 305 }
 306
 307 double hour24(hour)
 308 double hour;
 309 {
 310         int n;
 311
 312         if (hour < 0.0){
 313                 n = (int)(hour/24.0) - 1;
 314                 return(hour-n*24.0);
 315         }
 316         else if (hour > 24.0){
 317                 n = (int)(hour/24.0);
 318                 return(hour-n*24.0);
 319         }
 320         else{
 321                 return(hour);
 322         }
 323 }
 324
 325 double angle2pi(angle)
 326 double angle;
 327 {
 328         int n;
 329         double a;
 330         a = 2.0*M_PI;
 331
 332         if (angle < 0.0){
 333                 n = (int)(angle/a) - 1;
 334                 return(angle-n*a);
 335         }
 336         else if (angle > a){
 337                 n = (int)(angle/a);
 338                 return(angle-n*a);
 339         }
 340         else{
 341                 return(angle);
 342         }
 343 }
 344
 345 double angle360(angle)
 346 double angle;
 347 {
 348         int n;
 349
 350         if (angle < 0.0){
 351                 n = (int)(angle/360.0) - 1;
 352                 return(angle-n*360.0);
 353         }
 354         else if (angle > 360.0){
 355                 n = (int)(angle/360.0);
 356                 return(angle-n*360.0);
 357         }
 358         else{
 359                 return(angle);
 360         }
 361 }
 362
 363
 364 void Radec_to_Cart(ra, dec, r)
 365 double  ra, dec;        /* RA and DEC */
 366 Vector *r;              /* returns corresponding cartesian unit vector */
 367 {
 368
 369     /*
 370      *  Convert ra/dec from degrees to radians
 371      */
 372     ra  *= RadPerDeg;
 373     dec *= RadPerDeg;
 374
 375
 376     /*
 377      *  Compute cartesian coordinates (in GEI)
 378      */
 379     r->x = cos(dec) * cos(ra);
 380     r->y = cos(dec) * sin(ra);
 381     r->z = sin(dec);
 382
 383 }
 384
 385
 386
 387
 388
 389 int LeapYear(year)
 390 int year;
 391 {
 392     if ((year%100 == 0)&&(year%400 != 0)) return(0);
 393     else if (year%4 == 0) return(1);
 394     else return(0);
 395 }
 396
 397
 398
 399