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