home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part21 / nutation.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  2.1 KB  |  89 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4. #include "preferences.h"
  5.  
  6. /* given the modified JD, mjd, find the nutation in obliquity, *deps, and
  7.  * the nutation in longitude, *dpsi, each in radians.
  8.  */
  9. void
  10. nutation (mjd, deps, dpsi)
  11. double mjd;
  12. double *deps, *dpsi;
  13. {
  14.     static double lastmjd = -10000, lastdeps, lastdpsi;
  15.     double ls, ld;    /* sun's mean longitude, moon's mean longitude */
  16.     double ms, md;    /* sun's mean anomaly, moon's mean anomaly */
  17.     double nm;    /* longitude of moon's ascending node */
  18.     double t, t2;    /* number of Julian centuries of 36525 days since
  19.              * Jan 0.5 1900.
  20.              */
  21.     double tls, tnm, tld;    /* twice above */
  22.     double a, b;    /* temps */
  23.  
  24.     if (pref_get(PREF_ALGO) == PREF_FAST) {
  25.         *deps = 0.0;
  26.         *dpsi = 0.0;
  27.         return;
  28.     }
  29.  
  30.     if (mjd == lastmjd) {
  31.         *deps = lastdeps;
  32.         *dpsi = lastdpsi;
  33.         return;
  34.     }
  35.         
  36.     t = mjd/36525.;
  37.     t2 = t*t;
  38.  
  39.     a = 100.0021358*t;
  40.     b = 360.*(a-(long)a);
  41.     ls = 279.697+.000303*t2+b;
  42.  
  43.     a = 1336.855231*t;
  44.     b = 360.*(a-(long)a);
  45.     ld = 270.434-.001133*t2+b;
  46.  
  47.     a = 99.99736056000026*t;
  48.     b = 360.*(a-(long)a);
  49.     ms = 358.476-.00015*t2+b;
  50.  
  51.     a = 13255523.59*t;
  52.     b = 360.*(a-(long)a);
  53.     md = 296.105+.009192*t2+b;
  54.  
  55.     a = 5.372616667*t;
  56.     b = 360.*(a-(long)a);
  57.     nm = 259.183+.002078*t2-b;
  58.  
  59.     /* convert to radian forms for use with trig functions.
  60.      */
  61.     tls = 2*degrad(ls);
  62.     nm = degrad(nm);
  63.     tnm = 2*nm;
  64.     ms = degrad(ms);
  65.     tld = 2*degrad(ld);
  66.     md = degrad(md);
  67.  
  68.     /* find delta psi and eps, in arcseconds.
  69.      */
  70.     lastdpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls)
  71.            +.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms)
  72.            +.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms)
  73.            -.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms)
  74.            -.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md);
  75.     lastdeps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls)
  76.            -.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms)
  77.            +.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms)
  78.            -.0066*cos(tls-nm);
  79.  
  80.     /* convert to radians.
  81.      */
  82.     lastdpsi = degrad(lastdpsi/3600);
  83.     lastdeps = degrad(lastdeps/3600);
  84.  
  85.     lastmjd = mjd;
  86.     *deps = lastdeps;
  87.     *dpsi = lastdpsi;
  88. }
  89.