home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part14 / comet.c next >
C/C++ Source or Header  |  1993-05-15  |  3KB  |  96 lines
  1. #include <math.h>
  2. #include "astro.h"
  3.  
  4. #if defined(__STDC__) || defined(__cplusplus)
  5. #define P_(s) s
  6. #else
  7. #define P_(s) ()
  8. #endif
  9.  
  10. extern void range P_((double *v, double r));
  11. extern void sunpos P_((double Mjd, double *lsn, double *rsn));
  12.  
  13. void comet P_((double mjd, double ep, double inc, double ap, double qp, double om, double *lpd, double *psi, double *rp, double *rho, double *lam, double *bet));
  14.  
  15. #undef P_
  16.  
  17.  
  18. /* given a modified Julian date, mjd, and a set of heliocentric parabolic
  19.  * orbital elements referred to the epoch of date (mjd):
  20.  *   ep:   epoch of perihelion,
  21.  *   inc:  inclination,
  22.  *   ap:   argument of perihelion (equals the longitude of perihelion minus the
  23.  *       longitude of ascending node)
  24.  *   qp:   perihelion distance,
  25.  *   om:   longitude of ascending node;
  26.  * find:
  27.  *   lpd:  heliocentric longitude, 
  28.  *   psi:  heliocentric latitude,
  29.  *   rp:   distance from the sun to the planet, 
  30.  *   rho:  distance from the Earth to the planet,
  31.  *   lam:  geocentric ecliptic longitude, 
  32.  *   bet:  geocentric ecliptic latitude,
  33.  *         none are corrected for light time, ie, they are the true values for
  34.  *       the given instant.
  35.  *
  36.  * all angles are in radians, all distances in AU.
  37.  * mutual perturbation corrections with other solar system objects are not
  38.  * applied. corrections for nutation and abberation must be made by the caller.
  39.  * The RA and DEC calculated from the fully-corrected ecliptic coordinates are
  40.  * then the apparent geocentric coordinates. Further corrections can be made,
  41.  * if required, for atmospheric refraction and geocentric parallax.
  42.  */
  43. void
  44. comet (mjd, ep, inc, ap, qp, om, lpd, psi, rp, rho, lam, bet)
  45. double mjd;
  46. double ep, inc, ap, qp, om;
  47. double *lpd, *psi, *rp, *rho, *lam, *bet;
  48. {
  49.     double w, s, s2;
  50.     double l, sl, cl, y;
  51.     double spsi, cpsi;
  52.     double rd, lsn, rsn;
  53.     double lg, re, ll;
  54.     double cll, sll;
  55.     double nu;
  56.  
  57. #define    ERRLMT    0.0001
  58.         w = ((mjd-ep)*3.649116e-02)/(qp*sqrt(qp));
  59.         s = w/3;
  60.     for (;;) {
  61.         double d;
  62.         s2 = s*s;
  63.         d = (s2+3)*s-w;
  64.         if (fabs(d) <= ERRLMT)
  65.         break;
  66.         s = ((2*s*s2)+w)/(3*(s2+1));
  67.     }
  68.  
  69.         nu = 2*atan(s);
  70.     *rp = qp*(1+s2);
  71.     l = nu+ap;
  72.         sl = sin(l);
  73.     cl = cos(l);
  74.     spsi = sl*sin(inc);
  75.         *psi = asin(spsi);
  76.     y = sl*cos(inc);
  77.         *lpd = atan(y/cl)+om;
  78.     cpsi = cos(*psi);
  79.         if (cl<0) *lpd += PI;
  80.     range (lpd, 2*PI);
  81.         rd = *rp * cpsi;
  82.     sunpos (mjd, &lsn, &rsn);
  83.     lg = lsn+PI;
  84.         re = rsn;
  85.     ll = *lpd - lg;
  86.         cll = cos(ll);
  87.     sll = sin(ll);
  88.         *rho = sqrt((re * re)+(*rp * *rp)-(2*re*rd*cll));
  89.         if (rd<re) 
  90.             *lam = atan((-1*rd*sll)/(re-(rd*cll)))+lg+PI;
  91.     else
  92.         *lam = atan((re*sll)/(rd-(re*cll)))+*lpd;
  93.     range (lam, 2*PI);
  94.         *bet = atan((rd*spsi*sin(*lam-*lpd))/(cpsi*re*sll));
  95. }
  96.