/******************************************************************************* ** MOON-POS.C - Determine the position of the moon. ** ** ** ** Author: Jim Cobb ** ** jcobb@utah.edu.cs ** ** Computer Science Dept. ** ** University of Utah ** ** Date: Fri Mar 25 1988 ** ** ** ** Copyright (c) 1988, Jim Cobb, All rights reserved. ** ** ** ** (This copyright prevents you from selling the program - the author grants ** ** anyone the right to copy and install the program on any machine it will ** ** run on) ** ** ** ** Reference: Jean Meeus, "Astronomical Formulae for Calculators," chapter 30 ** ** Meeus states that this procedure has an accuracy of 10" in the longitude ** ** of the moon, 3" in its latitude and 0.2" in its parallax. ** *******************************************************************************/ #include #include #ifndef PI # define PI 3.14159265358979323846 /* Should be in */ #endif /* PI */ #ifndef RAD # define RAD 0.01745329251994329651 /* (PI / 180 deg) */ #endif /* RAD */ /******************************************************************************* ** Trigonometric functions in degrees. ** *******************************************************************************/ #define sind(x) sin(RAD*(x)) #define cosd(x) cos(RAD*(x)) #define MAGMOON -9.9 /******************************************************************************* ** "moon_pos" takes a time argument, t, the number of Julian centuries from ** ** 1900 January 0.5 ET, and the obliquity of the ecliptic, epli. ** ** It adds information about the moon's position to the logfile. ** *******************************************************************************/ void moon_pos(t,epli) double t,epli; { double e,e_2, m, /* Sun's anomaly. */ m_prime, /* Moon's anomaly. */ d, /* Moon's elongation. */ f, /* Distance of Moon from its ascending node. */ omega, /* Longitude of Moon's ascending node. */ venus_term, /* Great Venus term. */ lambda, /* Ecliptical longitude of the Moon's center. */ beta; /* Ecliptical latitude of the Moon's center. */ double phase; /* Phase of the moon */ #ifdef EUTSCH static char namestring[] = "Mond, 180."; #else static char namestring[] = "Moon, 180."; #endif double range(),poly(); int to_int(); void geo_trans(); /******************************************************************************* ** Compute the mean values for the terms. ** *******************************************************************************/ lambda = poly(270.434164,481267.8831,-0.001133,-0.0000019,t); m = poly(358.475833, 35999.0498,-0.000150, 0.0000033,t); m_prime = poly(296.104608,477198.8491, 0.009192, 0.0000144,t); d = poly(350.737486,445267.1142,-0.001436,-0.0000019,t); f = poly( 11.250889,483202.0251,-0.003211, 0.0000003,t); omega = poly(259.183275, -1934.1420,-0.002078,-0.0000022,t); /******************************************************************************* ** Additive terms. ** *******************************************************************************/ e_2 = sind( 51.2 + 20.2 * t ); lambda += 0.000233 * e_2; m += -0.001778 * e_2; m_prime += 0.000817 * e_2; d += 0.002011 * e_2; venus_term = 0.003964 * sind(poly(346.560,132.870,-0.0091731,0.,t)); lambda += venus_term; m_prime += venus_term; d += venus_term; f += venus_term; e_2 = sind( omega ); lambda += 0.001964 * e_2; m_prime += 0.002541 * e_2; d += 0.001964 * e_2; f += -0.024691 * e_2; f += -0.004328 * sind( omega + 275.05 - 2.30 * t ); e = poly(1.,-0.002495,-0.00000752,0.,t); e_2 = e*e; /******************************************************************************* ** Bring these angles within 0 to 360 degrees. ** *******************************************************************************/ m = range( m ); m_prime = range( m_prime ); d = range( d ); f = range( f ); omega = range( omega ); /******************************************************************************* ** Calculate lambda, the ecliptical longitude of the Moon's center. ** *******************************************************************************/ lambda += 6.288750 * sind( m_prime ) + 1.274018 * sind( 2.*d - m_prime ) + 0.658309 * sind( 2.*d ) + 0.213616 * sind( 2.*m_prime ) - e * 0.185596 * sind( m ) - 0.114336 * sind( 2.*f ) + 0.058793 * sind( 2.*d - 2.*m_prime ) + e * 0.057212 * sind( 2.*d - m_prime - m ) + 0.053320 * sind( 2.*d + m_prime ) + e * 0.045874 * sind( 2.*d - m ) + e * 0.041024 * sind( m_prime - m ) - 0.034718 * sind( d ) - e * 0.030465 * sind( m_prime + m ) + 0.015326 * sind( 2.*d - 2.*f ) - 0.012528 * sind( m_prime + 2.*f ) - 0.010980 * sind( - m_prime + 2.*f ) + 0.010674 * sind( 4.*d - m_prime ) + 0.010034 * sind( 3.*m_prime ) + 0.008548 * sind( 4.*d - 2.*m_prime ) - e * 0.007910 * sind( 2.*d - m_prime + m ) - e * 0.006783 * sind( 2.*d + m ) + 0.005162 * sind( - d + m_prime ) + e * 0.005000 * sind( d + m ) + e * 0.004049 * sind( 2.*d + m_prime - m ) + 0.003996 * sind( 2.*d + 2.*m_prime ) + 0.003862 * sind( 4.*d ) + 0.003665 * sind( 2.*d - 3.*m_prime ) + e * 0.002695 * sind( 2.*m_prime - m ) + 0.002602 * sind( -2.*d + m_prime - 2.*f ) + e * 0.002396 * sind( 2.*d - 2.*m_prime - m ) - 0.002349 * sind( d + m_prime ) + e_2 * 0.002249 * sind( 2.*d - 2.*m ) - e * 0.002125 * sind( 2.*m_prime + m ) - e_2 * 0.002079 * sind( 2.*m ) + e_2 * 0.002059 * sind( 2.*d - m_prime - 2.*m ) - 0.001773 * sind( 2.*d + m_prime - 2.*f ) - 0.001595 * sind( 2.*d + 2.*f ) + e * 0.001220 * sind( 4.*d - m_prime - m ) - 0.001110 * sind( 2.*m_prime + 2.*f ) + 0.000892 * sind( -3.*d + m_prime ) - e * 0.000811 * sind( 2.*d + m_prime + m ) + e * 0.000761 * sind( 4.*d - 2.*m_prime - m ) + e_2 * 0.000717 * sind( m_prime - 2.*m ) + e_2 * 0.000704 * sind( -2.*d + m_prime - 2.*m ) + e * 0.000693 * sind( 2.*d - 2.*m_prime + m ) + e * 0.000598 * sind( 2.*d - m - 2.*f ) + 0.000550 * sind( 4.*d + m_prime ) + 0.000538 * sind( 4.*m_prime ) + e * 0.000521 * sind( 4.*d - m ) + 0.000486 * sind( - d + 2.*m_prime ); lambda = range( lambda ); /******************************************************************************* ** Calculate beta, the ecliptical latitude of the Moon's center. ** *******************************************************************************/ beta = 5.128189 * sind( f ) + 0.280606 * sind( m_prime + f ) + 0.277693 * sind( m_prime - f ) + 0.173238 * sind( 2.*d - f ) + 0.055413 * sind( 2.*d - m_prime + f ) + 0.046272 * sind( 2.*d - m_prime - f ) + 0.032573 * sind( 2.*d + f ) + 0.017198 * sind( 2.*m_prime + f ) + 0.009267 * sind( 2.*d + m_prime - f ) + 0.008823 * sind( 2.*m_prime - f ) + e * 0.008247 * sind( 2.*d - m - f ) + 0.004323 * sind( 2.*d - 2.*m_prime - f ) + 0.004200 * sind( 2.*d + m_prime + f ) + e * 0.003372 * sind( -2.*d - m + f ) + e * 0.002472 * sind( 2.*d - m_prime - m + f ) + e * 0.002222 * sind( 2.*d - m + f ) + e * 0.002072 * sind( 2.*d - m_prime - m - f ) + e * 0.001877 * sind( m_prime - m + f ) + 0.001828 * sind( 4.*d - m_prime - f ) - e * 0.001803 * sind( m + f ) - 0.001750 * sind( 3.*f ) + e * 0.001570 * sind( m_prime - m - f ) - 0.001487 * sind( d + f ) - e * 0.001481 * sind( m_prime + m + f ) + e * 0.001417 * sind( - m_prime - m + f ) + e * 0.001350 * sind( - m + f ) + 0.001330 * sind( - d + f ) + 0.001106 * sind( 3.*m_prime + f ) + 0.001020 * sind( 4.*d - f ) + 0.000833 * sind( 4.*d - m_prime + f ) + 0.000781 * sind( m_prime - 3.*f ) + 0.000670 * sind( 4.*d - 2.*m_prime + f ) + 0.000606 * sind( 2.*d - 3.*f ) + 0.000597 * sind( 2.*d + 2.*m_prime - f ) + e * 0.000492 * sind( 2.*d + m_prime - m - f ) + 0.000450 * sind( -2.*d + 2.*m_prime - f ) + 0.000439 * sind( 3.*m_prime - f ) + 0.000423 * sind( 2.*d + 2.*m_prime + f ) + 0.000422 * sind( 2.*d - 3.*m_prime - f ) - e * 0.000367 * sind( 2.*d - m_prime + m + f ) - e * 0.000353 * sind( 2.*d + m + f ) + 0.000331 * sind( 4.*d + f ) + e * 0.000317 * sind( 2.*d + m_prime - m + f ) + e_2 * 0.000306 * sind( 2.*d - 2.*m - f ) - 0.000283 * sind( m_prime + 3.*f ); beta *= ( 1. - 0.0004664 * cosd( omega ) - 0.0000754 * cosd( omega + 275.05 - 2.30 * t )); /******************************************************************************* ** Roughly calculate the phase of the moon. ** ** Added by ccount ** ** Sun longitude is about 36000 * T0 + 279 ** ** Phase is 0 for full moon ** *******************************************************************************/ e_2 = 36000. * t + 279 - lambda; phase = 180. - (e_2 - ((double)(to_int(e_2/360.))*360.)); if(phase < 0.) phase += 360.; #ifdef EUTSCH sprintf(namestring,"Mond, %3.0f",phase); #else sprintf(namestring,"Moon, %3.0f",phase); #endif /******************************************************************************* ** Transform to right ascension and declination, and output the result. ** *******************************************************************************/ geo_trans(lambda,beta,epli,0.,MAGMOON,"PL",namestring ); }