home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part21 / moonnf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  1.9 KB  |  82 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. #if defined(__STDC__) || defined(__cplusplus)
  6. #define P_(s) s
  7. #else
  8. #define P_(s) ()
  9. #endif
  10.  
  11. extern void cal_mjd P_((int mn, double dy, int yr, double *Mjd));
  12. extern void mjd_cal P_((double Mjd, int *mn, double *dy, int *yr));
  13.  
  14. void moonnf P_((double mjd, double *mjdn, double *mjdf));
  15. static void m P_((double t, double k, double *mjd));
  16.  
  17. #undef P_
  18.  
  19. #define    unw(w,z)    ((w)-floor((w)/(z))*(z))
  20.  
  21. /* given a modified Julian date, mjd, return the mjd of the new
  22.  * and full moons about then, mjdn and mjdf.
  23.  * TODO: exactly which ones does it find? eg:
  24.  *   5/28/1988 yields 5/15 and 5/31
  25.  *   5/29             6/14     6/29
  26.  */
  27. void
  28. moonnf (mjd, mjdn, mjdf)
  29. double mjd;
  30. double *mjdn, *mjdf;
  31. {
  32.     int mo, yr;
  33.     double dy;
  34.     double mjd0;
  35.     double k, tn, tf, t;
  36.  
  37.     mjd_cal (mjd, &mo, &dy, &yr);
  38.     cal_mjd (1, 0., yr, &mjd0);
  39.     k = (yr-1900+((mjd-mjd0)/365))*12.3685;
  40.     k = floor(k+0.5);
  41.     tn = k/1236.85;
  42.     tf = (k+0.5)/1236.85;
  43.     t = tn;
  44.     m (t, k, mjdn);
  45.     t = tf;
  46.     k += 0.5;
  47.     m (t, k, mjdf);
  48. }
  49.  
  50. static void
  51. m (t, k, mjd)
  52. double t, k;
  53. double *mjd;
  54. {
  55.     double t2, a, a1, b, b1, c, ms, mm, f, ddjd;
  56.  
  57.     t2 = t*t;
  58.     a = 29.53*k;
  59.     c = degrad(166.56+(132.87-9.173e-3*t)*t);
  60.     b = 5.8868e-4*k+(1.178e-4-1.55e-7*t)*t2+3.3e-4*sin(c)+7.5933E-1;
  61.     ms = 359.2242+360*unw(k/1.236886e1,1)-(3.33e-5+3.47e-6*t)*t2;
  62.     mm = 306.0253+360*unw(k/9.330851e-1,1)+(1.07306e-2+1.236e-5*t)*t2;
  63.     f = 21.2964+360*unw(k/9.214926e-1,1)-(1.6528e-3+2.39e-6*t)*t2;
  64.     ms = unw(ms,360);
  65.     mm = unw(mm,360);
  66.     f = unw(f,360);
  67.     ms = degrad(ms);
  68.     mm = degrad(mm);
  69.     f = degrad(f);
  70.     ddjd = (1.734e-1-3.93e-4*t)*sin(ms)+2.1e-3*sin(2*ms)
  71.         -4.068e-1*sin(mm)+1.61e-2*sin(2*mm)-4e-4*sin(3*mm)
  72.         +1.04e-2*sin(2*f)-5.1e-3*sin(ms+mm)-7.4e-3*sin(ms-mm)
  73.         +4e-4*sin(2*f+ms)-4e-4*sin(2*f-ms)-6e-4*sin(2*f+mm)
  74.         +1e-3*sin(2*f-mm)+5e-4*sin(ms+2*mm);
  75.     a1 = (long)a;
  76.     b = b+ddjd+(a-a1);
  77.     b1 = (long)b;
  78.     a = a1+b1;
  79.     b = b-b1;
  80.     *mjd = a + b;
  81. }
  82.