/* code to compute lunar sunrise position and local sun angle. */ #include #include #if defined(__STDC__) #include #endif #include "P_.h" #include "astro.h" static void Librations P_((double RAD, double LAMH, double BH, double OM, double F, double L, double L1, double *L0, double *B0)); static void Moon P_((double RAD, double T, double T2, double LAM0, double R, double M, double *F, double *L1, double *OM, double *LAM, double *B, double *DR, double *LAMH, double *BH)); static void Sun P_((double RAD, double T, double T2, double *L, double *M, double *R, double *LAM0)); /* given a Julian date and a lunar location, find selenographic colongitude of * rising sun, lunar latitude of subsolar point, illuminated fraction, and alt * of sun at the given location. Any pointer may be 0 if not interested. * From Bruning and Talcott, October 1995 _Astronomy_, page 76. * N.B. lunar coordinates use +E, but selenograhic colongs are +W. */ void moon_colong (jd, lt, lg, cp, kp, ap, sp) double jd; /* jd */ double lt, lg; /* lat/long of location on moon, rads +N +E */ double *cp; /* selenographic colongitude (-lng of rising sun), rads */ double *kp; /* illuminated fraction of surface from Earth */ double *ap; /* sun altitude at location, rads */ double *sp; /* lunar latitude of subsolar point, rads */ { double RAD = .0174533; double T; double T2; double L, M, R, LAM0; double F, L1, OM, LAM, B, DR, LAMH, BH; double L0, B0; double TEMP; double C0; double PSI; double NUM, DEN; double I, K; double THETA, ETA; double H; T = (jd - 2451545)/36525.0; T2 = T * T; Sun(RAD, T, T2, &L, &M, &R, &LAM0); Moon(RAD, T, T2, LAM0, R, M, &F, &L1, &OM, &LAM, &B, &DR, &LAMH, &BH); Librations(RAD, LAMH, BH, OM, F, L, L1, &L0, &B0); if (sp) *sp = B0; TEMP = L0 / 360; L0 = ((TEMP) - (int)(TEMP)) * 360; if (L0 < 0) L0 = L0 + 360; if (L0 <= 90) C0 = 90 - L0; else C0 = 450 - L0; if (cp) { *cp = degrad(C0); range (cp, 2*PI); /* prefer 0..360 +W */ } if (kp) { TEMP = cos(B * RAD) * cos(LAM - LAM0 * RAD); PSI = acos(TEMP); NUM = R * sin(PSI); DEN = DR - R * TEMP; I = atan(NUM / DEN); if (NUM * DEN < 0) I = I + 3.14159; if (NUM < 0) I = I + 3.14159; K = (1 + cos(I)) / 2; *kp = K; } if (ap) { THETA = lt; ETA = lg; C0 = C0 * RAD; TEMP = sin(B0) * sin(THETA) + cos(B0) * cos(THETA) * sin(C0+ETA); H = asin(TEMP); *ap = H; } } static void Librations (RAD, LAMH, BH, OM, F, L, L1, L0, B0) double RAD; double LAMH; double BH; double OM; double F; double L; double L1; double *L0; double *B0; { double I, PSI, W, NUM, DEN, A, TEMP; /* inclination of lunar equator */ I = 1.54242 * RAD; /* nutation in longitude, in arcseconds */ PSI = -17.2 * sin(OM) - 1.32 * sin(2 * L) - .23 * sin(2 * L1) + .21 * sin(2 * OM); PSI = PSI * RAD / 3600; /* optical librations */ W = (LAMH - PSI) - OM; NUM = sin(W) * cos(BH) * cos(I) - sin(BH) * sin(I); DEN = cos(W) * cos(BH); A = atan(NUM / DEN); if (NUM * DEN < 0) A = A + 3.14159; if (NUM < 0) A = A + 3.14159; *L0 = (A - F) / RAD; TEMP = -sin(W) * cos(BH) * sin(I) - sin(BH) * cos(I); *B0 = asin(TEMP); } static void Moon (RAD, T, T2, LAM0, R, M, F, L1, OM, LAM, B, DR, LAMH, BH) double RAD; double T; double T2; double LAM0; double R; double M; double *F; double *L1; double *OM; double *LAM; double *B; double *DR; double *LAMH; double *BH; { double T3, M1, D2, SUMR, SUML, DIST; T3 = T * T2; /* argument of the latitude of the Moon */ *F = (93.2721 + 483202 * T - .003403 * T2 - T3 / 3526000) * RAD; /* mean longitude of the Moon */ *L1 = (218.316 + 481268. * T) * RAD; /* longitude of the ascending node of Moon's mean orbit */ *OM = (125.045 - 1934.14 * T + .002071 * T2 + T3 / 450000) * RAD; /* Moon's mean anomaly */ M1 = (134.963 + 477199 * T + .008997 * T2 + T3 / 69700) * RAD; /* mean elongation of the Moon */ D2 = (297.85 + 445267 * T - .00163 * T2 + T3 / 545900) * 2 * RAD; /* Lunar distance */ SUMR = -20954 * cos(M1) - 3699 * cos(D2 - M1) - 2956 * cos(D2); *DR = 385000 + SUMR; /* geocentric latitude */ *B = 5.128 * sin(*F) + .2806 * sin(M1 + *F) + .2777 * sin(M1 - *F) + .1732 * sin(D2 - *F); SUML = 6.289 * sin(M1) + 1.274 * sin(D2 - M1) + .6583 * sin(D2) + .2136 * sin(2 * M1) - .1851 * sin(M) - .1143 * sin(2 * *F); *LAM = *L1 + SUML * RAD; DIST = *DR / R; *LAMH = (LAM0 + 180 + DIST * cos(*B) * sin(LAM0 * RAD - *LAM) / RAD) * RAD; *BH = DIST * *B * RAD; } static void Sun (RAD, T, T2, L, M, R, LAM0) double RAD; double T; double T2; double *L; double *M; double *R; double *LAM0; { double T3, C, V, E, THETA, OM; T3 = T2 * T; /* mean longitude of the Sun */ *L = 280.466 + 36000.8 * T; /* mean anomaly of the Sun */ *M = 357.529 + 35999 * T - .0001536 * T2 + T3 / 24490000; *M = *M * RAD; /* correction for Sun's elliptical orbit */ C = (1.915 - .004817 * T - .000014 * T2) * sin(*M) + (.01999 - .000101 * T) * sin(2 * *M) + .00029 * sin(3 * *M); /* true anomaly of the Sun */ V = *M + C * RAD; /* eccentricity of Earth's orbit */ E = .01671 - .00004204 * T - .0000001236 * T2; /* Sun-Earth distance */ *R = .99972 / (1 + E * cos(V)) * 145980000; /* true geometric longitude of the Sun */ THETA = *L + C; /* apparent longitude of the Sun */ OM = 125.04 - 1934.1 * T; *LAM0 = THETA - .00569 - .00478 * sin(OM * RAD); } #ifdef TESTCOLONG /* insure 0 <= *v < r. */ void range (v, r) double *v, r; { *v -= r*floor(*v/r); } /* To be sure the program is functioning properly, try the test case * 2449992.5 (1 Oct 1995): the colongitude should be 3.69 degrees. */ int main (int ac, char *av[]) { double jd, lt, lg; double c, k, a; if (ac != 2) { fprintf (stderr, "%s: JD\n", av[0]); exit (1); } jd = atof(av[1]); printf ("Latitude of lunar feature: "); fscanf (stdin, "%lf", <); lt = degrad(lt); printf ("Longitude: "); fscanf (stdin, "%lf", &lg); lg = degrad(lg); moon_colong (jd, lt, lg, &c, &k, &a); printf ("Selenographic colongitude is %g\n", raddeg(c)); printf ("The illuminated fraction of the Moon is %g\n", k); printf ("Altitude of Sun above feature is %g\n", raddeg(a)); return (0); } #endif /* For RCS Only -- Do Not Edit */ static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: mooncolong.c,v $ $Date: 1996/08/26 19:42:16 $ $Revision: 1.1 $ $Name: $"};