home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume14 / ephem-4.21 / part02 / moon.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-08-30  |  5.0 KB  |  136 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. /* given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
  6.  * bet, and horizontal parallax, hp for the moon.
  7.  * N.B. series for long and lat are good to about 10 and 3 arcseconds. however,
  8.  *   math errors cause up to 100 and 30 arcseconds error, even if use double.
  9.  *   why?? suspect highly sensitive nature of difference used to get m1..6.
  10.  * N.B. still need to correct for nutation. then for topocentric location
  11.  *   further correct for parallax and refraction.
  12.  */
  13. moon (mjd, lam, bet, hp)
  14. double mjd;
  15. double *lam, *bet, *hp;
  16. {
  17.     double t, t2;
  18.     double ld;
  19.     double ms;
  20.     double md;
  21.     double de;
  22.     double f;
  23.     double n;
  24.     double a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2;
  25.     double m1, m2, m3, m4, m5, m6;
  26.  
  27.     t = mjd/36525.;
  28.     t2 = t*t;
  29.  
  30.     m1 = mjd/27.32158213;
  31.     m1 = 360.0*(m1-(long)m1);
  32.     m2 = mjd/365.2596407;
  33.     m2 = 360.0*(m2-(long)m2);
  34.     m3 = mjd/27.55455094;
  35.     m3 = 360.0*(m3-(long)m3);
  36.     m4 = mjd/29.53058868;
  37.     m4 = 360.0*(m4-(long)m4);
  38.     m5 = mjd/27.21222039;
  39.     m5 = 360.0*(m5-(long)m5);
  40.     m6 = mjd/6798.363307;
  41.     m6 = 360.0*(m6-(long)m6);
  42.  
  43.     ld = 270.434164+m1-(.001133-.0000019*t)*t2;
  44.     ms = 358.475833+m2-(.00015+.0000033*t)*t2;
  45.     md = 296.104608+m3+(.009192+.0000144*t)*t2;
  46.     de = 350.737486+m4-(.001436-.0000019*t)*t2;
  47.     f = 11.250889+m5-(.003211+.0000003*t)*t2;
  48.     n = 259.183275-m6+(.002078+.000022*t)*t2;
  49.  
  50.     a = degrad(51.2+20.2*t);
  51.     sa = sin(a);
  52.     sn = sin(degrad(n));
  53.     b = 346.56+(132.87-.0091731*t)*t;
  54.     sb = .003964*sin(degrad(b));
  55.     c = degrad(n+275.05-2.3*t);
  56.     sc = sin(c);
  57.     ld = ld+.000233*sa+sb+.001964*sn;
  58.     ms = ms-.001778*sa;
  59.     md = md+.000817*sa+sb+.002541*sn;
  60.     f = f+sb-.024691*sn-.004328*sc;
  61.     de = de+.002011*sa+sb+.001964*sn;
  62.     e = 1-(.002495+7.52e-06*t)*t;
  63.     e2 = e*e;
  64.  
  65.     ld = degrad(ld);
  66.     ms = degrad(ms);
  67.     n = degrad(n);
  68.     de = degrad(de);
  69.     f = degrad(f);
  70.     md = degrad(md);
  71.  
  72.     l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+
  73.         .213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+
  74.         .058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+
  75.         .05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms);
  76.     l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))-
  77.         .012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+
  78.         .010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)-
  79.         e*.006783*sin(2*de+ms);
  80.     l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+
  81.         e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+
  82.         .003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+
  83.         .002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)-
  84.         .002349*sin(md+de);
  85.     l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)-
  86.         e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)-
  87.         .001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+
  88.         e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de);
  89.     l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+
  90.          e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+
  91.          e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+
  92.          e*.000521*sin(4*de-ms)+.000486*sin(2*md-de);
  93.     l = l+e2*.000717*sin(md-2*ms);
  94.     *lam = ld+degrad(l);
  95.     range (lam, 2*PI);
  96.  
  97.     g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+
  98.         .173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+
  99.         .032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+
  100.         .008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f);
  101.     g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+
  102.         e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+
  103.         e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+
  104.         e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)-
  105.         .00175*sin(3*f);
  106.     g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+
  107.          e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+
  108.          .001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+
  109.          .000781*sin(md-3*f)+.00067*sin(f+4*de-2*md);
  110.     g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+
  111.         e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+
  112.         .000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+
  113.         .000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)-
  114.         e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de);
  115.     g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)-
  116.         .000283*sin(md+3*f);
  117.     w1 = .0004664*cos(n);
  118.     w2 = .0000754*cos(c);
  119.     *bet = degrad(g)*(1-w1-w2);
  120.  
  121.     *hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+
  122.           .002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+
  123.           e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)-
  124.           e*.000264*cos(ms+md)-.000198*cos(2*f-md);
  125.     *hp = *hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+
  126.          .000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)-
  127.          e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+
  128.          e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+
  129.          e*.000041*cos(ms+de);
  130.     *hp = *hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)-
  131.          .00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+
  132.          e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+
  133.          e*.000019*cos(4*de-ms-md);
  134.     *hp = degrad(*hp);
  135. }
  136.