/* * 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 /* Trigonometric functions in degrees. */ #ifndef PI # define PI 3.1415926535897932385 /* Should be in math.h . */ #endif /* PI */ #define RAD (PI/180.0) #define sind( x ) ( sin( RAD * (x) ) ) #define cosd( x ) ( cos( RAD * (x) ) ) #define MAGMOON -9.9 extern double range(); /* 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. */ moon_pos( t, epli ) double t, epli; { double t_2 = t*t, t_3 = t*t*t, e, e_2, l_prime, /* Moon's longitude. */ 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. */ omega1, omega2, venus_term, /* Great Venus term. */ lambda, /* Ecliptical longitude of the Moon's center. */ b, beta, /* Ecliptical latitude of the Moon's center. */ parallax, /* Equatorial horizontal parallax of the Moon. */ tmp; double phase; /* phase of the moon */ static char namestring[] = "Moon, 180."; /* Compute the mean values for the terms. */ l_prime = 270.434164 + 481267.8831 * t - 0.001133 * t_2 + 0.0000019 * t_3; m = 358.475833 + 35999.0498 * t - 0.000150 * t_2 - 0.0000033 * t_3; m_prime = 296.104608 + 477198.8491 * t + 0.009192 * t_2 + 0.0000144 * t_3; d = 350.737486 + 445267.1142 * t - 0.001436 * t_2 + 0.0000019 * t_3; f = 11.250889 + 483202.0251 * t - 0.003211 * t_2 - 0.0000003 * t_3; omega = 259.183275 - 1934.1420 * t + 0.002078 * t_2 + 0.0000022 * t_3; /* Additive terms. */ tmp = sind( 51.2 + 20.2 * t ); l_prime += 0.000233 * tmp; m += -0.001778 * tmp; m_prime += 0.000817 * tmp; d += 0.002011 * tmp; venus_term = 0.003964 * sind( 346.560 + 132.870 * t - 0.0091731 * t_2 ); l_prime += venus_term; m_prime += venus_term; d += venus_term; f += venus_term; tmp = sind( omega ); l_prime += 0.001964 * tmp; m_prime += 0.002541 * tmp; d += 0.001964 * tmp; f += -0.024691 * tmp; f += -0.004328 * sind( omega + 275.05 - 2.30 * t ); e = 1 - 0.002495 * t - 0.00000752 * t_2; 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 ); /* broke into pieces so a micro compute compiler could handel it -- Bob L */ /* Calculate lambda, the ecliptical longitude of the Moon's center. */ lambda = l_prime + 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 - m_prime ); lambda += 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 + m_prime ) + 0.015326 * sind( 2*d - 2*f ) - 0.012528 * sind( 2*f + m_prime ) - 0.010980 * sind( 2*f - m_prime ); lambda += 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( m - m_prime + 2*d ) - e * 0.006783 * sind( 2*d + m ) + 0.005162 * sind( m_prime - d ) + e * 0.005000 * sind( m + d ) + e * 0.004049 * sind( m_prime - m + 2*d ) + 0.003996 * sind( 2*m_prime + 2*d ) + 0.003862 * sind( 4*d ); lambda += 0.003665 * sind( 2*d - 3*m_prime ) + e * 0.002695 * sind( 2*m_prime - m ) + 0.002602 * sind( m_prime - 2*f - 2*d ) + e * 0.002396 * sind( 2*d - m - 2*m_prime ) - 0.002349 * sind( m_prime + d ) + 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( m_prime + 2*d - 2*f ) - 0.001595 * sind( 2*f + 2*d ); lambda += e * 0.001220 * sind( 4*d - m - m_prime ) - 0.001110 * sind( 2*m_prime + 2*f ) + 0.000892 * sind( m_prime - 3*d ) - e * 0.000811 * sind( m + m_prime + 2*d ) + e * 0.000761 * sind( 4*d - m - 2*m_prime ) + e_2 * 0.000717 * sind( m_prime - 2*m ) + e_2 * 0.000704 * sind( m_prime - 2*m -2*d ) + e * 0.000693 * sind( m - 2*m_prime + 2*d ) + e * 0.000598 * sind( 2*d - m - 2*f ) + 0.000550 * sind( m_prime + 4*d ) + 0.000538 * sind( 4*m_prime ) + e * 0.000521 * sind( 4*d - m ) + 0.000486 * sind( 2*m_prime - d ); lambda = range( lambda ); /* Calculate beta, the ecliptical latitude of the Moon's center. */ b = 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 + f - m_prime ) + 0.046272 * sind( 2*d - f - m_prime ) + 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 - f - 2*m_prime ); b += 0.004200 * sind( 2*d + f + m_prime ) + e * 0.003372 * sind( f - m - 2*d ) + e * 0.002472 * sind( 2*d + f - m - m_prime ) + e * 0.002222 * sind( 2*d + f - m ) + e * 0.002072 * sind( 2*d - f - m - m_prime ) + e * 0.001877 * sind( f - m + m_prime ) + 0.001828 * sind( 4*d - f - m_prime ) - e * 0.001803 * sind( f + m ) - 0.001750 * sind( 3*f ) + e * 0.001570 * sind( m_prime - m - f ) - 0.001487 * sind( f + d ) - e * 0.001481 * sind( f + m + m_prime ); b += e* 0.001417 * sind( f - m - m_prime ) + e * 0.001350 * sind( f - m ) + 0.001330 * sind( f - d ) + 0.001106 * sind( f + 3*m_prime ) + 0.001020 * sind( 4*d - f ) + 0.000833 * sind( f + 4*d - m_prime ) + 0.000781 * sind( m_prime - 3*f ) + 0.000670 * sind( f + 4*d - 2*m_prime ) + 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 ); b += 0.000450 * sind( 2*m_prime - f - 2*d ) + 0.000439 * sind( 3*m_prime - f ) + 0.000423 * sind( f + 2*d + 2*m_prime ) + 0.000422 * sind( 2*d - f - 3*m_prime ) - e * 0.000367 * sind( m + f + 2*d - m_prime ) - e * 0.000353 * sind( m + f + 2*d ) + 0.000331 * sind( f + 4*d ) + e * 0.000317 * sind( 2*d + f - m + m_prime ) + e_2 * 0.000306 * sind( 2*d - 2*m - f ) - 0.000283 * sind( m_prime + 3*f ); omega1 = 0.0004664 * cosd( omega ); omega2 = 0.0000754 * cosd( omega + 275.05 - 2.30 * t ); beta = b * ( 1 - omega1 - omega2 ); /* roughly calculate the phase of the moon. Added by ccount Sun longi is about 36000 * T0 + 279 Phase is 0 for full moon */ tmp = 36000.0 * t +279 - lambda; phase = 180.0 - (tmp - ((double)((int)(tmp/360.0))*360.0)); phase = (phase < 0.0) ? 360 + phase:phase; sprintf(namestring, "Moon, %3.0f", phase); /* Transform to right ascension and declination, and output the * result. */ geo_trans( lambda, beta, epli, 0.0, MAGMOON, "PL", namestring ); }