wmmoonclock/src/CalcEphem.c

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