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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
#include "astro.h"
void
venus(void)
{
double pturbl, pturbb, pturbr;
double lograd;
double dele, enom, vnom, nd, sl;
double v0, t0, m0, j0, s0;
double lsun, elong, ci, dlong;
/*
* here are the mean orbital elements
*/
ecc = .00682069 - .00004774*capt + 0.091e-6*capt2;
incl = 3.393631 + .0010058*capt - 0.97e-6*capt2;
node = 75.779647 + .89985*capt + .00041*capt2;
argp = 130.163833 + 1.408036*capt - .0009763*capt2;
mrad = .7233316;
anom = 212.603219 + 1.6021301540*eday + .00128605*capt2;
motion = 1.6021687039;
/*
* mean anomalies of perturbing planets
*/
v0 = 212.60 + 1.602130154*eday;
t0 = 358.63 + .985608747*eday;
m0 = 319.74 + 0.524032490*eday;
j0 = 225.43 + .083090842*eday;
s0 = 175.8 + .033459258*eday;
v0 *= radian;
t0 *= radian;
m0 *= radian;
j0 *= radian;
s0 *= radian;
incl *= radian;
node *= radian;
argp *= radian;
anom = fmod(anom, 360.)*radian;
/*
* computation of long period terms affecting the mean anomaly
*/
anom +=
(2.761-0.022*capt)*radsec*sin(
13.*t0 - 8.*v0 + 43.83*radian + 4.52*radian*capt)
+ 0.268*radsec*cos(4.*m0 - 7.*t0 + 3.*v0)
+ 0.019*radsec*sin(4.*m0 - 7.*t0 + 3.*v0)
- 0.208*radsec*sin(s0 + 1.4*radian*capt);
/*
* computation of elliptic orbit
*/
enom = anom + ecc*sin(anom);
do {
dele = (anom - enom + ecc * sin(enom)) /
(1 - ecc*cos(enom));
enom += dele;
} while(fabs(dele) > converge);
vnom = 2*atan2(sqrt((1+ecc)/(1-ecc))*sin(enom/2),
cos(enom/2));
rad = mrad*(1 - ecc*cos(enom));
lambda = vnom + argp;
/*
* perturbations in longitude
*/
icosadd(venfp, vencp);
pturbl = cosadd(4, v0, t0, m0, j0);
pturbl *= radsec;
/*
* perturbations in latidude
*/
pturbb = cosadd(3, v0, t0, j0);
pturbb *= radsec;
/*
* perturbations in log radius vector
*/
pturbr = cosadd(4, v0, t0, m0, j0);
/*
* reduction to the ecliptic
*/
lambda += pturbl;
nd = lambda - node;
lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
sl = sin(incl)*sin(nd);
beta = atan2(sl, pyth(sl)) + pturbb;
lograd = pturbr*2.30258509;
rad *= 1 + lograd;
motion *= radian*mrad*mrad/(rad*rad);
/*
* computation of magnitude
*/
lsun = 99.696678 + 0.9856473354*eday;
lsun *= radian;
elong = lambda - lsun;
ci = (rad - cos(elong))/sqrt(1 + rad*rad - 2*rad*cos(elong));
dlong = atan2(pyth(ci), ci)/radian;
mag = -4 + .01322*dlong + .0000004247*dlong*dlong*dlong;
semi = 8.41;
helio();
geo();
}
|