home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part21 / sun.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  2.2 KB  |  80 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 anomaly P_((double ma, double s, double *nu, double *ea));
  12. extern void range P_((double *v, double r));
  13.  
  14. void sunpos P_((double mjd, double *lsn, double *rsn));
  15.  
  16. #undef P_
  17.  
  18.  
  19. /* given the modified JD, mjd, return the true geocentric ecliptic longitude
  20.  *   of the sun for the mean equinox of the date, *lsn, in radians, and the
  21.  *   sun-earth distance, *rsn, in AU. (the true ecliptic latitude is never more
  22.  *   than 1.2 arc seconds and so may be taken to be a constant 0.)
  23.  * if the APPARENT ecliptic longitude is required, correct the longitude for
  24.  *   nutation to the true equinox of date and for aberration (light travel time,
  25.  *   approximately  -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
  26.  */
  27. void
  28. sunpos (mjd, lsn, rsn)
  29. double mjd;
  30. double *lsn, *rsn;
  31. {
  32.     static double last_mjd = -3691, last_lsn, last_rsn;
  33.     double t, t2;
  34.     double ls, ms;    /* mean longitude and mean anomoay */
  35.     double s, nu, ea; /* eccentricity, true anomaly, eccentric anomaly */
  36.     double a, b, a1, b1, c1, d1, e1, h1, dl, dr;
  37.  
  38.     if (mjd == last_mjd) {
  39.         *lsn = last_lsn;
  40.         *rsn = last_rsn;
  41.         return;
  42.     }
  43.  
  44.     t = mjd/36525.;
  45.     t2 = t*t;
  46.     a = 100.0021359*t;
  47.     b = 360.*(a-(long)a);
  48.     ls = 279.69668+.0003025*t2+b;
  49.     a = 99.99736042000039*t;
  50.     b = 360*(a-(long)a);
  51.     ms = 358.47583-(.00015+.0000033*t)*t2+b;
  52.     s = .016751-.0000418*t-1.26e-07*t2;
  53.     anomaly (degrad(ms), s, &nu, &ea);
  54.     a = 62.55209472000015*t;
  55.     b = 360*(a-(long)a);
  56.     a1 = degrad(153.23+b);
  57.     a = 125.1041894*t;
  58.     b = 360*(a-(long)a);
  59.     b1 = degrad(216.57+b);
  60.     a = 91.56766028*t;
  61.     b = 360*(a-(long)a);
  62.     c1 = degrad(312.69+b);
  63.     a = 1236.853095*t;
  64.     b = 360*(a-(long)a);
  65.     d1 = degrad(350.74-.00144*t2+b);
  66.     e1 = degrad(231.19+20.2*t);
  67.     a = 183.1353208*t;
  68.     b = 360*(a-(long)a);
  69.     h1 = degrad(353.4+b);
  70.     dl = .00134*cos(a1)+.00154*cos(b1)+.002*cos(c1)+.00179*sin(d1)+
  71.                                 .00178*sin(e1);
  72.     dr = 5.43e-06*sin(a1)+1.575e-05*sin(b1)+1.627e-05*sin(c1)+
  73.                         3.076e-05*cos(d1)+9.27e-06*sin(h1);
  74.     *lsn = nu+degrad(ls-ms+dl);
  75.     range (lsn, 2*PI);
  76.     last_lsn = *lsn;
  77.     *rsn = last_rsn = 1.0000002*(1-s*cos(ea))+dr;
  78.     last_mjd = mjd;
  79. }
  80.