/* Expansions for the geocentric ecliptic longitude, latitude, * and equatorial horizontal parallax of the Moon referred to * the mean equinox and ecliptic of date. * * Coefficients are derived from * M. Chapront-Touze' and J. Chapront, "ELP2000-85: a semi-analytical * lunar ephemeris adequate for historical times," Astronomy and * Astrophysics 190, 342-352 (1988). * In this program, their main problem coefficients have been rounded * off to the nearest .0001". * The authors give the precision as about 20" from 2000 A.D. to 1500 B.C. * I have spot checked this program against the DE102 at 400-day * intervals with the following results for the maximum errors in 500-year * intervals: * * Years * From To Long Lat Radius * -1410 -1000 33.3" 6.4" 14.0 x10^-8 au * -1000 -500 29.4 7.4 10.4 * -500 0 20.2 7.1 9.0 * 0 500 15.3 5.7 6.9 * 500 1000 8.4 4.9 7.7 * 1000 1500 7.6 4.5 5.2 * 1500 2000 6.0 5.7 5.6 * 2000 2500 6.7 4.9 5.6 * 2500 3002 14.4 5.4 6.8 * * For this comparison the DE102 was corrected by subtracting * 0.024" per century squared from the longitude. It should * be noted that the real longitude error due to uncertainty * in the tidal acceleration of the Moon grows by at least 0.1" per * century squared. This means that the real longitude of the * Moon in ancient times is uncertain by an amount that is * considerably greater than the theoretical error of the program. * * The program has two subroutine entry points. * Entry gmoon() returns the geometric position of the Moon * relative to the Earth. Its calling procedure is as follows: * * double JD; input Julian Ephemeris Date * double rect[3]; output equatorial x, y, z, in astronomical units * double pol[3]; output ecliptic polar coordinatees in radians and au * gmoon( JD, rect, pol ); * * The second entry, domoon(), is intended to be called by the AA.ARC * ephemeris calculator. It calculates the Moon's apparent position. * Corrections for nutation and light time are included, but the * required routines for this are given as separate files in AA.ARC. * Note that the difference between Ephemeris and Universal * time is significant here, since the Moon moves in right * ascension at the rate of about 2s/minute. * * * - S. L. Moshier, August, 1991 */ #include "kep.h" /* Perturbation tables */ #define NLR 73 static short LR[8*NLR] = { /* Longitude Radius D l' l F 1" .0001" 1km .0001km */ 0, 0, 1, 0, 22639, 5858,-20905,-3550, 2, 0,-1, 0, 4586, 4383, -3699,-1109, 2, 0, 0, 0, 2369, 9139, -2955,-9676, 0, 0, 2, 0, 769, 257, -569,-9251, 0, 1, 0, 0, -666,-4171, 48, 8883, 0, 0, 0, 2, -411,-5957, -3,-1483, 2, 0,-2, 0, 211, 6556, 246, 1585, 2,-1,-1, 0, 205, 4358, -152,-1377, 2, 0, 1, 0, 191, 9562, -170,-7331, 2,-1, 0, 0, 164, 7285, -204,-5860, 0, 1,-1, 0, -147,-3213, -129,-6201, 1, 0, 0, 0, -124,-9881, 108, 7427, 0, 1, 1, 0, -109,-3803, 104, 7552, 2, 0, 0,-2, 55, 1771, 10, 3211, 0, 0, 1, 2, -45, -996, 0, 0, 0, 0, 1,-2, 39, 5333, 79, 6606, 4, 0,-1, 0, 38, 4298, -34,-7825, 0, 0, 3, 0, 36, 1238, -23,-2104, 4, 0,-2, 0, 30, 7726, -21,-6363, 2, 1,-1, 0, -28,-3971, 24, 2085, 2, 1, 0, 0, -24,-3582, 30, 8238, 1, 0,-1, 0, -18,-5847, -8,-3791, 1, 1, 0, 0, 17, 9545, -16,-6747, 2,-1, 1, 0, 14, 5303, -12,-8314, 2, 0, 2, 0, 14, 3797, -10,-4448, 4, 0, 0, 0, 13, 8991, -11,-6500, 2, 0,-3, 0, 13, 1941, 14, 4027, 0, 1,-2, 0, -9,-6791, -7, -27, 2, 0,-1, 2, -9,-3659, 0, 0, 2,-1,-2, 0, 8, 6055, 10, 562, 1, 0, 1, 0, -8,-4531, 6, 3220, 2,-2, 0, 0, 8, 502, -9,-8845, 0, 1, 2, 0, -7,-6302, 5, 7509, 0, 2, 0, 0, -7,-4475, 1, 657, 2,-2,-1, 0, 7, 3712, -4,-9501, 2, 0, 1,-2, -6,-3832, 4, 1311, 2, 0, 0, 2, -5,-7416, 0, 0, 4,-1,-1, 0, 4, 3740, -3,-9580, 0, 0, 2, 2, -3,-9976, 0, 0, 3, 0,-1, 0, -3,-2097, 3, 2582, 2, 1, 1, 0, -2,-9145, 2, 6164, 4,-1,-2, 0, 2, 7319, -1,-8970, 0, 2,-1, 0, -2,-5679, -2,-1171, 2, 2,-1, 0, -2,-5212, 2, 3536, 2, 1,-2, 0, 2, 4889, 0, 0, 2,-1, 0,-2, 2, 1461, 0, 0, 4, 0, 1, 0, 1, 9777, -1,-4226, 0, 0, 4, 0, 1, 9337, -1,-1169, 4,-1, 0, 0, 1, 8708, -1,-5714, 1, 0,-2, 0, -1,-7530, -1,-7385, 2, 1, 0,-2, -1,-4372, 0, 0, 0, 0, 2,-2, -1,-3726, -4,-4212, 1, 1, 1, 0, 1, 2618, 0, 0, 3, 0,-2, 0, -1,-2241, 0, 0, 4, 0,-3, 0, 1, 1868, 0, 0, 2,-1, 2, 0, 1, 1770, 0, 0, 0, 2, 1, 0, -1,-1617, 1, 1655, 1, 1,-1, 0, 1, 777, 0, 0, 2, 0, 3, 0, 1, 595, 0, 0, 2, 0, 1, 2, 0,-9902, 0, 0, 2, 0,-4, 0, 0, 9483, 0, 0, 2,-2, 1, 0, 0, 7517, 0, 0, 0, 1,-3, 0, 0,-6694, 0, 0, 4, 1,-1, 0, 0,-6352, 0, 0, 1, 0, 2, 0, 0,-5840, 0, 0, 1, 0, 0,-2, 0,-5833, 0, 0, 6, 0,-2, 0, 0, 5716, 0, 0, 2, 0,-2,-2, 0,-5606, 0, 0, 1,-1, 0, 0, 0,-5569, 0, 0, 0, 1, 3, 0, 0,-5459, 0, 0, 2, 0,-2, 2, 0,-5357, 0, 0, 2, 0,-1,-2, 0, 0, 8, 7516, 3, 0, 0, 0, 0, 0, -1,-4189, }; #define NMB 56 static short MB[6*NMB] = { /* Latitude D l' l F 1" .0001" */ 0, 0, 0, 1,18461, 2387, 0, 0, 1, 1, 1010, 1671, 0, 0, 1,-1, 999, 6936, 2, 0, 0,-1, 623, 6524, 2, 0,-1, 1, 199, 4837, 2, 0,-1,-1, 166, 5741, 2, 0, 0, 1, 117, 2607, 0, 0, 2, 1, 61, 9120, 2, 0, 1,-1, 33, 3572, 0, 0, 2,-1, 31, 7597, 2,-1, 0,-1, 29, 5766, 2, 0,-2,-1, 15, 5663, 2, 0, 1, 1, 15, 1216, 2, 1, 0,-1, -12, -941, 2,-1,-1, 1, 8, 8681, 2,-1, 0, 1, 7, 9586, 2,-1,-1,-1, 7, 4346, 0, 1,-1,-1, -6,-7314, 4, 0,-1,-1, 6, 5796, 0, 1, 0, 1, -6,-4601, 0, 0, 0, 3, -6,-2965, 0, 1,-1, 1, -5,-6324, 1, 0, 0, 1, -5,-3684, 0, 1, 1, 1, -5,-3113, 0, 1, 1,-1, -5, -759, 0, 1, 0,-1, -4,-8396, 1, 0, 0,-1, -4,-8057, 0, 0, 3, 1, 3, 9841, 4, 0, 0,-1, 3, 6745, 4, 0,-1, 1, 2, 9985, 0, 0, 1,-3, 2, 7986, 4, 0,-2, 1, 2, 4139, 2, 0, 0,-3, 2, 1863, 2, 0, 2,-1, 2, 1462, 2,-1, 1,-1, 1, 7660, 2, 0,-2, 1, -1,-6244, 0, 0, 3,-1, 1, 5813, 2, 0, 2, 1, 1, 5198, 2, 0,-3,-1, 1, 5156, 2, 1,-1, 1, -1,-3178, 2, 1, 0, 1, -1,-2643, 4, 0, 0, 1, 1, 1919, 2,-1, 1, 1, 1, 1346, 2,-2, 0,-1, 1, 859, 0, 0, 1, 3, -1, -194, 2, 1, 1,-1, 0,-8227, 1, 1, 0,-1, 0, 8042, 1, 1, 0, 1, 0, 8026, 0, 1,-2,-1, 0,-7932, 2, 1,-1,-1, 0,-7910, 1, 0, 1, 1, 0,-6674, 2,-1,-2,-1, 0, 6502, 0, 1, 2, 1, 0,-6388, 4, 0,-2,-1, 0, 6337, 4,-1,-1,-1, 0, 5958, 1, 0, 1,-1, 0,-5889, }; #define NLRT 15 static short LRT[8*NLRT] = { /* Multiply by T Longitude Radius D l' l F .1" .00001" .1km .00001km */ 0, 1, 0, 0, 16, 7680, -1,-2302, 2,-1,-1, 0, -5,-1642, 3, 8245, 2,-1, 0, 0, -4,-1383, 5, 1395, 0, 1,-1, 0, 3, 7115, 3, 2654, 0, 1, 1, 0, 2, 7560, -2,-6396, 2, 1,-1, 0, 0, 7118, 0,-6068, 2, 1, 0, 0, 0, 6128, 0,-7754, 1, 1, 0, 0, 0,-4516, 0, 4194, 2,-2, 0, 0, 0,-4048, 0, 4970, 0, 2, 0, 0, 0, 3747, 0, 0, 2,-2,-1, 0, 0,-3707, 0, 0, 2,-1, 1, 0, 0,-3649, 0, 3222, 0, 1,-2, 0, 0, 2438, 0, 0, 2,-1,-2, 0, 0,-2165, 0, 0, 0, 1, 2, 0, 0, 1923, 0, 0, }; #define NBT 7 static short BT[5*NBT] = { /* Multiply by T Latitude D l' l F .00001" */ 2,-1, 0,-1, -7430, 2, 1, 0,-1, 3043, 2,-1,-1, 1, -2229, 2,-1, 0, 1, -1999, 2,-1,-1,-1, -1869, 0, 1,-1,-1, 1696, 0, 1, 0, 1, 1623, }; #define NLRT2 5 static short LRT2[6*NLRT2] = { /* Multiply by T^2 Longitude Radius D l' l F .00001" .00001km */ 0, 1, 0, 0, 487, 0, 2,-1,-1, 0, -150, 111, 2,-1, 0, 0, -120, 149, 0, 1,-1, 0, 108, 0, 0, 1, 1, 0, 80, 0, }; /* Rate of motion in R.A. and Dec., computed by this program. */ extern double dradt, ddecdt; /* Obliquity of the ecliptic, computed by epsiln(): */ extern double coseps, sineps; /* Answers posted by angles(): */ extern double pq, qe, ep; /* Nutation, computed by nutlo(): */ extern double nutl, nuto; /* The following times are set up by update() and refer * to the same instant. The distinction between them * is required by altaz(). */ extern double TDT; /* Ephemeris time */ extern double UT; /* Universal time */ extern double ss[5][8]; /* see nutate.c */ extern double cc[5][8]; #define DEBUG 0 static double l = 0.0; /* Moon's ecliptic longitude */ static double B = 0.0; /* Ecliptic latitude */ static double p = 0.0; /* Parallax */ static double ra = 0.0; /* Right Ascension */ static double dec = 0.0; /* Declination */ double moonpol[3]; double moonpp[3]; /* Beware of atan2()! */ double floor(), sin(), cos(), asin(), mods3600(); /* equatorial radius of the earth, in au */ #define Rearth 4.263523e-5 /* in kilometers */ #define Kearth 6378.14 /* Program begins: */ domoon() { int i; double ra0, dec0, r; double x, y, z, lon0; double pp[3], qq[3]; double acos(); /* Compute obliquity of the ecliptic, coseps, and sineps */ epsiln(TDT); /* Run the orbit calculation twice, at two different times, * in order to find the rate of change of R.A. and Dec. */ /* Calculate for 0.001 day ago */ i = prtflg; prtflg = 0; /* disable display */ moonll(TDT-0.001); lonlat( rearth, TDT, eapolar ); /* precess earth to date */ prtflg = i; /* reenable display */ ra0 = ra; dec0 = dec; lon0 = l; /* Calculate for present instant. */ moonll(TDT); /* The rates of change. These are used by altaz() to * correct the time of rising, transit, and setting. */ dradt = 1000.0*(ra-ra0); ddecdt = 1000.0*(dec-dec0); /* Post the ecliptic longitude and latitude, in radians, * and the radius in au. */ obpolar[0] = l; obpolar[1] = B; r = 1.0/sin(p); /* distance in earth-radii */ ra0 = 4.263523e-5*r; /* factor is radius of earth in au */ obpolar[2] = ra0; /* Rate of change in longitude, degrees per day * used for phase of the moon */ lon0 = 1000.0*RTD*(l - lon0); /* convert to ecliptic rectangular coordinates */ z = ra0 * cos(B); x = z * cos(l); y = z * sin(l); z = ra0 * sin(B); /* convert to equatorial coordinates */ pp[0] = x; pp[1] = y * coseps - z * sineps; pp[2] = y * sineps + z * coseps; /* Find sun-moon-earth angles */ precess( rearth, TDT, -1 ); for( i=0; i<3; i++ ) qq[i] = rearth[i] + pp[i]; angles( pp, qq, rearth ); /* Display answers */ if( prtflg ) { printf( "Apparent geocentric longitude %.3lf deg", RTD*(l+nutl) ); printf( " latitude %.3lf deg\n", RTD*B ); printf( "Distance %.5lf Earth-radii\n", r ); printf( "Horizontal parallax" ); dms( p ); printf( "Semidiameter" ); x = 0.272453 * p + 0.0799/RTS; /* AA page L6 */ dms( x ); x = RTD * acos(-ep); printf( "\nElongation from sun %.2lf deg,", x ); x = 0.5 * (1.0 + pq); printf( " Illuminated fraction %.2lf\n", x ); /* Find phase of the Moon by comparing Moon's longitude * with Earth's longitude. * * The number of days before or past indicated phase is * estimated by assuming the true longitudes change linearly * with time. These rates are estimated for the date, but * do not stay constant. The error can exceed 0.15 day in 4 days. */ /* Apparent longitude of sun. * Moon's longitude was already corrected for light time. */ y = eapolar[0] - 20.496/(RTS*eapolar[2]); x = obpolar[0] - y; x = modtp( x ) * RTD; /* difference in longitude */ i = x/90; /* number of quarters */ x = (x - i*90.0); /* phase angle mod 90 degrees */ /* days per degree of phase angle */ z = 1.0/(lon0 - (0.9856/eapolar[2])); if( x > 45.0 ) { y = -(x - 90.0)*z; if( y > 1.0 ) printf( "Phase %.1lf days before ", y ); else printf( "Phase %.2lf days before ", y ); i = (i+1) & 3; } else { y = x*z; if( y > 1.0 ) printf( "Phase %.1lf days past ", y ); else printf( "Phase %.2lf days past ", y ); } switch(i) { case 0: printf( "Full Moon\n" ); break; case 1: printf( "Third Quarter\n" ); break; case 2: printf( "New Moon\n" ); break; case 3: printf( "First Quarter\n" ); break; } } /* if prtflg */ printf( " Apparent: R.A." ); hms(ra); printf( "Declination" ); dms(dec); printf( "\n" ); /* Compute and display topocentric position (altaz.c) */ pp[0] = ra; pp[1] = dec; pp[2] = r * Rearth; altaz( pp, UT ); } /* Orbit calculation begins. */ static double f = 0.0; static double g = 0.0; static double LP = 0.0; static double M = 0.0; static double MP = 0.0; static double D = 0.0; static double NF = 0.0; static double OM = 0.0; static double Ve = 0.0; static double Ea = 0.0; static double Ma = 0.0; static double Ju = 0.0; static double T = 0.0; static double T2 = 0.0; /* Calculate apparent latitude, longitude, and horizontal parallax * of the Moon at Julian date J. */ moonll(J) double J; { /* Note, the time scale is supposed to be TDB but the difference * TDT - TDB is negligible compared to the truncation error * of the trigonometric series. */ T = (J-J2000)/36525.0; T2 = T*T; /* Program is broken up to get it thru micro compiler */ moon1(); moon2(); moon3(); moon4(1); /* Correct for nutation at date TDT. */ nutate( TDT, moonpp ); ra = zatan2(moonpp[0],moonpp[1]); dec = asin(moonpp[2]); } /* Calculate geometric coordinates of Moon * without light time or nutation correction. */ gmoon(J, rect, pol) double J; double rect[], pol[]; { double r; int i; epsiln(J); T = (J-J2000)/36525.0; T2 = T*T; moon1(); moon2(); moon3(); moon4(0); r = moonpol[2]; for( i=0; i<3; i++ ) { rect[i] = moonpp[i] * r; pol[i] = moonpol[i]; } } /* Reduce arc seconds modulo 360 degrees * answer in arc seconds */ double mods3600(x) double x; { double y; double floor(); y = x - 1296000. * floor( x/1296000.); return(y); } moon1() { double a; /* Mean longitude of moon */ LP = mods3600( 1732564372.83264 * T + 785939.95571 ); /* Mean distance of moon from its ascending node = F */ NF = mods3600( 1739527263.0983 * T + 335779.55755 ); /* Mean anomaly of sun = l' */ M = mods3600( 129596581.0474 * T + 1287104.79306 ); /* Mean anomaly of moon = l */ MP = mods3600( 1717915923.4728 * T + 485868.28096 ); /* Mean elongation of moon = D */ D = mods3600( 1602961601.4603 * T + 1072260.73512 ); /* Longitude of ascending node of moon = Omega */ /*OM = mods3600( -6962890.2656 * T + 450160.39816 );*/ /* Mean longitudes of planets */ Ve = mods3600( 210664136.43355 * T + 655127.28305 ); Ea = mods3600( 129597742.2758 * T + 361679.22059 ); Ma = mods3600( 68905077.59284 * T + 1279559.78866 ); Ju = mods3600( 10925660.42861 * T + 123665.34212 ); sscc( 0, STR*D, 6 ); sscc( 1, STR*M, 4 ); sscc( 2, STR*MP, 4 ); sscc( 3, STR*NF, 4 ); moonpol[0] = 0.0; moonpol[1] = 0.0; moonpol[2] = 0.0; /* terms in T^2 */ chewm( LRT2, NLRT2, 4, 2, moonpol ); moonpol[0] *= T; moonpol[2] *= T; /* terms in T */ chewm( BT, NBT, 4, 4, moonpol ); chewm( LRT, NLRT, 4, 1, moonpol ); f = 18 * Ve - 16 * Ea; g = STR*(f - MP + 412435.8); moonpol[0] += 25425. * sin(g); a = T / 10.0; /* set amplitude scale of 1.0 = 10^-4 arcsec */ moonpol[0] *= a; moonpol[1] *= a; moonpol[2] *= a; } moon2() { /* terms in T^0 */ g = STR*(f - MP + 95553.396); l = 14.24883 * sin(g); g = STR*(LP-NF+3.384); l += 7.06304 * sin(g); g = STR*(2*(Ea-Ju+D)-MP+648431.172); l += 1.14307 * sin(g); g = STR*( 4*Ea - 8*Ma + 3*Ju + 1029553.452); l += 0.90114 * sin(g); g = STR*(Ve-Ea+648035.568); l += 0.82155 * sin(g); g = STR*(f - 2*MP + 95555.664); l += 0.78811 * sin(g); g = STR*(f + 95564.16); l += 0.73930 * sin(g); g = STR*(3*(Ve-Ea)+2*D-MP+647933.184); l += 0.64371 * sin(g); g = STR*(Ea-Ju+4424.04); l += 0.63880 * sin(g); g = STR*(10*Ve-3*Ea-MP+1199899.836); l += 0.56341 * sin(g); g = STR*(LP+648002.556); B = 8.04508 * sin(g); g = STR*(Ea+D+996048.252); B += 1.51021 * sin(g); g = STR*(f - MP + NF + 95554.332); B += 0.63037 * sin(g); g = STR*(f - MP - NF + 95553.792); B += 0.63014 * sin(g); moonpol[0] += 10000.0 * l; moonpol[1] += 10000.0 * B; } moon3() { LP += ((-0.00005522 * T + 0.006681)*T - 4.7763)*T2; MP += ((-0.00024470 * T + 0.051651)*T + 32.3893)*T2; M += ( 0.000147 * T - 0.5529)*T2; NF += (( 0.00000417 * T - 0.001021)*T - 12.2505)*T2; D += ((-0.00003184 * T + 0.006595)*T - 5.8681)*T2; sscc( 0, STR*D, 6 ); sscc( 1, STR*M, 4 ); sscc( 2, STR*MP, 4 ); sscc( 3, STR*NF, 4 ); /* terms in T^0 */ chewm( LR, NLR, 4, 1, moonpol ); chewm( MB, NMB, 4, 3, moonpol ); moonpol[0] /= 10000.0; /* change scale to 1.0 = 1 arcsec */ moonpol[0] += LP; moonpol[1] /= 10000.0; moonpol[2] = moonpol[2] / 10000.0 + 385000.52899; /* kilometers */ } /* Comput final ecliptic polar coordinates * and convert to equatorial rectangular coordinates */ moon4(ltflag) int ltflag; { double cosB, sinB, cosL, sinL, sp; sp = Kearth/moonpol[2]; p = asin( sp ); moonpol[2] /= 1.4959787066e8; /* Kilometers per au */ l = STR * mods3600( moonpol[0] ); /* Light time correction to longitude, * about 0.7". */ if( ltflag ) l -= DTR * 0.0118 * sp; moonpol[0] = l; B = STR * moonpol[1]; moonpol[1] = B; /* printf( "%.7e %.7e %.7e\n", l, B, p ); */ cosB = cos(B); sinB = sin(B); cosL = cos(l); sinL = sin(l); moonpp[0] = cosB*cosL; moonpp[1] = coseps*cosB*sinL - sineps*sinB; moonpp[2] = sineps*cosB*sinL + coseps*sinB; } /* Program to step through the perturbation table */ chewm( p, nlines, nangles, typflg, ans ) register short *p; int nlines, nangles; int typflg; double ans[]; { int i, j, k, k1, m; double cu, su, cv, sv, f; for( i=0; i 0 */ /* sin, cos (k*angle) from lookup table */ su = ss[m][k-1]; cu = cc[m][k-1]; if( j < 0 ) su = -su; /* negative angle factor */ if( k1 == 0 ) { /* Set sin, cos of first angle. */ sv = su; cv = cu; k1 = 1; } else { /* Combine angles by trigonometry. */ f = su*cv + cu*sv; cv = cu*cv - su*sv; sv = f; } } } /* Accumulate */ switch( typflg ) { /* large longitude and radius */ case 1: j = *p++; k = *p++; ans[0] += (10000.0 * j + k) * sv; j = *p++; k = *p++; if( k ) { ans[2] += (10000.0 * j + k) * cv; } break; /* longitude and radius */ case 2: ans[0] += *p++ * sv; ans[2] += *p++ * cv; break; /* large latitude */ case 3: j = *p++; k = *p++; ans[1] += ( 10000.0*j + k)*sv; break; /* latitude */ case 4: ans[1] += *p++ * sv; break; } } }