home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 383.lha / Planets / moon_pos.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-25  |  8.2 KB  |  237 lines
  1. /* 
  2.  * moon-pos.c - Determine the position of the moon.
  3.  * 
  4.  * Author:    Jim Cobb
  5.  *         jcobb@utah.edu.cs
  6.  *         Computer Science Dept.
  7.  *         University of Utah
  8.  * Date:    Fri Mar 25 1988
  9.  * 
  10.  * Copyright (c) 1988, Jim Cobb
  11.  * All rights reserved
  12.  *
  13.  * (This copyright prevents you from selling the program - the
  14.  *  author grants anyone the right to copy and install the program
  15.  *  on any machine it will run on)
  16.  */
  17.  
  18. /* Reference: Jean Meeus, " Astronomical Formulae for Calculators,"
  19.  * chapter 30.  Meeus states that this procedure has an accuracy of 10"
  20.  * in the longitude of the moon, 3" in its latitude and 0.2" in its
  21.  * parallax.
  22.  */
  23.  
  24. #include <stdio.h>
  25. #include <math.h>
  26.  
  27. /* Trigonometric functions in degrees. */
  28. #ifndef PI
  29. #   define PI 3.1415926535897932385    /* Should be in math.h . */
  30. #endif /* PI */
  31. #define RAD    (PI/180.0)
  32. #define sind( x )    ( sin( RAD * (x) ) )
  33. #define cosd( x )    ( cos( RAD * (x) ) )
  34.  
  35. #define MAGMOON -9.9
  36.  
  37. extern double range();
  38.  
  39.  
  40. /* Moon_pos takes a time argument, t, the number of Julian centuries
  41.  * from 1900 January 0.5 ET, and the obliquity of the ecliptic, epli.
  42.  * It adds information about the moon's position to the logfile.
  43.  */
  44. moon_pos( t, epli )
  45. double t, epli;
  46. {
  47.     double t_2 = t*t, 
  48.        t_3 = t*t*t,
  49.        e, e_2,
  50.        l_prime,        /* Moon's longitude. */
  51.        m,            /* Sun's anomaly. */
  52.        m_prime,        /* Moon's anomaly. */
  53.        d,            /* Moon's elongation. */
  54.        f,            /* Distance of Moon from its ascending node. */
  55.        omega,        /* Longitude of Moon's ascending node. */
  56.        omega1, omega2,
  57.        venus_term,        /* Great Venus term. */
  58.        lambda,           /* Ecliptical longitude of the Moon's center. */
  59.        b, beta,            /* Ecliptical latitude of the Moon's center. */
  60.        parallax,          /* Equatorial horizontal parallax of the Moon. */
  61.        tmp;
  62.  
  63.     double phase;             /* phase of the moon */
  64.     static char namestring[] = "Moon, 180.";
  65.  
  66.     /* Compute the mean values for the terms. */
  67.     l_prime = 270.434164 + 481267.8831 * t - 0.001133 * t_2 + 0.0000019 * t_3;
  68.     m = 358.475833 + 35999.0498 * t - 0.000150 * t_2 - 0.0000033 * t_3;
  69.     m_prime = 296.104608 + 477198.8491 * t + 0.009192 * t_2 + 0.0000144 * t_3;
  70.     d = 350.737486 + 445267.1142 * t - 0.001436 * t_2 + 0.0000019 * t_3;
  71.     f = 11.250889 + 483202.0251 * t - 0.003211 * t_2 - 0.0000003 * t_3;
  72.     omega = 259.183275 - 1934.1420 * t + 0.002078 * t_2 + 0.0000022 * t_3;
  73.  
  74.     /* Additive terms. */
  75.     tmp = sind( 51.2 + 20.2 * t );
  76.     l_prime += 0.000233 * tmp;
  77.     m += -0.001778 * tmp;
  78.     m_prime += 0.000817 * tmp;
  79.     d += 0.002011 * tmp;
  80.  
  81.     venus_term = 0.003964 * 
  82.          sind( 346.560 + 132.870 * t - 0.0091731 * t_2 );
  83.     l_prime += venus_term;
  84.     m_prime += venus_term;
  85.     d += venus_term;
  86.     f += venus_term;
  87.  
  88.     tmp = sind( omega );
  89.     l_prime += 0.001964 * tmp;
  90.     m_prime += 0.002541 * tmp;
  91.     d += 0.001964 * tmp;
  92.     f += -0.024691 * tmp;
  93.  
  94.     f += -0.004328 * sind( omega + 275.05 - 2.30 * t );
  95.  
  96.     e = 1 - 0.002495 * t - 0.00000752 * t_2;
  97.     e_2 = e*e;
  98.  
  99.     /* Bring these angles within 0 to 360 degrees. */
  100.     m = range( m );
  101.     m_prime = range( m_prime );
  102.     d = range( d );
  103.     f = range( f );
  104.     omega = range( omega );
  105.  
  106.  
  107. /* broke into pieces so a micro compute compiler could handel it -- Bob L */
  108.     /* Calculate lambda, the ecliptical longitude of the Moon's center. */
  109.     lambda = l_prime +       6.288750 * sind( m_prime )
  110.              +       1.274018 * sind( 2*d - m_prime )
  111.              +       0.658309 * sind( 2*d )
  112.              +       0.213616 * sind( 2 * m_prime )
  113.              - e   * 0.185596 * sind( m )
  114.              -       0.114336 * sind( 2*f )
  115.              +       0.058793 * sind( 2*d - 2*m_prime )
  116.              + e   * 0.057212 * sind( 2*d - m - m_prime );
  117.  
  118.     lambda +=        0.053320 * sind( 2*d + m_prime )
  119.              + e   * 0.045874 * sind( 2*d - m )
  120.              + e   * 0.041024 * sind( m_prime - m )
  121.              -       0.034718 * sind( d )
  122.              - e   * 0.030465 * sind( m + m_prime )
  123.              +       0.015326 * sind( 2*d - 2*f )
  124.              -       0.012528 * sind( 2*f + m_prime )
  125.              -       0.010980 * sind( 2*f - m_prime );
  126.  
  127.     lambda +=        0.010674 * sind( 4*d - m_prime )
  128.              +       0.010034 * sind( 3*m_prime )
  129.              +       0.008548 * sind( 4*d - 2*m_prime )
  130.              - e   * 0.007910 * sind( m - m_prime + 2*d )
  131.              - e   * 0.006783 * sind( 2*d + m )
  132.              +       0.005162 * sind( m_prime - d )
  133.              + e   * 0.005000 * sind( m + d )
  134.              + e   * 0.004049 * sind( m_prime - m + 2*d )
  135.              +       0.003996 * sind( 2*m_prime + 2*d )
  136.              +       0.003862 * sind( 4*d );
  137.  
  138.     lambda +=        0.003665 * sind( 2*d - 3*m_prime )
  139.              + e   * 0.002695 * sind( 2*m_prime - m )
  140.              +       0.002602 * sind( m_prime - 2*f - 2*d ) 
  141.              + e   * 0.002396 * sind( 2*d - m - 2*m_prime )
  142.              -       0.002349 * sind( m_prime + d )
  143.              + e_2 * 0.002249 * sind( 2*d - 2*m )
  144.              - e   * 0.002125 * sind( 2*m_prime + m )
  145.              - e_2 * 0.002079 * sind( 2*m )
  146.              + e_2 * 0.002059 * sind( 2*d - m_prime - 2*m )
  147.              -       0.001773 * sind( m_prime + 2*d - 2*f )
  148.              -       0.001595 * sind( 2*f + 2*d );
  149.  
  150.     lambda +=  e   * 0.001220 * sind( 4*d - m - m_prime )
  151.              -       0.001110 * sind( 2*m_prime + 2*f )
  152.              +       0.000892 * sind( m_prime - 3*d )
  153.              - e   * 0.000811 * sind( m + m_prime + 2*d )
  154.              + e   * 0.000761 * sind( 4*d - m - 2*m_prime )
  155.              + e_2 * 0.000717 * sind( m_prime - 2*m )
  156.              + e_2 * 0.000704 * sind( m_prime - 2*m -2*d )
  157.              + e   * 0.000693 * sind( m - 2*m_prime + 2*d )
  158.              + e   * 0.000598 * sind( 2*d - m - 2*f )
  159.              +       0.000550 * sind( m_prime + 4*d )
  160.              +       0.000538 * sind( 4*m_prime )
  161.              + e   * 0.000521 * sind( 4*d - m )
  162.              +       0.000486 * sind( 2*m_prime - d );
  163.  
  164.     lambda = range( lambda );
  165.  
  166.     /* Calculate beta, the ecliptical latitude of the Moon's center. */
  167.     b =         5.128189 * sind( f )
  168.     +       0.280606 * sind( m_prime + f )
  169.     +       0.277693 * sind( m_prime - f )
  170.     +       0.173238 * sind( 2*d - f )
  171.     +       0.055413 * sind( 2*d + f - m_prime )
  172.     +       0.046272 * sind( 2*d - f - m_prime )
  173.     +       0.032573 * sind( 2*d + f )
  174.     +       0.017198 * sind( 2*m_prime + f )
  175.     +       0.009267 * sind( 2*d + m_prime - f )
  176.     +       0.008823 * sind( 2*m_prime - f )
  177.     + e   * 0.008247 * sind( 2*d - m - f )
  178.     +       0.004323 * sind( 2*d - f - 2*m_prime );
  179.  
  180.     b +=    0.004200 * sind( 2*d + f + m_prime )
  181.     + e   * 0.003372 * sind( f - m - 2*d )
  182.     + e   * 0.002472 * sind( 2*d + f - m - m_prime )
  183.     + e   * 0.002222 * sind( 2*d + f - m )
  184.     + e   * 0.002072 * sind( 2*d - f - m - m_prime )
  185.     + e   * 0.001877 * sind( f - m + m_prime )
  186.     +       0.001828 * sind( 4*d - f - m_prime )
  187.     - e   * 0.001803 * sind( f + m )
  188.     -       0.001750 * sind( 3*f )
  189.     + e   * 0.001570 * sind( m_prime - m - f )
  190.     -       0.001487 * sind( f + d )
  191.     - e   * 0.001481 * sind( f + m + m_prime );
  192.  
  193.     b += e* 0.001417 * sind( f - m - m_prime )
  194.     + e   * 0.001350 * sind( f - m )
  195.     +       0.001330 * sind( f - d )
  196.     +       0.001106 * sind( f + 3*m_prime )
  197.     +       0.001020 * sind( 4*d - f )
  198.     +       0.000833 * sind( f + 4*d - m_prime )
  199.     +       0.000781 * sind( m_prime - 3*f )
  200.     +       0.000670 * sind( f + 4*d - 2*m_prime )
  201.     +       0.000606 * sind( 2*d - 3*f )
  202.     +       0.000597 * sind( 2*d + 2*m_prime - f )
  203.     + e   * 0.000492 * sind( 2*d + m_prime - m - f );
  204.  
  205.     b +=    0.000450 * sind( 2*m_prime - f - 2*d )
  206.     +       0.000439 * sind( 3*m_prime - f )
  207.     +       0.000423 * sind( f + 2*d + 2*m_prime )
  208.     +       0.000422 * sind( 2*d - f - 3*m_prime )
  209.     - e   * 0.000367 * sind( m + f + 2*d - m_prime )
  210.     - e   * 0.000353 * sind( m + f + 2*d )
  211.     +       0.000331 * sind( f + 4*d )
  212.     + e   * 0.000317 * sind( 2*d + f - m + m_prime )
  213.     + e_2 * 0.000306 * sind( 2*d - 2*m - f )
  214.     -       0.000283 * sind( m_prime + 3*f );
  215.  
  216.     omega1 = 0.0004664 * cosd( omega );
  217.     omega2 = 0.0000754 * cosd( omega + 275.05 - 2.30 * t );
  218.  
  219.     beta = b * ( 1 - omega1 - omega2 );
  220. /* roughly calculate the phase of the moon.
  221. Added by ccount
  222. Sun longi is about 36000 * T0 + 279
  223. Phase is 0 for full moon */
  224.  
  225.     tmp = 36000.0 * t +279 - lambda;
  226.     phase = 180.0 - (tmp - ((double)((int)(tmp/360.0))*360.0));
  227.  
  228.     phase = (phase < 0.0) ? 360 + phase:phase;
  229.     sprintf(namestring, "Moon, %3.0f", phase);
  230.  
  231.     /* Transform to right ascension and declination, and output the
  232.      * result.
  233.      */
  234.     geo_trans( lambda, beta, epli, 0.0, MAGMOON, "PL", namestring );
  235. }
  236.  
  237.