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