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