1 #include <math.h>
2
3 #include "astro.h"
4
5 /* given a modified Julian date, mj, and a set of heliocentric parabolic
6 * orbital elements referred to the epoch of date (mj):
7 * ep: epoch of perihelion,
8 * inc: inclination,
9 * ap: argument of perihelion (equals the longitude of perihelion minus the
10 * longitude of ascending node)
11 * qp: perihelion distance,
12 * om: longitude of ascending node;
13 * find:
14 * lpd: heliocentric longitude,
15 * psi: heliocentric latitude,
16 * rp: distance from the sun to the planet,
17 * rho: distance from the Earth to the planet,
18 * lam: geocentric ecliptic longitude,
19 * bet: geocentric ecliptic latitude,
20 * none are corrected for light time, ie, they are the true values for
21 * the given instant.
22 *
23 * all angles are in radians, all distances in AU.
24 * mutual perturbation corrections with other solar system objects are not
25 * applied. corrections for nutation and abberation must be made by the caller.
26 * The RA and DEC calculated from the fully-corrected ecliptic coordinates are
27 * then the apparent geocentric coordinates. Further corrections can be made,
28 * if required, for atmospheric refraction and geocentric parallax.
29 */
30 void
comet(double mj,double ep,double inc,double ap,double qp,double om,double * lpd,double * psi,double * rp,double * rho,double * lam,double * bet)31 comet (double mj, double ep, double inc, double ap, double qp, double om,
32 double *lpd, double *psi, double *rp, double *rho, double *lam, double *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 = ((mj-ep)*3.649116e-02)/(qp*sqrt(qp));
44 s = w/3;
45 for (;;) {
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 (mj, &lsn, &rsn, 0);
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
82