1 // SPDX-License-Identifier: LGPL-2.1-or-later
2 //
3 // SPDX-FileCopyrightText: 2014 Gerhard Holtkamp
4 //
5
6 /* =========================================================================
7 Procedures for calculating the positions of planets.
8 The procedures in this unit are taken from Montenbruck, Pfleger
9 "Astronomie mit dem Personal Computer", Springer Verlag, 1989
10 and modified correspondingly.
11 The library Astrolib has to be included.
12
13 Copyright :Gerhard HOLTKAMP 26-MAR-2014
14 ========================================================================= */
15
16 #include "astr2lib.h"
17 #include <cmath>
18 using namespace std;
19
20 #include "astrolib.h"
21
22 extern double frac (double f);
23 extern double atan21 (double y, double x);
24
25
26 /*-------------------- Class Plan200 --------------------------------------*/
27 /*
28 Ecliptic coordinates (in A.U.) and velocity (in A.U./day)
29 of the Planets for Equinox of Date given in Julian centuries since J2000.
30 ======================================================================
31 */
Plan200()32 Plan200::Plan200 ()
33 { }
34
addthe(double c1,double s1,double c2,double s2,double & cc,double & ss)35 void Plan200::addthe (double c1, double s1, double c2, double s2, double& cc, double& ss)
36 {
37 cc=c1*c2-s1*s2;
38 ss=s1*c2+c1*s2;
39 }
40
term(int i1,int i,int it,double dlc,double dls,double drc,double drs,double dbc,double dbs)41 void Plan200::term (int i1, int i, int it, double dlc, double dls, double drc,
42 double drs, double dbc, double dbs)
43 {
44 if (it == 0) addthe (c3[i1],s3[i1],c[i],s[i],u,v);
45 else
46 {
47 u=u*tt;
48 v=v*tt;
49 }
50 dl = dl + dlc*u + dls*v;
51 dr = dr + drc*u + drs*v;
52 db = db + dbc*u + dbs*v;
53 }
54
velocity()55 Vec3 Plan200::velocity() // return last calculated planet velocity
56 {
57 return vp;
58 }
59
state(Vec3 & rs,Vec3 & vs)60 void Plan200::state (Vec3& rs, Vec3& vs)
61 {
62 /* State vector rs (position) and vs (velocity) of the Sun in
63 ecliptic of date coordinates for last calculated planet
64 */
65 rs = rp;
66 vs = vp;
67 }
68
posvel()69 void Plan200::posvel ()
70 {
71 /* auxiliary program to calculate position and velocity
72 vectors rp, vp.
73 to be called by the various planet procedures.
74 */
75
76 double cl, sl, cb, sb;
77
78 cl = cos(l); sl = sin(l); cb = cos(b); sb = sin(b);
79 rp[0] = r*cl*cb; rp[1] = r*sl*cb; rp[2] = r*sb;
80
81 // velocity
82 vp[0] = (dr*cl*cb - dl*r*sl*cb - db*r*cl*sb) * 1E-4;
83 vp[1] = (dr*sl*cb + dl*r*cl*cb - db*r*sl*sb) * 1E-4;
84 vp[2] = (dr*sb + db*r*cb) * 1E-4;
85 }
86
87 /*--------------------- Mercury ------------------------------------------*/
88
Mercury(double t)89 Vec3 Plan200::Mercury (double t) // position of Mercury at time t
90 {
91 /*
92 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
93 of Mercury for Equinox of Date.
94 t is the time in Julian centuries since J2000.
95 */
96
97 const double arc = 206264.81; // 3600*180/pi = ''/r
98 const double p2 = M_PI * 2.0;
99
100 int i;
101
102 tt = t;
103 dl = 0.0; dr = 0.0; db = 0.0;
104 m1 = p2 * frac(0.4855407 + 415.2014314*t);
105 m2 = p2 * frac(0.1394222 + 162.5490444*t);
106 m3 = p2 * frac(0.9937861 + 99.9978139*t);
107 m5 = p2 * frac(0.0558417 + 8.4298417*t);
108 m6 = p2 * frac(0.8823333 + 3.3943333*t);
109 c3[1] = 1.0; s3[1]= 0.0;
110 c3[2] = cos(m1); s3[2] = sin(m1); c3[0] = c3[2]; s3[0] = -s3[2];
111 for (i=3; i<11; i++)
112 addthe(c3[i-1],s3[i-1],c3[2],s3[2],c3[i],s3[i]);
113
114 // perturbations by Venus
115 c[5] = 1.0; s[5] = 0.0; c[4] = cos(m2); s[4] = -sin(m2);
116 for (i=4; i>0; i--)
117 addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
118 term(2, 5, 0, 259.74, 84547.39, -78342.34,0.01,11683.22,21203.79);
119 term(2, 5, 1, 2.3, 5.04, -7.52, 0.02, 138.55, -71.01);
120 term(2, 5, 2, 0.01, -0.01, 0.01, 0.01, -0.19, -0.54);
121 term(3, 5, 0, -549.71, 10394.44, -7955.45, 0.0, 2390.29, 4306.79);
122 term(3, 5, 1, -4.77, 8.97, -1.53, 0.0, 28.49, -14.18);
123 term(3, 5, 2, 0.0, 0.0, 0.0, 0.0, -0.04, -0.11);
124 term(4, 5, 0, -234.04, 1748.74, -1212.86, 0.0, 535.41, 984.33);
125 term(4, 5, 1, -2.03, 3.48, -0.35, 0.0, 6.56, -2.91);
126 term(5, 5, 0, -77.64, 332.63, -219.23, 0.0, 124.4, 237.03);
127 term(5, 5, 1, -0.7, 1.1, -0.08, 0.0, 1.59, -0.59);
128 term(6, 5, 0, -23.59, 67.28, -43.54, 0.0, 29.44, 58.77);
129 term(6, 5, 1, -0.23, 0.32, -0.02, 0.0, 0.39, -0.11);
130 term(7, 5, 0, -6.86, 14.06, -9.18, 0.0, 7.03, 14.84);
131 term(7, 5, 1, -0.07, 0.09, -0.01, 0.0, 0.1, -0.02);
132 term(8, 5, 0, -1.94, 2.98, -2.02, 0.0, 1.69, 3.8);
133 term(9, 5, 0, -0.54, 0.63, -0.46, 0.0, 0.41, 0.98);
134 term(10, 5, 0, -0.15, 0.13, -0.11, 0.0, 0.1, 0.25);
135 term(0, 3, 0, -0.17, -0.06, -0.05, 0.14, -0.06, -0.07);
136 term(1, 4, 0, 0.24, -0.16, -0.11, -0.16, 0.04, -0.01);
137 term(1, 3, 0, -0.68, -0.25, -0.26, 0.73, -0.16, -0.18);
138 term(1, 0, 0, 0.37, 0.08, 0.06, -0.28, 0.13, 0.12);
139 term(2, 4, 0, 0.58, -0.41, 0.26, 0.36, 0.01, -0.01);
140 term(2, 3, 0, -3.51, -1.23, 0.23, -0.63, -0.05, -0.06);
141 term(2, 2, 0, 0.08, 0.53, -0.11, 0.04, 0.02, -0.09);
142 term(2, 0, 0, 1.44, 0.31, 0.3, -1.39, 0.34, 0.29);
143 term(3, 4, 0, 0.15, -0.11, 0.09, 0.12, 0.02, -0.04);
144 term(3, 3, 0, -1.99, -0.68, 0.65, -1.91, -0.2, 0.03);
145 term(3, 2, 0, -0.34, -1.28, 0.97, -0.26, 0.03, 0.03);
146 term(3, 1, 0, -0.33, 0.35, -0.13, -0.13, -0.01, 0.0);
147 term(3, 0, 0, 7.19, 1.56, -0.05, 0.12, 0.06, 0.05);
148 term(4, 3, 0, -0.52, -0.18, 0.13, -0.39, -0.16, 0.03);
149 term(4, 2, 0, -0.11, -0.42, 0.36, -0.1, -0.05, -0.05);
150 term(4, 1, 0, -0.19, 0.22, -0.23, -0.2, -0.01, 0.02);
151 term(4, 0, 0, 2.77, 0.49, -0.45, 2.56, 0.4, -0.12);
152 term(5, 0, 0, 0.67, 0.12, -0.09, 0.47, 0.24, -0.08);
153 term(6, 0, 0, 0.18, 0.03, -0.02, 0.12, 0.09, -0.03);
154
155 // perturbations by Earth
156 c[4] = cos(m3); s[4] = -sin(m3);
157 for (i=4; i>1; i--)
158 addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
159 term(1, 1, 0, -0.11, -0.07, -0.08, 0.11, -0.02, -0.04);
160 term(2, 4, 0, 0.1, -0.2, 0.15, 0.07, 0.0, 0.0);
161 term(2, 3, 0, -0.35, 0.28, -0.13, -0.17, -0.01, 0.0);
162 term(2, 1, 0, -0.67, -0.45, 0.0, 0.01, -0.01, -0.01);
163 term(3, 3, 0, -0.2, 0.16, -0.16, -0.2, -0.01, 0.02);
164 term(3, 2, 0, 0.13, -0.02, 0.02, 0.14, 0.01, 0.0);
165 term(3, 1, 0, -0.33, -0.18, 0.17, -0.31, -0.04, 0.0);
166
167 // perturbations by Jupiter
168 c[4] = cos(m5); s[4] = -sin(m5);
169 for (i=4; i>2; i--)
170 addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
171 term(0, 4, 0, -0.08, 0.16, 0.15, 0.08, -0.04, 0.01);
172 term(0, 3, 0, 0.1, -0.06, -0.07, -0.12, 0.07, -0.01);
173 term(1, 4, 0, -0.31, 0.48, -0.02, 0.13, -0.03, -0.02);
174 term(1, 3, 0, 0.42, -0.26, -0.38, -0.5, 0.2, -0.03);
175 term(2, 4, 0, -0.7, 0.01, -0.02, -0.63, 0.0, 0.03);
176 term(2, 3, 0, 2.61, -1.97, 1.74, 2.32, 0.01, 0.01);
177 term(2, 2, 0, 0.32, -0.15, 0.13, 0.28, 0.0, 0.0);
178 term(3, 4, 0, -0.18, 0.01, 0.0, -0.13, -0.03, 0.03);
179 term(3, 3, 0, 0.75, -0.56, 0.45, 0.6, 0.08, -0.17);
180 term(4, 3, 0, 0.2, -0.15, 0.1, 0.14, 0.04, -0.08);
181
182 // perturbations by Saturn
183 c[3] = cos(2*m6); s[3] = -sin(2*m6);
184 term(2, 3, 0, -0.19, 0.33, 0.0, 0.0, 0.0, 0.0);
185
186 dl = dl + (2.8 + 3.2*t);
187 l = p2 * frac(0.2151379 + m1/p2 + ((5601.7+1.1*t)*t+dl)/1296.0e3);
188 b = (-2522.15+(-30.18+0.04*t)*t+db) / arc;
189 r = 0.3952829 +0.0000016*t + dr*1.0e-6;
190 dl = 714.0+292.66*cos(m1)+71.96*cos(2*m1)+18.16*cos(3*m1)
191 +4.61*cos(4*m1)+3.81*sin(2*m1)+2.43*sin(3*m1)+1.08*sin(4*m1);
192 dr = 55.94*sin(m1)+11.36*sin(2*m1)+2.6*sin(3*m1);
193 db = 73.4*cos(m1)+29.82*cos(2*m1)+10.22*cos(3*m1)+3.28*cos(4*m1)
194 -40.44*sin(m1)-16.55*sin(2*m1)-5.56*sin(3*m1)-1.72*sin(4*m1);
195
196 posvel();
197
198 return rp;
199 }
200
201 /*--------------------- Venus ------------------------------------------*/
202
Venus(double t)203 Vec3 Plan200::Venus (double t) // position of Venus at time t
204 {
205 /*
206 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
207 of Venus for Equinox of Date.
208 t is the time in Julian centuries since J2000.
209 */
210
211 const double arc = 206264.81; // 3600*180/pi = ''/r
212 const double p2 = M_PI * 2.0;
213
214 int i;
215
216 tt = t;
217 dl = 0.0; dr = 0.0; db = 0.0;
218 m1 = p2 * frac(0.4861431+415.2018375*t);
219 m2 = p2 * frac(0.1400197+162.5494552*t);
220 m3 = p2 * frac(0.9944153+99.9982208*t);
221 m4 = p2 * frac(0.0556297+53.1674631*t);
222 m5 = p2 * frac(0.0567028+8.4305083*t);
223 m6 = p2 * frac(0.8830539+3.3947206*t);
224
225 c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m2); s3[1] = sin(m2);
226 for (i=2; i<9; i++)
227 addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);
228
229 // perturbations by Mercury
230 c[8] = 1.0; s[8] = 0.0; c[7] = cos(m1); s[7] = -sin(m1);
231 addthe(c[7],s[7],c[7],s[7],c[6],s[6]);
232 term(1, 7, 0, 0.0, 0.0, 0.06, -0.09, 0.01, 0.0);
233 term(2, 7, 0, 0.25, -0.09, -0.09, -0.27, 0.0, 0.0);
234 term(4, 6, 0, -0.07, -0.08, -0.14, 0.14, -0.01, -0.01);
235 term(5, 6, 0, -0.35, 0.08, 0.02, 0.09, 0.0, 0.0);
236
237 // perturbations by Earth
238 c[7] = cos(m3); s[7] = -sin(m3);
239 for (i=7; i>0; i--)
240 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
241 term(1, 8, 0, 2.37, 2793.23, -4899.07, 0.11, 9995.27, 7027.22);
242 term(1, 8, 1, 0.1, -19.65, 34.4, 0.22, 64.95, -86.1);
243 term(1, 8, 2, 0.06, 0.04, -0.07, 0.11, -0.55, -0.07);
244 term(2, 8, 0, -170.42, 73.13, -16.59, 0.0, 67.71, 47.56);
245 term(2, 8, 1, 0.93, 2.91, 0.23, 0.0, -0.03, -0.92);
246 term(3, 8, 0, -2.31, 0.9, -0.08, 0.0, 0.04, 2.09);
247 term(1, 7, 0, -2.38, -4.27, 3.27, -1.82, 0.0, 0.0);
248 term(1, 6, 0, 0.09, 0.0, -0.08, 0.05, -0.02, -0.25);
249 term(2, 6, 0, -9.57, -5.93, 8.57, -13.83, -0.01, -0.01);
250 term(2, 5, 0, -2.47, -2.4, 0.83, -0.95, 0.16, 0.24);
251 term(3, 6, 0, -0.09, -0.05, 0.08, -0.13, -0.28, 0.12);
252 term(3, 5, 0, 7.12, 0.32, -0.62, 13.76, -0.07, 0.01);
253 term(3, 4, 0, -0.65, -0.17, 0.18, -0.73, 0.1, 0.05);
254 term(3, 3, 0, -1.08, -0.95, -0.17, 0.22, -0.03, -0.03);
255 term(4, 5, 0, 0.06, 0.0, -0.01, 0.08, 0.14, -0.18);
256 term(4, 4, 0, 0.93, -0.46, 1.06, 2.13, -0.01, 0.01);
257 term(4, 3, 0, -1.53, 0.38, -0.64, -2.54, 0.27, 0.0);
258 term(4, 2, 0, -0.17, -0.05, 0.03, -0.11, 0.02, 0.0);
259 term(5, 3, 0, 0.18, -0.28, 0.71, 0.47, -0.02, 0.04);
260 term(5, 2, 0, 0.15, -0.14, 0.3, 0.31, -0.04, 0.03);
261 term(5, 1, 0, -0.08, 0.02, -0.03, -0.11, 0.01, 0.0);
262 term(5, 0, 0, -0.23, 0.0, 0.01, -0.04, 0.0, 0.0);
263 term(6, 2, 0, 0.01, -0.14, 0.39, 0.04, 0.0, -0.01);
264 term(6, 1, 0, 0.02, -0.05, 0.12, 0.04, -0.01, 0.01);
265 term(6, 0, 0, 0.1, -0.1, 0.19, 0.19, -0.02, 0.02);
266 term(7, 1, 0, -0.03, -0.06, 0.18, -0.08, 0.0, 0.0);
267 term(8, 0, 0, -0.03, -0.02, 0.06, -0.08, 0.0, 0.0);
268
269 // perturbations by Mars
270 c[7] = cos(m4); s[7] = -sin(m4);
271 for (i=7; i>5; i--)
272 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
273 term(1, 5, 0, -0.65, 1.02, -0.04, -0.02, -0.02, 0.0);
274 term(2, 6, 0, -0.05, 0.04, -0.09, -0.1, 0.0, 0.0);
275 term(2, 5, 0, -0.5, 0.45, -0.79, -0.89, 0.01, 0.03);
276
277 // perturbations by Jupiter
278 c[7] = cos(m5); s[7] = -sin(m5);
279 for (i=7; i>5; i--)
280 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
281 term(0, 7, 0, -0.05, 1.56, 0.16, 0.04, -0.08, -0.04);
282 term(1, 7, 0, -2.62, 1.4, -2.35, -4.4, 0.02, 0.03);
283 term(1, 6, 0, -0.47, -0.08, 0.12, -0.76, 0.04, -0.18);
284 term(2, 6, 0, -0.73, -0.51, 1.27, -1.82, -0.01, 0.01);
285 term(2, 5, 0, -0.14, -0.1, 0.25, -0.34, 0.0, 0.0);
286 term(3, 5, 0, -0.01, 0.04, -0.11, -0.02, 0.0, 0.0);
287
288 // perturbations by Saturn
289 c[7] = cos(m6); s[7] = -sin(m6);
290 term(0, 7, 0, 0.0, 0.21, 0.0, 0.0, 0.0, -0.01);
291 term(1, 7, 0, -0.11, -0.14, 0.24, -0.2, 0.01, 0.0);
292
293 dl = dl+2.74*sin(p2*(0.0764+0.4174*t))+0.27*sin(p2*(0.9201+0.3307*t));
294 dl = dl+(1.9+1.8*t);
295
296 l = p2 * frac(0.3654783 + m2/p2 + ((5071.2+1.1*t)*t+dl)/1296.0e3);
297 b = (-67.7+(0.04+0.01*t)*t+db)/arc;
298 r = 0.7233482-0.0000002*t+dr*1.0e-6;
299
300 dl = 280.0+3.79*cos(m2);
301 dr = 1.37*sin(m2);
302 db = 9.54*cos(m2)-13.57*sin(m2);
303
304 posvel();
305
306 return rp;
307 }
308
309 /*--------------------- Mars ------------------------------------------*/
310
Mars(double t)311 Vec3 Plan200::Mars (double t) // position of Mars at time t
312 {
313 /*
314 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
315 of Mars for Equinox of Date.
316 t is the time in Julian centuries since J2000.
317 */
318
319 const double arc = 206264.81; // 3600*180/pi = ''/r
320 const double p2 = M_PI * 2.0;
321
322 int i;
323
324 tt = t;
325 dl = 0.0; dr = 0.0; db = 0.0;
326 m2 = p2 * frac(0.1382208+162.5482542*t);
327 m3 = p2 * frac(0.9926208+99.9970236*t);
328 m4 = p2 * frac(0.0538553+53.1662736*t);
329 m5 = p2 * frac(0.0548944+8.4290611*t);
330 m6 = p2 * frac(0.8811167+3.3935250*t);
331 c3[2] = 1.0; s3[2] = 0.0; c3[3] = cos(m4); s3[3] = sin(m4);
332 for (i=4; i<19; i++)
333 addthe(c3[i-1],s3[i-1],c3[3],s3[3],c3[i],s3[i]);
334 c3[0] = c3[4];
335 s3[0] = -s3[4];
336 c3[1] = c3[3];
337 s3[1] = -s3[3];
338
339 // perturbations by Venus
340 c[9] = 1.0; s[9] = 0.0; c[8] = cos(m2); s[8] = -sin(m2);
341 addthe (c[8],s[8],c[8],s[8],c[7],s[7]);
342 term(2, 8, 0, -0.01, -0.03, 0.1, -0.04, 0.0, 0.0);
343 term(3, 8, 0, 0.05, 0.1, -2.08, 0.75, 0.0, 0.0);
344 term(4, 8, 0, -0.25, -0.57, -2.58, 1.18, 0.05, -0.04);
345 term(4, 7, 0, 0.02, 0.02, 0.13, -0.14, 0.0, 0.0);
346 term(5, 8, 0, 3.41, 5.38, 1.87, -1.15, 0.01, -0.01);
347 term(5, 7, 0, 0.02, 0.02, 0.11, -0.13, 0.0, 0.0);
348 term(6, 8, 0, 0.32, 0.49, -1.88, 1.21, -0.07, 0.07);
349 term(6, 7, 0, 0.03, 0.03, 0.12, -0.14, 0.0, 0.0);
350 term(7, 8, 0, 0.04, 0.06, -0.17, 0.11, -0.01, 0.01);
351 term(7, 7, 0, 0.11, 0.09, 0.35, -0.43, -0.01, 0.01);
352 term(8, 7, 0, -0.36, -0.28, -0.2, 0.25, 0.0, 0.0);
353 term(9, 7, 0, -0.03, -0.03, 0.11, -0.13, 0.0, 0.0);
354
355 // perturbations by Earth
356 c[8] = cos(m3); s[8] = -sin(m3);
357 for (i=8; i>0; i--)
358 addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
359 term(3, 9, 0, -5.32, 38481.97, -141856.04, 0.4, -6321.67, 1876.89);
360 term(3, 9, 1, -1.12, 37.98, -138.67, -2.93, 37.28, 117.48);
361 term(3, 9, 2, -0.32, -0.03, 0.12, -1.19, 1.04, -0.4);
362 term(4, 9, 0, 28.28, 2285.8, -6608.37, 0.0, -589.35, 174.81);
363 term(4, 9, 1, 1.64, 3.37, -12.93, 0.0, 2.89, 11.1);
364 term(4, 9, 2, 0.0, 0.0, 0.0, 0.0, 0.1, -0.03);
365 term(5, 9, 0, 5.31, 189.29, -461.81, 0.0, -61.98, 18.53);
366 term(5, 9, 1, 0.31, 0.35, -1.36, 0.0, 0.25, 1.19);
367 term(6, 9, 0, 0.81, 17.96, -38.26, 0.0, -6.88, 2.08);
368 term(6, 9, 1, 0.05, 0.04, -0.15, 0.0, 0.02, 0.14);
369 term(7, 9, 0, 0.11, 1.83, -3.48, 0.0, -0.79, 0.24);
370 term(8, 9, 0, 0.02, 0.2, -0.34, 0.0, -0.09, 0.3);
371 term(1, 8, 0, 0.09, 0.06, 0.14, -0.22, 0.02, -0.02);
372 term(2, 8, 0, 0.72, 0.49, 1.55, -2.31, 0.12, -0.1);
373 term(3, 8, 0, 7.0, 4.92, 13.93, -20.48, 0.08, -0.13);
374 term(4, 8, 0, 13.08, 4.89, -4.53, 10.01, -0.05, 0.13);
375 term(4, 7, 0, 0.14, 0.05, -0.48, -2.66, 0.01, 0.14);
376 term(5, 8, 0, 1.38, 0.56, -2.0, 4.85, -0.01, 0.19);
377 term(5, 7, 0, -6.85, 2.68, 8.38, 21.42, 0.0, 0.03);
378 term(5, 6, 0, -0.08, 0.2, 1.2, 0.46, 0.0, 0.0);
379 term(6, 8, 0, 0.16, 0.07, -0.19, 0.47, -0.01, 0.05);
380 term(6, 7, 0, -4.41, 2.14, -3.33, -7.21, -0.07, -0.09);
381 term(6, 6, 0, -0.12, 0.33, 2.22, 0.72, -0.03, -0.02);
382 term(6, 5, 0, -0.04, -0.06, -0.36, 0.23, 0.0, 0.0);
383 term(7, 7, 0, -0.44, 0.21, -0.7, -1.46, -0.06, -0.07);
384 term(7, 6, 0, 0.48, -2.6, -7.25, -1.37, 0.0, 0.0);
385 term(7, 5, 0, -0.09, -0.12, -0.66, 0.5, 0.0, 0.0);
386 term(7, 4, 0, 0.03, 0.0, 0.01, -0.17, 0.0, 0.0);
387 term(8, 7, 0, -0.05, 0.03, -0.07, -0.15, -0.01, -0.01);
388 term(8, 6, 0, 0.1, -0.96, 2.36, 0.3, 0.04, 0.0);
389 term(8, 5, 0, -0.17, -0.2, -1.09, 0.94, 0.02, -0.02);
390 term(8, 4, 0, 0.05, 0.0, 0.0, -0.3, 0.0, 0.0);
391 term(9, 6, 0, 0.01, -0.1, 0.32, 0.04, 0.02, 0.0);
392 term(9, 5, 0, 0.86, 0.77, 1.86, -2.01, 0.01, -0.01);
393 term(9, 4, 0, 0.09, -0.01, -0.05, -0.44, 0.0, 0.0);
394 term(9, 3, 0, -0.01, 0.02, 0.1, 0.08, 0.0, 0.0);
395 term(10, 5, 0, 0.2, 0.16, -0.53, 0.64, -0.01, 0.02);
396 term(10, 4, 0, 0.17, -0.03, -0.14, -0.84, 0.0, 0.01);
397 term(10, 3, 0, -0.02, 0.03, 0.16, 0.09, 0.0, 0.0);
398 term(11, 4, 0, -0.55, 0.15, 0.3, 1.1, 0.0, 0.0);
399 term(11, 3, 0, -0.02, 0.04, 0.2, 0.1, 0.0, 0.0);
400 term(12, 4, 0, -0.09, 0.03, -0.1, -0.33, 0.0, -0.01);
401 term(12, 3, 0, -0.05, 0.11, 0.48, 0.21, -0.01, 0.0);
402 term(13, 3, 0, 0.1, -0.35, -0.52, -0.15, 0.0, 0.0);
403 term(13, 2, 0, -0.01, -0.02, -0.1, 0.07, 0.0, 0.0);
404 term(14, 3, 0, 0.01, -0.04, 0.18, 0.4, 0.01, 0.0);
405 term(14, 2, 0, -0.05, -0.07, -0.29, 0.2, 0.01, 0.0);
406 term(15, 2, 0, 0.23, 0.27, 0.25, -0.21, 0.0, 0.0);
407 term(16, 2, 0, 0.02, 0.03, -0.1, 0.09, 0.0, 0.0);
408 term(16, 1, 0, 0.05, 0.01, 0.03, -0.23, 0.0, 0.3);
409 term(17, 1, 0, -1.53, 0.27, 0.06, 0.42, 0.0, 0.0);
410 term(18, 1, 0, -0.14, 0.02, -0.1, -0.55, -0.01, -0.02);
411 term(18, 0, 0, 0.03, -0.06, -0.25, -0.11, 0.0, 0.0);
412
413 // perturbation by Jupiter
414 c[8] = cos(m5); s[8] = -sin(m5);
415 for (i=8; i>4; i--)
416 addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
417 term(0, 8, 0, 0.05, 0.03, 0.08, -0.14, 0.01, -0.01);
418 term(1, 8, 0, 0.39, 0.27, 0.92, -1.5, -0.03, -0.06);
419 term(1, 7, 0, -0.16, 0.03, 0.13, 0.67, -0.01, 0.06);
420 term(1, 6, 0, -0.02, 0.01, 0.05, 0.09, 0.0, 0.01);
421 term(2, 8, 0, 3.56, 1.13, -5.41, -7.18, -0.25, -0.24);
422 term(2, 7, 0, -1.44, 0.25, 1.24, 7.96, 0.02, 0.31);
423 term(2, 6, 0, -0.21, 0.11, 0.55, 1.04, 0.01, 0.05);
424 term(2, 5, 0, -0.02, 0.02, 0.11, 0.11, 0.0, 0.01);
425 term(3, 8, 0, 16.67, -19.15, 61.0, 53.36, -0.06, -0.07);
426 term(3, 7, 0, -21.64, 3.18, -7.77, -54.64, -0.31, 0.5);
427 term(3, 6, 0, -2.82, 1.45, -2.53, -5.73, 0.01, 0.07);
428 term(3, 5, 0, -0.31, 0.28, -0.34, -0.51, 0.0, 0.0);
429 term(4, 8, 0, 2.15, -2.29, 7.04, 6.94, 0.33, 0.19);
430 term(4, 7, 0, -15.69, 3.31, -15.7, -73.17, -0.17, -0.25);
431 term(4, 6, 0, -1.73, 1.95, -9.19, -7.2, 0.02, -0.03);
432 term(4, 5, 0, -0.01, 0.33, -1.42, 0.08, 0.01, -0.01);
433 term(4, 4, 0, 0.03, 0.03, -0.13, 0.12, 0.0, 0.0);
434 term(5, 8, 0, 0.26, -0.28, 0.73, 0.71, 0.08, 0.04);
435 term(5, 7, 0, -2.06, 0.46, -1.61, -6.72, -0.13, -0.25);
436 term(5, 6, 0, -1.28, -0.27, 2.21, -6.9, -0.04, -0.02);
437 term(5, 5, 0, -0.22, 0.08, -0.44, -1.25, 0.0, 0.01);
438 term(5, 4, 0, -0.02, 0.03, -0.15, -0.08, 0.0, 0.0);
439 term(6, 8, 0, 0.03, -0.03, 0.08, 0.08, 0.01, 0.01);
440 term(6, 7, 0, -0.26, 0.06, -0.17, -0.7, -0.03, -0.05);
441 term(6, 6, 0, -0.2, -0.05, 0.22, -0.79, -0.01, -0.02);
442 term(6, 5, 0, -0.11, -0.14, 0.93, -0.6, 0.0, 0.0);
443 term(6, 4, 0, -0.04, -0.02, 0.09, -0.23, 0.0, 0.0);
444 term(7, 5, 0, -0.02, -0.03, 0.13, -0.09, 0.0, 0.0);
445 term(7, 4, 0, 0.0, -0.03, 0.21, 0.01, 0.0, 0.0);
446
447 // perturbations by Saturn
448 c[8] = cos(m6); s[8] = -sin(m6);
449 for (i=8; i>5; i--)
450 addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
451 term(1, 8, 0, 0.03, 0.13, 0.48, -0.13, 0.02, 0.0);
452 term(2, 8, 0, 0.27, 0.84, 0.4, -0.43, 0.01, -0.01);
453 term(2, 7, 0, 0.12, -0.04, -0.33, -0.55, -0.01, -0.02);
454 term(2, 6, 0, 0.02, -0.01, -0.07, -0.08, 0.0, 0.0);
455 term(3, 8, 0, 1.12, 0.76, -2.66, 3.91, -0.01, 0.01);
456 term(3, 7, 0, 1.49, -0.95, 3.07, 4.83, 0.04, -0.05);
457 term(3, 6, 0, 0.21, -0.18, 0.55, 0.64, 0.0, 0.0);
458 term(4, 8, 0, 0.12, 0.1, -0.29, 0.34, -0.01, 0.02);
459 term(4, 7, 0, 0.51, -0.36, 1.61, 2.25, 0.03, 0.01);
460 term(4, 6, 0, 0.1, -0.1, 0.5, 0.43, 0.0, 0.0);
461 term(4, 5, 0, 0.01, -0.02, 0.11, 0.05, 0.0, 0.0);
462 term(5, 7, 0, 0.07, -0.05, 0.16, 0.22, 0.01, 0.01);
463
464 dl = dl+52.49*sin(p2*(0.1868+0.0549*t))+0.61*sin(p2*(0.922+0.3307*t))
465 +0.32*sin(p2*(0.4731+2.1485*t))+0.28*sin(p2*(0.9467+1.1133*t));
466 dl = dl + (0.14+0.87*t-0.11*t*t);
467 l = p2 * frac(0.9334591+m4/p2+((6615.5+1.1*t)*t+dl)/1296.0e3);
468 r = 1.5303352 + 0.0000131*t + dr*1.0e-6;
469 b = (596.32+(-2.92-0.1*t)*t+db) / arc;
470
471 dl = 91.50 + 17.07*cos(m4) + 2.03*cos(2*m4);
472 dr = 12.98*sin(m4) + 1.21*cos(2*m4);
473 db = 0.83*cos(m4) + 2.8*sin(m4);
474
475 posvel();
476
477 return rp;
478 }
479
480 /*--------------------- Jupiter ------------------------------------------*/
481
Jupiter(double t)482 Vec3 Plan200::Jupiter (double t) // position of Jupiter at time t
483 {
484 /*
485 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
486 of Jupiter for Equinox of Date.
487 t is the time in Julian centuries since J2000.
488 */
489
490 const double arc = 206264.81; // 3600*180/pi = ''/r
491 const double p2 = M_PI * 2.0;
492
493 int i;
494
495 tt = t;
496 dl = 0.0; dr = 0.0; db = 0.0;
497 m5 = p2 * frac(0.0565314+8.4302963*t);
498 m6 = p2 * frac(0.8829867+3.3947688*t);
499 m7 = p2 * frac(0.3969537+1.1902586*t);
500 c3[1] = 1.0; s3[1] = 0.0;
501 c3[2] = cos(m5); s3[2] = sin(m5); c3[0] = c3[2]; s3[0] = -s3[2];
502 for (i=3; i<7; i++)
503 addthe(c3[i-1],s3[i-1],c3[2],s3[2],c3[i],s3[i]);
504
505 // perturbations by Saturn
506 c[10] = 1.0; s[10] = 0.0; c[9] = cos(m6); s[9] = -sin(m6);
507 for (i=9; i>0; i--)
508 addthe(c[i],s[i],c[9],s[9],c[i-1],s[i-1]);
509 term(0, 9, 0, -0.2, 1.4, 2.0, 0.6, 0.1, -0.2);
510 term(1, 9, 0, 9.4, 8.9, 3.9, -8.3, -0.4, -1.4);
511 term(1, 8, 0, 5.6, -3.0, -5.4, -5.7, -2.0, 0.0);
512 term(1, 7, 0, -4.0, -0.1, 0.0, 5.5, 0.0, 0.0);
513 term(1, 5, 0, 3.3, -1.6, -1.6, -3.1, -0.5, -1.2);
514 term(2, 10, 0, -113.1, 19998.6, -25208.2, -142.2, -4670.7, 288.9);
515 term(2, 10, 1, -76.1, 66.9, -84.2, -95.8, 21.6, 29.4);
516 term(2, 10, 2, -0.5, -0.3, 0.4, -0.7, 0.1, -0.1);
517 term(2, 9, 0, 78.8, -14.5, 11.5, 64.4, -0.2, 0.2);
518 term(2, 8, 0, -2.0, -132.4, 28.8, 4.3, -1.7, 0.4);
519 term(2, 8, 1, -1.1, -0.7, 0.2, -0.3, 0.0, 0.0);
520 term(2, 7, 0, -7.5, -6.8, -0.4, -1.1, 0.6, -0.9);
521 term(2, 6, 0, 0.7, 0.7, 0.6, -1.1, 0.0, -0.2);
522 term(2, 5, 0, 51.5, -26.0, -32.5, -64.4, -4.9, -12.4);
523 term(2, 5, 1, -1.2, -2.2, -2.7, 1.5, -0.4, 0.3);
524 term(3, 10, 0, -3.4, 632.0, -610.6, -6.5, -226.8, 12.7);
525 term(3, 10, 1, -4.2, 3.8, -4.1, -4.5, 0.2, 0.6);
526 term(3, 9, 0, 5.3, -0.7, 0.7, 6.1, 0.2, 1.1);
527 term(3, 8, 0, -76.4, -185.1, 260.2, -108.0, 1.6, 0.0);
528 term(3, 7, 0, 66.7, 47.8, -51.4, 69.8, 0.9, 0.3);
529 term(3, 7, 1, 0.6, -1.0, 1.0, 0.6, 0.0, 0.0);
530 term(3, 6, 0, 17.0, 1.4, -1.8, 9.6, 0.0, -0.1);
531 term(3, 5, 0, 1066.2, -518.3, -1.3, -23.9, 1.8, -0.3);
532 term(3, 5, 1, -25.4, -40.3, -0.9, 0.3, 0.0, 0.0);
533 term(3, 5, 2, -0.7, 0.5, 0.0, 0.0, 0.0, 0.0);
534 term(4, 10, 0, -0.1, 28.0, -22.1, -0.2, -12.5, 0.7);
535 term(4, 8, 0, -5.0, -11.5, 11.7, -5.4, 2.1, -1.0);
536 term(4, 7, 0, 16.9, -6.4, 13.4, 26.9, -0.5, 0.8);
537 term(4, 6, 0, 7.2, -13.3, 20.9, 10.5, 0.1, -0.1);
538 term(4, 5, 0, 68.5, 134.3, -166.9, 86.5, 7.1, 15.2);
539 term(4, 5, 1, 3.5, -2.7, 3.4, 4.3, 0.5, -0.4);
540 term(4, 4, 0, 0.6, 1.0, -0.9, 0.5, 0.0, 0.0);
541 term(4, 3, 0, -1.1, 1.7, -0.4, -0.2, 0.0, 0.0);
542 term(5, 10, 0, 0.0, 1.4, -1.0, 0.0, -0.6, 0.0);
543 term(5, 8, 0, -0.3, -0.7, 0.4, -0.2, 0.2, -0.1);
544 term(5, 7, 0, 1.1, -0.6, 0.9, 1.2, 0.1, 0.2);
545 term(5, 6, 0, 3.2, 1.7, -4.1, 5.8, 0.2, 0.1);
546 term(5, 5, 0, 6.7, 8.7, -9.3, 8.7, -1.1, 1.6);
547 term(5, 4, 0, 1.5, -0.3, 0.6, 2.4, 0.0, 0.0);
548 term(5, 3, 0, -1.9, 2.3, -3.2, -2.7, 0.0, -0.1);
549 term(5, 2, 0, 0.4, -1.8, 1.9, 0.5, 0.0, 0.0);
550 term(5, 1, 0, -0.2, -0.5, 0.3, -0.1, 0.0, 0.0);
551 term(5, 0, 0, -8.6, -6.8, -0.4, 0.1, 0.0, 0.0);
552 term(5, 0, 1, -0.5, 0.6, 0.0, 0.0, 0.0, 0.0);
553 term(6, 5, 0, -0.1, 1.5, -2.5, -0.8, -0.1, 0.1);
554 term(6, 4, 0, 0.1, 0.8, -1.6, 0.1, 0.0, 0.0);
555 term(6, 1, 0, -0.5, -0.1, 0.1, -0.8, 0.0, 0.0);
556 term(6, 0, 0, 2.5, -2.2, 2.8, 3.1, 0.1, -0.2);
557
558 // perturbations by Uranus
559 c[9] = cos(m7); s[9] = -sin(m7);
560 addthe(c[9],s[9],c[9],s[9],c[8],s[8]);
561 term(2, 9, 0, 0.4, 0.9, 0.0, 0.0, 0.0, 0.0);
562 term(2, 8, 0, 0.4, 0.4, -0.4, 0.3, 0.0, 0.0);
563
564 // perturbations by Saturn and Uranus
565 double phi, x, y;
566
567 phi = (2*m5-6*m6+3*m7); x = cos(phi); y = sin(phi);
568 dl = dl-0.8*x+8.5*y; dr = dr-0.1*x;
569 addthe(x,y,c3[2],s3[2],x,y);
570 dl = dl+0.4*x+0.5*y; dr = dr-0.7*x+0.5*y; db = db-0.1*x;
571
572 l = p2 * frac(0.0388910 + m5/p2 + ((5025.2+0.8*t)*t+dl)/1296.0e3);
573 b = (227.3 - 0.3*t + db) / arc;
574 r = 5.208873 + 0.000041*t + dr*1.0e-5;
575
576 dl = 14.5+1.41*cos(m5); dr = 3.66*sin(m5); db = 0.33*sin(m5);
577
578 posvel();
579
580 return rp;
581 }
582
583
584 /*--------------------- Saturn ------------------------------------------*/
585
Saturn(double t)586 Vec3 Plan200::Saturn (double t) // position of Saturn at time t
587 {
588 /*
589 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
590 of Saturn for Equinox of Date.
591 t is the time in Julian centuries since J2000.
592 */
593
594 const double arc = 206264.81; // 3600*180/pi = ''/r
595 const double p2 = M_PI * 2.0;
596
597 int i;
598
599 tt = t;
600 dl = 0.0; dr = 0.0; db = 0.0;
601 m5 = p2 * frac(0.0565314+8.4302963*t);
602 m6 = p2 * frac(0.8829867+3.3947688*t);
603 m7 = p2 * frac(0.3969537+1.1902586*t);
604 m8 = p2 * frac(0.7208473+0.6068623*t);
605 c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m6); s3[1] = sin(m6);
606 for (i=2; i<12; i++)
607 addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);
608
609 // perturbations by Jupiter
610 c[6] = 1.0; s[6] = 0.0; c[7] = cos(m5); s[7] = sin(m5);
611 for (i=6; i>0; i--)
612 addthe(c[i],s[i],c[7],-s[7],c[i-1],s[i-1]);
613 term(0, 5, 0, 12.0, -1.4, -13.9, 6.4, 1.2, -1.8);
614 term(0, 4, 0, 0.0, -0.2, -0.9, 1.0, 0.0, -0.1);
615 term(1, 7, 0, 0.9, 0.4, -1.8, 1.9, 0.2, 0.2);
616 term(1, 6, 0, -348.3, 22907.7, -52915.5, -752.2, -3266.5, 8314.4);
617 term(1, 6, 1, -225.2, -146.2, 337.7, -521.3, 79.6, 17.4);
618 term(1, 6, 2, 1.3, -1.4, 3.2, 2.9, 0.1, -0.4);
619 term(1, 5, 0, -1.0, -30.7, 108.6, -815.0, -3.6, -9.3);
620 term(1, 4, 0, -2.0, -2.7, -2.1, -11.9, -0.1, -0.4);
621 term(2, 7, 0, 0.1, 0.2, -1.0, 0.3, 0.0, 0.0);
622 term(2, 6, 0, 44.2, 724.0, -1464.3, -34.7, -188.7, 459.1);
623 term(2, 6, 1, -17.0, -11.3, 18.9, -28.9, 1.0, -3.7);
624 term(2, 5, 0, -3.5, -426.6, -546.5, -26.5, -1.6, -2.7);
625 term(2, 5, 1, 3.5, -2.2, -2.6, -4.3, 0.0, 0.0);
626 term(2, 4, 0, 10.5, -30.9, -130.5, -52.3, -1.9, 0.2);
627 term(2, 3, 0, -0.2, -0.4, -1.2, -0.1, -0.1, 0.0);
628 term(3, 6, 0, 6.5, 30.5, -61.1, 0.4, -11.6, 28.1);
629 term(3, 6, 1, -1.2, -0.7, 1.1, -1.8, -0.2, -0.6);
630 term(3, 5, 0, 29.0, -40.2, 98.2, 45.3, 3.2, -9.4);
631 term(3, 5, 1, 0.6, 0.6, -1.0, 1.3, 0.0, 0.0);
632 term(3, 4, 0, -27.0, -21.1, -68.5, 8.1, -19.8, 5.4);
633 term(3, 4, 1, 0.9, -0.5, -0.4, -2.0, -0.1, -0.8);
634 term(3, 3, 0, -5.4, -4.1, -19.1, 26.2, -0.1, -0.1);
635 term(4, 6, 0, 0.6, 1.4, -3.0, -0.2, -0.6, 1.6);
636 term(4, 5, 0, 1.5, -2.5, 12.4, 4.7, 1.0, -1.1);
637 term(4, 4, 0, -821.9, -9.6, -26.0, 1873.6, -70.5, -4.4);
638 term(4, 4, 1, 4.1, -21.9, -50.3, -9.9, 0.7, -3.0);
639 term(4, 3, 0, -2.0, -4.7, -19.3, 8.2, -0.1, -0.3);
640 term(4, 2, 0, -1.5, 1.3, 6.5, 7.3, 0.0, 0.0);
641 term(5, 4, 0, -2627.6, -1277.3, 117.4, -344.1, -13.8, -4.3);
642 term(5, 4, 1, 63.0, -98.6, 12.7, 6.7, 0.1, -0.2);
643 term(5, 4, 2, 1.7, 1.2, -0.2, 0.3, 0.0, 0.0);
644 term(5, 3, 0, 0.4, -3.6, -11.3, -1.6, 0.0, -0.3);
645 term(5, 2, 0, -1.4, 0.3, 1.5, 6.3, -0.1, 0.0);
646 term(5, 1, 0, 0.3, 0.6, 3.0, -1.7, 0.0, 0.0);
647 term(6, 4, 0, -146.7, -73.7, 166.4, -334.3, -43.6, -46.7);
648 term(6, 4, 1, 5.2, -6.8, 15.1, 11.4, 1.7, -1.0);
649 term(6, 3, 0, 1.5, -2.9, -2.2, -1.3, 0.1, -0.1);
650 term(6, 2, 0, -0.7, -0.2, -0.7, 2.8, 0.0, 0.0);
651 term(6, 1, 0, 0.0, 0.5, 2.5, -0.1, 0.0, 0.0);
652 term(6, 0, 0, 0.3, -0.1, -0.3, -1.2, 0.0, 0.0);
653 term(7, 4, 0, -9.6, -3.9, 9.6, -18.6, -4.7, -5.3);
654 term(7, 4, 1, 0.4, -0.5, 1.0, 0.9, 0.3, -0.1);
655 term(7, 3, 0, 3.0, 5.3, 7.5, -3.5, 0.0, 0.0);
656 term(7, 2, 0, 0.2, 0.4, 1.6, -1.3, 0.0, 0.0);
657 term(7, 1, 0, -0.1, 0.2, 1.0, 0.5, 0.0, 0.0);
658 term(7, 0, 0, 0.2, 0.0, 0.2, -1.0, 0.0, 0.0);
659 term(8, 4, 0, -0.7, -0.2, 0.6, -1.2, -0.4, -0.4);
660 term(8, 3, 0, 0.5, 1.0, -2.0, 1.5, 0.1, 0.2);
661 term(8, 2, 0, 0.4, 1.3, 3.6, -0.9, 0.0, -0.1);
662 term(9, 2, 0, 4.0, -8.7, -19.9, -9.9, 0.2, -0.4);
663 term(9, 2, 1, 0.5, 0.3, 0.8, -1.8, 0.0, 0.0);
664 term(10, 2, 0, 21.3, -16.8, 3.3, 3.3, 0.2, -0.2);
665 term(10, 2, 1, 1.0, 1.7, -0.4, 0.4, 0.0, 0.0);
666 term(11, 2, 0, 1.6, -1.3, 3.0, 3.7, 0.8, -0.2);
667
668 // perturbations by Uranus
669 c[5] = cos(m7); s[5] = -sin(m7);
670 for (i=5; i>1; i--)
671 addthe(c[i],s[i],c[5],s[5],c[i-1],s[i-1]);
672 term(0, 5, 0, 1.0, 0.7, 0.4, -1.5, 0.1, 0.0);
673 term(0, 4, 0, 0.0, -0.4, -1.1, 0.1, -0.1, -0.1);
674 term(0, 3, 0, -0.9, -1.2, -2.7, 2.1, -0.5, -0.3);
675 term(1, 5, 0, 7.8, -1.5, 2.3, 12.7, 0.0, 0.0);
676 term(1, 4, 0, -1.1, -8.1, 5.2, -0.3, -0.3, -0.3);
677 term(1, 3, 0, -16.4, -21.0, -2.1, 0.0, 0.4, 0.0);
678 term(2, 5, 0, 0.6, -0.1, 0.1, 1.2, 0.1, 0.0);
679 term(2, 4, 0, -4.9, -11.7, 31.5, -13.3, 0.0, -0.2);
680 term(2, 3, 0, 19.1, 10.0, -22.1, 42.1, 0.1, -1.1);
681 term(2, 2, 0, 0.9, -0.1, 0.1, 1.4, 0.0, 0.0);
682 term(3, 4, 0, -0.4, -0.9, 1.7, -0.8, 0.0, -0.3);
683 term(3, 3, 0, 2.3, 0.0, 1.0, 5.7, 0.3, 0.3);
684 term(3, 2, 0, 0.3, -0.7, 2.0, 0.7, 0.0, 0.0);
685 term(3, 1, 0, -0.1, -0.4, 1.1, -0.3, 0.0, 0.0);
686
687 // perturbations by Neptune
688 c[5] = cos(m8); s[5] = -sin(m8);
689 addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
690 term(1, 5, 0, -1.3, -1.2, 2.3, -2.5, 0.0, 0.0);
691 term(1, 4, 0, 1.0, -0.1, 0.1, 1.4, 0.0, 0.0);
692 term(2, 4, 0, 1.1, -0.1, 0.2, 3.3, 0.0, 0.0);
693
694 // perturbations by Jupiter and Uranus
695 double phi, x, y;
696
697 phi = (-2*m5+5*m6-3*m7); x = cos(phi); y = sin(phi);
698 dl = dl-0.8*x-0.1*y; dr = dr-0.2*x+1.8*y; db = db+0.3*x+0.5*y;
699 addthe(x,y,c3[1],s3[1],x,y);
700 dl = dl+(2.4-0.7*t)*x+(27.8-0.4*t)*y;
701 dr = dr+2.1*x-0.2*y;
702 addthe(x,y,c3[1],s3[1],x,y);
703 dl = dl+0.1*x+1.6*y; dr = dr-3.6*x+0.3*y; db = db-0.2*x+0.6*y;
704
705 l = p2 * frac(0.2561136+m6/p2+((5018.6+t*1.9)*t+dl)/1296.0e3);
706 b = (175.1 - 10.2*t + db) / arc;
707 r = 9.557584 - 0.000186*t + dr*1.0e-5;
708
709 dl = 5.84+0.65*cos(m6); dr = 3.09*sin(m6); db = 0.24*cos(m6);
710
711 posvel();
712
713 return rp;
714 }
715
716
717 /*--------------------- Uranus ------------------------------------------*/
718
Uranus(double t)719 Vec3 Plan200::Uranus (double t) // position of Uranus at time t
720 {
721 /*
722 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
723 of Uranus for Equinox of Date.
724 t is the time in Julian centuries since J2000.
725 */
726
727 const double arc = 206264.81; // 3600*180/pi = ''/r
728 const double p2 = M_PI * 2.0;
729
730 int i;
731
732 tt = t;
733 dl = 0.0; dr = 0.0; db = 0.0;
734 m5 = p2 * frac(0.0564472+8.4302889*t);
735 m6 = p2 * frac(0.8829611+3.3947583*t);
736 m7 = p2 * frac(0.3967117+1.1902849*t);
737 m8 = p2 * frac(0.7216833+0.6068528*t);
738 c3[2] = 1.0; s3[2] = 0.0; c3[3] = cos(m7); s3[3] = sin(m7);
739 for (i=4; i<10; i++)
740 addthe(c3[i-1],s3[i-1],c3[3],s3[3],c3[i],s3[i]);
741 c3[1] = c3[3];
742 s3[1] = -s3[3];
743 c3[0] = c3[4];
744 s3[0] = -s3[4];
745
746 // perturbations by Jupiter
747 c[8] = 1.0; s[8] = 0.0; c[7] = cos(m5); s[7] = -sin(m5);
748 addthe(c[7],s[7],c[7],s[7],c[6],s[6]);
749 term(1, 7, 0, 0.0, 0.0, -0.1, 1.7, -0.1, 0.0);
750 term(2, 7, 0, 0.5, -1.2, 18.9, 9.1, -0.9, 0.1);
751 term(3, 7, 0, -21.2, 48.7, -455.5, -198.8, 0.0, 0.0);
752 term(3, 6, 0, -0.5, 1.2, -10.9, -4.8, 0.0, 0.0);
753 term(4, 7, 0, -1.3, 3.2, -23.2, -11.1, 0.3, 0.1);
754 term(4, 6, 0, -0.2, 0.2, 1.1, 1.5, 0.0, 0.0);
755 term(5, 7, 0, 0.0, 0.2, -1.8, 0.4, 0.0, 0.0);
756
757 // perturbations by Saturn
758 c[7] = cos(m6); s[7] = -sin(m6);
759 for (i=7; i>4; i--)
760 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
761 term(2, 7, 0, 1.4, -0.5, -6.4, 9.0, -0.4, -0.8);
762 term(3, 7, 0, -18.6, -12.6, 36.7, -336.8, 1.0, 0.3);
763 term(3, 6, 0, -0.7, -0.3, 0.5, -7.5, 0.1, 0.0);
764 term(4, 7, 0, 20.0, -141.6, -587.1, -107.0, 3.1, -0.8);
765 term(4, 7, 1, 1.0, 1.4, 5.8, -4.0, 0.0, 0.0);
766 term(4, 6, 0, 1.6, -3.8, -35.6, -16.0, 0.0, 0.0);
767 term(5, 7, 0, 75.3, -100.9, 128.9, 77.5, -0.8, 0.1);
768 term(5, 7, 1, 0.2, 1.8, -1.9, 0.3, 0.0, 0.0);
769 term(5, 6, 0, 2.3, -1.3, -9.5, -17.9, 0.0, 0.1);
770 term(5, 5, 0, -0.7, -0.5, -4.9, 6.8, 0.0, 0.0);
771 term(6, 7, 0, 3.4, -5.0, 21.6, 14.3, -0.8, -0.5);
772 term(6, 6, 0, 1.9, 0.1, 1.2, -12.1, 0.0, 0.0);
773 term(6, 5, 0, -0.1, -0.4, -3.9, 1.2, 0.0, 0.0);
774 term(6, 4, 0, -0.2, 0.1, 1.6, 1.8, 0.0, 0.0);
775 term(7, 7, 0, 0.2, -0.3, 1.0, 0.6, -0.1, 0.0);
776 term(7, 6, 0, -2.2, -2.2, -7.7, 8.5, 0.0, 0.0);
777 term(7, 5, 0, 0.1, -0.2, -1.4, -0.4, 0.0, 0.0);
778 term(7, 4, 0, -0.1, 0.0, 0.1, 1.2, 0.0, 0.0);
779 term(8, 6, 0, -0.2, -0.6, 1.4, -0.7, 0.0, 0.0);
780
781 // perturbations by Neptune
782 c[7] = cos(m8); s[7] = -sin(m8);
783 for (i=7; i>0; i--)
784 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
785 term(3, 8, 0, -78.1, 19518.1, -90718.2, -334.7, 2759.5, -311.9);
786 term(3, 8, 1, -81.6, 107.7, -497.4, -379.5, -2.8, -43.7);
787 term(3, 8, 2, -6.6, -3.1, 14.4, -30.6, -0.4, -0.5);
788 term(3, 8, 3, 0.0, -0.5, 2.4, 0.0, 0.0, 0.0);
789 term(4, 8, 0, -2.4, 586.1, -2145.2, -15.3, 130.6, -14.3);
790 term(4, 8, 1, -4.5, 6.6, -24.2, -17.8, 0.7, -1.6);
791 term(4, 8, 2, -0.4, 0.0, 0.1, -1.4, 0.0, 0.0);
792 term(5, 8, 0, 0.0, 24.5, -76.2, -0.6, 7.0, -0.7);
793 term(5, 8, 1, -0.2, 0.4, -1.4, -0.8, 0.1, -0.1);
794 term(6, 8, 0, 0.0, 1.1, -3.0, 0.1, 0.4, 0.0);
795 term(1, 7, 0, -0.2, 0.2, 0.7, 0.7, -0.1, 0.0);
796 term(2, 7, 0, -2.8, 2.5, 8.7, 10.5, -0.4, -0.1);
797 term(3, 7, 0, -28.4, 20.3, -51.4, -72.0, 0.0, 0.0);
798 term(3, 6, 0, -0.6, -0.1, 4.2, -14.6, 0.2, 0.4);
799 term(3, 5, 0, 0.2, 0.5, 3.4, -1.6, -0.1, 0.1);
800 term(4, 7, 0, -1.8, 1.3, -5.5, -7.7, 0.0, 0.3);
801 term(4, 6, 0, 29.4, 10.2, -29.0, 83.2, 0.0, 0.0);
802 term(4, 5, 0, 8.8, 17.8, -41.9, 21.5, -0.1, -0.3);
803 term(4, 4, 0, 0.0, 0.1, -2.1, -0.9, 0.1, 0.0);
804 term(5, 6, 0, 1.5, 0.5, -1.7, 5.1, 0.1, -0.2);
805 term(5, 5, 0, 4.4, 14.6, -84.3, 25.2, 0.1, -0.1);
806 term(5, 4, 0, 2.4, -4.5, 12.0, 6.2, 0.0, 0.0);
807 term(5, 3, 0, 2.9, -0.9, 2.1, 6.2, 0.0, 0.0);
808 term(6, 5, 0, 0.3, 1.0, -4.0, 1.1, 0.1, -0.1);
809 term(6, 4, 0, 2.1, -2.7, 17.9, 14.0, 0.0, 0.0);
810 term(6, 3, 0, 3.0, -0.4, 2.3, 17.6, -0.1, -0.1);
811 term(6, 2, 0, -0.6, -0.5, 1.1, -1.6, 0.0, 0.0);
812 term(7, 4, 0, 0.2, -0.2, 1.0, 0.8, 0.0, 0.0);
813 term(7, 3, 0, -0.9, -0.1, 0.6, -7.1, 0.0, 0.0);
814 term(7, 2, 0, -0.5, -0.6, 3.8, -3.6, 0.0, 0.0);
815 term(7, 1, 0, 0.0, -0.5, 3.0, 0.1, 0.0, 0.0);
816 term(8, 2, 0, 0.2, 0.3, -2.7, 1.6, 0.0, 0.0);
817 term(8, 1, 0, -0.1, 0.2, -2.0, -0.4, 0.0, 0.0);
818 term(9, 1, 0, 0.1, -0.2, 1.3, 0.5, 0.0, 0.0);
819 term(9, 0, 0, 0.1, 0.0, 0.4, 0.9, 0.0, 0.0);
820
821 // perturbations by Jupiter and Saturn
822 c[7] = cos(m6); s[7] = -sin(m6);
823 c[4] = cos(-4*m6+2*m5); s[4] = sin(-4*m6+2*m5);
824 for (i=4; i>2; i--)
825 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
826 term(0, 4, 0, -0.7, 0.4, -1.5, -2.5, 0.0, 0.0);
827 term(1, 4, 0, -0.1, -0.1, -2.2, 1.0, 0.0, 0.0);
828 term(3, 3, 0, 0.1, -0.4, 1.4, 0.2, 0.0, 0.0);
829 term(3, 2, 0, 0.4, 0.5, -0.8, -0.8, 0.0, 0.0);
830 term(4, 2, 0, 5.7, 6.3, 28.5, -25.5, 0.0, 0.0);
831 term(4, 2, 1, 0.1, -0.2, -1.1, -0.6, 0.0, 0.0);
832 term(5, 2, 0, -1.4, 29.2, -11.4, 1.1, 0.0, 0.0);
833 term(5, 2, 1, 0.8, -0.4, 0.2, 0.3, 0.0, 0.0);
834 term(6, 2, 0, 0.0, 1.3, -6.0, -0.1, 0.0, 0.0);
835
836 l = p2 * frac(0.4734843+m7/p2+((5082.3+34.2*t)*t+dl)/1296.0E3);
837 b = (-130.61+(-0.54+0.04*t)*t+db)/arc;
838 r = 19.211991+(-0.000333-0.000005*t)*t+dr*1.0e-5;
839
840 dl = 2.05+0.19*cos(m7); dr = 1.86*sin(m7); db = -0.03*sin(m7);
841
842 posvel();
843
844 return rp;
845 }
846
847
848 /*--------------------- Neptune ------------------------------------------*/
849
Neptune(double t)850 Vec3 Plan200::Neptune (double t) // position of Neptune at time t
851 {
852 /*
853 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
854 of Neptune for Equinox of Date.
855 t is the time in Julian centuries since J2000.
856 */
857
858 const double arc = 206264.81; // 3600*180/pi = ''/r
859 const double p2 = M_PI * 2.0;
860
861 int i;
862
863 tt = t;
864 dl = 0.0; dr = 0.0; db = 0.0;
865 m5 = p2 * frac(0.0563867+8.4298907*t);
866 m6 = p2 * frac(0.8825086+3.3957748*t);
867 m7 = p2 * frac(0.3965358+1.1902851*t);
868 m8 = p2 * frac(0.7214906+0.6068526*t);
869 c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m8); s3[1] = sin(m8);
870 for (i=2; i<7; i++)
871 addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);
872
873 // perturbations by Jupiter
874 c[6] = 1.0; s[6] = 0.0; c[5] = cos(m5); s[5] = -sin(m5);
875 addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
876 term(0, 5, 0, 0.1, 0.1, -3.0, 1.8, -0.3, -0.3);
877 term(1, 6, 0, 0.0, 0.0, -15.9, 9.0, 0.0, 0.0);
878 term(1, 5, 0, -17.6, -29.3, 416.1, -250.0, 0.0, 0.0);
879 term(1, 4, 0, -0.4, -0.7, 10.4, -6.2, 0.0, 0.0);
880 term(2, 5, 0, -0.2, -0.4, 2.4, -1.4, 0.4, -0.3);
881
882 // perturbations by Saturn
883 c[6] = 1.0; s[6] = 0.0; c[5] = cos(m6); s[5] = -sin(m6);
884 addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
885 term(0, 5, 0, -0.1, 0.0, 0.2, -1.8, -0.1, -0.5);
886 term(1, 6, 0, 0.0, 0.0, -8.3, -10.4, 0.0, 0.0);
887 term(1, 5, 0, 13.6, -12.7, 187.5, 201.1, 0.0, 0.0);
888 term(1, 4, 0, 0.4, -0.4, 4.5, 4.5, 0.0, 0.0);
889 term(2, 5, 0, 0.4, -0.1, 1.7, -3.2, 0.2, 0.2);
890 term(2, 4, 0, -0.1, 0.0, -0.2, 2.7, 0.0, 0.0);
891
892 // perturbations by Uranus
893 c[6] = 1.0; s[6] = 0.0; c[5] = cos(m7); s[5] = -sin(m7);
894 for (i=5; i>0; i--)
895 addthe(c[i],s[i],c[5],s[5],c[i-1],s[i-1]);
896 term(1, 6, 0, 32.3, 3549.5, -25880.2, 235.8, -6360.5, 374.0);
897 term(1, 6, 1, 31.2, 34.4, -251.4, 227.4, 34.9, 29.3);
898 term(1, 6, 2, -1.4, 3.9, -28.6, -10.1, 0.0, -0.9);
899 term(2, 6, 0, 6.1, 68.0, -111.4, 2.0, -54.7, 3.7);
900 term(2, 6, 1, 0.8, -0.2, -2.1, 2.0, -0.2, 0.8);
901 term(3, 6, 0, 0.1, 1.0, -0.7, 0.0, -0.8, 0.1);
902 term(0, 5, 0, -0.1, -0.3, -3.6, 0.0, 0.0, 0.0);
903 term(1, 6, 0, 0.0, 0.0, 5.5, -6.9, 0.1, 0.0);
904 term(1, 5, 0, -2.2, -1.6, -116.3, 163.6, 0.0, -0.1);
905 term(1, 4, 0, 0.2, 0.1, -1.2, 0.4, 0.0, -0.1);
906 term(2, 5, 0, 4.2, -1.1, -4.4, -34.6, -0.2, 0.1);
907 term(2, 4, 0, 8.6, -2.9, -33.4, -97.0, 0.2, 0.1);
908 term(3, 5, 0, 0.1, -0.2, 2.1, -1.2, 0.0, 0.1);
909 term(3, 4, 0, -4.6, 9.3, 38.2, 19.8, 0.1, 0.1);
910 term(3, 3, 0, -0.5, 1.7, 23.5, 7.0, 0.0, 0.0);
911 term(4, 4, 0, 0.2, 0.8, 3.3, -1.5, -0.2, -0.1);
912 term(4, 3, 0, 0.9, 1.7, 17.9, -9.1, -0.1, 0.0);
913 term(4, 2, 0, -0.4, -0.4, -6.2, 4.8, 0.0, 0.0);
914 term(5, 3, 0, -1.6, -0.5, -2.2, 7.0, 0.0, 0.0);
915 term(5, 2, 0, -0.4, -0.1, -0.7, 5.5, 0.0, 0.0);
916 term(5, 1, 0, 0.2, 0.0, 0.0, -3.5, 0.0, 0.0);
917 term(6, 2, 0, -0.3, 0.2, 2.1, 2.7, 0.0, 0.0);
918 term(6, 1, 0, 0.1, -0.1, -1.4, -1.4, 0.0, 0.0);
919 term(6, 0, 0, -0.1, 0.1, 1.4, 0.7, 0.0, 0.0);
920
921 l = p2 * frac(0.1254046+m8/p2+((4982.8-21.3*t)*t+dl)/1296.0e3);
922 b = (54.77+(0.26+0.06*t)*t+db)/arc;
923 r = 30.072984+(0.001234+0.000003*t)*t+dr*1.0e-5;
924
925 dl = 1.04+0.02*cos(m8); dr = 0.27*sin(m8); db = 0.03*sin(m8);
926
927 posvel();
928
929 return rp;
930 }
931
932 /*--------------------- Pluto ------------------------------------------*/
933
Pluto(double t)934 Vec3 Plan200::Pluto (double t) // position of Pluto at time t
935 {
936 /*
937 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
938 of Pluto for Equinox of Date.
939 t is the time in Julian centuries since J2000.
940 */
941
942 const double arc = 206264.81; // 3600*180/pi = ''/r
943 const double p2 = M_PI * 2.0;
944
945 int i;
946
947 tt = t;
948 dl = 0.0; dr = 0.0; db = 0.0;
949 m5 = p2 * frac(0.0565314+8.4302963*t);
950 m6 = p2 * frac(0.8829867+3.3947688*t);
951 m1 = p2 * frac(0.0385795+0.4026667*t);
952 c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m1); s3[1] = sin(m1);
953 for (i=2; i<7; i++)
954 addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);
955
956 // Kepler terms and perturbations by Jupiter
957 c[3] = 1.0; s[3] = 0.0; c[4] = cos(m5); s[4] = sin(m5);
958 for (i=3; i>1; i--)
959 addthe(c[i],s[i],c[4],-s[4],c[i-1],s[i-1]);
960 addthe(c[4],s[4],c[4],s[4],c[5],s[5]);
961 term(1, 3, 0, 0.06,100924.08, -960396.0,15965.1, 51987.68,-24288.76);
962 term(2, 3, 0, 3274.74,17835.12, -118252.2,3632.4, 12687.49,-6049.72);
963 term(3, 3, 0, 1543.52,4631.99, -21446.6,1167.0, 3504.0,-1853.1);
964 term(4, 3, 0, 688.99,1227.08, -4823.4,213.5, 1048.19,-648.26);
965 term(5, 3, 0, 242.27, 415.93, -1075.4, 140.6, 302.33, -209.76);
966 term(6, 3, 0, 138.41, 110.91, -308.8, -55.3, 109.52, -93.82);
967 term(3, 2, 0, -0.99, 5.06, -25.6, 19.8, 1.26, -1.96);
968 term(2, 2, 0, 7.15, 5.61, -96.7, 57.2, 1.64, -2.16);
969 term(1, 2, 0, 10.79, 23.13, -390.4, 236.4, -0.33, 0.86);
970 term(0, 4, 0, -0.23, 4.43, 102.8, 63.2, 3.15, 0.34);
971 term(1, 4, 0, -1.1, -0.92, 11.8, -2.3, 0.43, 0.14);
972 term(2, 4, 0, 0.62, 0.84, 2.3, 0.7, 0.05, -0.04);
973 term(3, 4, 0, -0.38, -0.45, 1.2, -0.8, 0.04, 0.05);
974 term(4, 4, 0, 0.17, 0.25, 0.0, 0.2, -0.01, -0.01);
975 term(3, 1, 0, 0.06, 0.07, -0.6, 0.3, 0.03, -0.03);
976 term(2, 1, 0, 0.13, 0.2, -2.2, 1.5, 0.03, -0.07);
977 term(1, 1, 0, 0.32, 0.49, -9.4, 5.7, -0.01, 0.03);
978 term(0, 1, 0, -0.04, -0.07, 2.6, -1.5, 0.07, -0.02);
979
980 // perturbations by Saturn
981 c[4] = cos(m6); s[4] = sin(m6);
982 for (i=3; i>1; i--)
983 addthe(c[i],s[i],c[4],-s[4],c[i-1],s[i-1]);
984 term(1, 2, 0, -29.47, 75.97, -106.4, -204.9, -40.71, -17.55);
985 term(0, 4, 0, -13.88, 18.2, 42.6, -46.1, 1.13, 0.43);
986 term(1, 4, 0, 5.81, -23.48, 15.0, -6.8, -7.48, 3.07);
987 term(2, 4, 0, -10.27, 14.16, -7.9, 0.4, 2.43, -0.09);
988 term(3, 4, 0, 6.86, -10.66, 7.3, -0.3, -2.25, 0.69);
989 term(2, 1, 0, 4.23, 2.0, 0.0, -2.2, -0.24, 0.12);
990 term(1, 1, 0, -5.04, -0.83, -9.2, -3.1, 0.79, -0.24);
991 term(0, 1, 0, 4.25, 2.48, -5.9, -3.3, 0.58, 0.02);
992
993 // perturbations by Jupiter and Saturn
994 double x, y;
995 y = (m5-m6); x = cos(y); y = sin(y);
996 dl = dl-9.11*x+0.12*y; dr = dr-3.4*x-3.3*y; db = db+0.81*x+0.78*y;
997 addthe(x,y,c3[1],s3[1],x,y);
998 dl = dl+5.92*x+0.25*y; dr = dr+2.3*x-3.8*y; db = db-0.67*x-0.51*y;
999
1000 l = p2 * frac(0.6232469+m1/p2+dl/1296.0e3);
1001 b = -6.8232495189e-2 + db/arc;
1002 r = 40.7247248+dr*1.0e-5;
1003
1004 dl = 0.69+0.34*cos(m1)+0.12*cos(2*m1)+0.05*cos(3*m1);
1005 dr = 6.66*sin(m1)+1.64*sin(2*m1);
1006 db = -0.08*cos(m1)-0.17*sin(m1)-0.09*sin(2*m1);
1007
1008 // position
1009 double cl, sl, cb, sb;
1010 Mat3 mx;
1011
1012 cl = cos(l); sl = sin(l); cb = cos(b); sb = sin(b);
1013 rp[0] = r*cl*cb; rp[1] = r*sl*cb; rp[2] = r*sb;
1014 mx = pmatecl (-0.5000021, t); // 1950.0 -> Equinox of Date
1015 rp = mxvct (mx, rp);
1016
1017 // velocity
1018 vp[0] = (dr*cl*cb - dl*r*sl*cb - db*r*cl*sb) * 1E-4;
1019 vp[1] = (dr*sl*cb + dl*r*cl*cb - db*r*sl*sb) * 1E-4;
1020 vp[2] = (dr*sb + db*r*cb) * 1E-4;
1021
1022 return rp;
1023 }
1024
1025 /***************************************************************************/
1026 /* =========================================================================
1027 Procedures for calculating the positions of various moons in the
1028 solar system. The algorithms for the procedures of this unit are
1029 taken from the Explanatory Supplement to the Astronomical Almanac
1030 (edited by Kenneth Seidelmann, University Science Books, Mill Valley,
1031 California, 1992. A number of errors in chapter 6 of this book have
1032 been identified and corrected in the code used here.
1033
1034 Copyright :Gerhard HOLTKAMP 25-MAR-2014
1035 ========================================================================= */
1036
1037
1038 /*--------------------- MarPhobos -----------------------------------------*/
MarPhobos(double t,Vec3 & rs,Vec3 & vs)1039 void MarPhobos (double t, Vec3& rs, Vec3& vs)
1040 {
1041 /*
1042 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1043 of the Mars moon Phobos with respect to Mars for Equinox of Date.
1044 t is the time in Julian centuries since J2000.
1045 ====================================================================*/
1046
1047 const double p2 = M_PI * 2.0;
1048 const double ax =6.26974E-5;
1049 const double nm = 1128.844556;
1050 const double ee = 0.0150;
1051 const double gam = 1.10*M_PI/180.0;
1052
1053 double thet, ll, pp, na, ja;
1054 double al, dy, yr;
1055 Mat3 mt1;
1056
1057 // time from 11-NOV-1971
1058 dy = t * 36525.0 + 10278.5;
1059 yr = dy / 365.25;
1060
1061 // normalized angles
1062 thet = p2 * frac ((327.9 - 0.43533*dy) / 360.0);
1063 ll = p2 * frac ((232.41 + nm*dy + 0.00124*yr*yr) / 360.0);
1064 pp = p2 * frac ((278.96 + 0.43526*dy) / 360.0);
1065 na = p2 * frac ((47.39 - 0.0014*yr) / 360.0);
1066 ja = p2 * frac ((37.27 + 0.0008*yr) / 360.0);
1067
1068 // in-orbit-plane position vector
1069 al = ll - pp; // mean anomaly
1070 al = eccanom (al, ee); // eccentric anomaly
1071 rs[0] = ax * (cos(al) - ee);
1072 rs[1] = ax * sqrt(1.0 - ee*ee) * sin(al);
1073 rs[2] = 0.0;
1074
1075 // in-orbit-plane velocity vector
1076 dy = cos(al);
1077 yr = 1.0 - ee*dy;
1078 yr = 1.235083014E-3 / yr;
1079 vs.assign(-yr*sin(al),yr*sqrt(1.0-ee*ee)*dy,0.0);
1080
1081 // convert to equatorial coordinates
1082 al = pp - thet - na;
1083 mt1 = zrot (-al);
1084 rs = mxvct (mt1, rs);
1085 vs = mxvct (mt1, vs);
1086 mt1 = xrot (-gam);
1087 rs = mxvct (mt1, rs);
1088 vs = mxvct (mt1, vs);
1089 mt1 = zrot (-thet);
1090 rs = mxvct (mt1, rs);
1091 vs = mxvct (mt1, vs);
1092 mt1 = xrot (-ja);
1093 rs = mxvct (mt1, rs);
1094 vs = mxvct (mt1, vs);
1095 mt1 = zrot (-na);
1096 rs = mxvct (mt1, rs);
1097 vs = mxvct (mt1, vs);
1098
1099 // convert to ecliptic coordinates of date
1100 mt1 = pmatequ (-0.500002096, t);
1101 rs = mxvct (mt1, rs);
1102 vs = mxvct (mt1, vs);
1103 rs = equecl (t, rs);
1104 vs = equecl (t, vs);
1105
1106 }
1107
1108 /*--------------------- MarDeimos ------------------------------------------*/
MarDeimos(double t,Vec3 & rs,Vec3 & vs)1109 void MarDeimos (double t, Vec3& rs, Vec3& vs)
1110 {
1111 /*
1112 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1113 of the Mars moon Deimos with respect to Mars for Equinox of Date.
1114 t is the time in Julian centuries since J2000.
1115 ====================================================================== */
1116
1117 const double p2 = M_PI * 2.0;
1118 const double ax =1.56828E-4;
1119 const double nm = 285.161888;
1120 const double ee = 0.0004;
1121 const double gam = 1.79*M_PI/180.0;
1122
1123 double thet, ll, pp, na, ja;
1124 double al, dy, yr;
1125 Mat3 mt1;
1126
1127 // time from 11-NOV-1971
1128 dy = t * 36525.0 + 10278.5;
1129 yr = dy / 365.25;
1130
1131 // normalized angles
1132 thet = p2 * frac ((240.38 - 0.01801*dy) / 360.0);
1133 ll = p2 * frac ((196.55 - 0.01801*dy) / 360.0);
1134 ll = p2 * frac ((28.96 + nm*dy - 0.27*sin(ll)) / 360.0);
1135 pp = p2 * frac ((111.7 + 0.01798*dy) / 360.0);
1136 na = p2 * frac ((46.37 - 0.0014*yr) / 360.0);
1137 ja = p2 * frac ((36.62 + 0.0008*yr) / 360.0);
1138
1139 // in-orbit- plane position vector
1140 al = ll - pp; // mean anomaly
1141 al = eccanom (al, ee); // eccentric anomaly
1142 rs[0] = ax * (cos(al) - ee);
1143 rs[1] = ax * sqrt(1.0 - ee*ee) * sin(al);
1144 rs[2] = 0.0;
1145
1146 // in-orbit-plane velocity vector
1147 dy = cos(al);
1148 yr = 1.0 - ee*dy;
1149 yr = 7.8046400669E-4 / yr;
1150 vs.assign(-yr*sin(al),yr*sqrt(1.0-ee*ee)*dy,0.0);
1151
1152 // convert to equatorial coordinates
1153 al = pp - thet - na;
1154 mt1 = zrot (-al);
1155 rs = mxvct (mt1, rs);
1156 vs = mxvct (mt1, vs);
1157 mt1 = xrot (-gam);
1158 rs = mxvct (mt1, rs);
1159 vs = mxvct (mt1, vs);
1160 mt1 = zrot (-thet);
1161 rs = mxvct (mt1, rs);
1162 vs = mxvct (mt1, vs);
1163 mt1 = xrot (-ja);
1164 rs = mxvct (mt1, rs);
1165 vs = mxvct (mt1, vs);
1166 mt1 = zrot (-na);
1167 rs = mxvct (mt1, rs);
1168 vs = mxvct (mt1, vs);
1169
1170 // convert to ecliptic coordinates of date
1171 mt1 = pmatequ (-0.500002096, t);
1172 rs = mxvct (mt1, rs);
1173 vs = mxvct (mt1, vs);
1174 rs = equecl (t, rs);
1175 vs = equecl (t, vs);
1176
1177 }
1178
1179 /*--------------------- PosJIo ------------------------------------------*/
PosJIo(double t)1180 Vec3 PosJIo (double t)
1181 {
1182 /* ---------------------------------------------------------------------- }
1183 Ecliptic coordinates (in A.U.) rs
1184 of Io with respect to Jupiter for Equinox of Date.
1185 t is the time in Julian centuries since J2000.
1186 ======================================================================*/
1187
1188 const double pi2 = M_PI * 2.0;
1189 const double ax = 0.002819347;
1190 const double omj = 99.99754*M_PI/180.0;
1191 const double jj = 1.30691*M_PI/180.0;
1192 const double ii = 3.10401*M_PI/180.0;
1193
1194 double l1, l2, phl, p1, p2, p3, p4;
1195 double pij, th1, th2, psi, gg;
1196 double xi, nu, ze;
1197 double d; // : EXTENDED;
1198 Vec3 rs;
1199 Mat3 mt1;
1200
1201 // general jupiter moon data
1202 d = t*36525.0 + 8544.5;
1203 l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
1204 l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
1205 phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
1206 p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
1207 p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
1208 p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
1209 p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
1210 pij = 0.23510274429; // 13.470395
1211 th1 = pi2* frac((308.365749 - 0.1328061*d)/360.0);
1212 th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
1213 psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
1214 gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);
1215
1216 // io specific lines
1217 xi = -41279.0 * cos(2.0*(l1-l2))*1E-7;
1218 nu = -5596.0*sin(p3-p4) - 2198.0*sin(p1+p2-2.0*(pij+gg))
1219 +1321.0*sin(phl) - 1157.0*sin(l1-2.0*l2+p4)
1220 -1940.0*sin(l1-2.0*l2+p3) + 791.0*sin(l1-2.0*l2+p2)
1221 + 82363.0*sin(2.0*(l1-l2));
1222 nu *= 1E-7;
1223 ze = 7038.0*sin(l1-th1+nu) + 1835.0*sin(l1-th2+nu);
1224 ze *= 1E-7;
1225
1226 // convert to rectangular coordinates in Jovian equatorial
1227
1228 rs.assign (ax * (1+xi) * cos(l1-psi+nu),
1229 ax * (1+xi) * sin(l1-psi+nu),
1230 ax * ze);
1231
1232 // convert into ecliptic of 1950.0
1233 mt1 = xrot (-ii);
1234 rs = mxvct (mt1, rs);
1235 mt1 = zrot (-(psi-omj));
1236 rs = mxvct (mt1, rs);
1237 mt1 = xrot (-jj);
1238 rs = mxvct (mt1, rs);
1239 mt1 = zrot (-omj);
1240 rs = mxvct (mt1, rs);
1241
1242 // convert into ecliptic of date
1243 mt1 = pmatecl (-0.500002096, t);
1244 rs = mxvct (mt1, rs);
1245
1246 return rs;
1247
1248 }
1249
1250 /*--------------------- PosEuropa ------------------------------------------*/
PosEuropa(double t)1251 Vec3 PosEuropa (double t)
1252 {
1253 /*
1254 Ecliptic coordinates (in A.U.) rs
1255 of Europa with respect to Jupiter for Equinox of Date.
1256 t is the time in Julian centuries since J2000.
1257 =====================================================================*/
1258
1259 const double pi2 = M_PI * 2.0;
1260 const double ax = 0.004485872;
1261 const double omj = 99.99754*M_PI/180.0;
1262 const double jj = 1.30691*M_PI/180.0;
1263 const double ii = 3.10401*M_PI/180.0;
1264
1265 double l1, l2, l3, phl, p2, p3, p4;
1266 double pij, th2, th3, psi, gg;
1267 double xi, nu, ze;
1268 double d; // : EXTENDED;
1269 Vec3 rs;
1270 Mat3 mt1;
1271
1272
1273 // general jupiter moon data
1274 d = t*36525.0 + 8544.5;
1275 l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
1276 l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
1277 l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
1278 phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
1279 // p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
1280 p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
1281 p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
1282 p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
1283 pij = 0.23510274429; // 13.470395
1284 th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
1285 th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
1286 psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
1287 gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);
1288
1289 // Europa specific lines
1290 xi = -3187*cos(l2-p3) - 1738*cos(l2-p4)
1291 +93748*cos(l1-l2);
1292 xi *= 1E-7;
1293 nu = -1158*sin(2*(psi-pij)) + 1715*sin(-2*(pij+gg)+th3+psi)
1294 -1846*sin(gg) + 2397*sin(p3-p4) - 3172*sin(phl)
1295 -1993*sin(l2-l3) + 1844*sin(l2-p2)
1296 +6394*sin(l2-p3) + 3451*sin(l2-p4)
1297 +4159*sin(l1-2*l2+p4) + 7571*sin(l1-2*l2+p3)
1298 -1491*sin(l1-2*l2+p3) - 185640*sin(l1-l2)
1299 -803*sin(l1-l3) + 915*sin(2*(l1-l2));
1300 nu *= 1E-7;
1301 ze = 81575*sin(l2-th2+nu) + 4512*sin(l2-th3+nu)
1302 -3286*sin(l2-psi+nu);
1303 ze *= 1E-7;
1304
1305 // convert to rectangular coordinates in Jovian equatorial
1306
1307 rs.assign (ax * (1+xi) * cos(l2-psi+nu),
1308 ax * (1+xi) * sin(l2-psi+nu),
1309 ax * ze);
1310
1311 // convert into ecliptic of 1950.0
1312
1313 mt1 = xrot (-ii);
1314 rs = mxvct (mt1, rs);
1315 mt1 = zrot (-(psi-omj));
1316 rs = mxvct (mt1, rs);
1317 mt1 = xrot (-jj);
1318 rs = mxvct (mt1, rs);
1319 mt1 = zrot (-omj);
1320 rs = mxvct (mt1, rs);
1321
1322 // convert into ecliptic of date
1323 mt1 = pmatecl (-0.500002096, t);
1324 rs = mxvct (mt1, rs);
1325
1326 return rs;
1327
1328 }
1329
1330 /*--------------------- PosGanymede --------------------------------------*/
PosGanymede(double t)1331 Vec3 PosGanymede (double t)
1332 {
1333 /*
1334 Ecliptic coordinates (in A.U.)
1335 of Ganymede with respect to Jupiter for Equinox of Date.
1336 t is the time in Julian centuries since J2000.
1337 =====================================================================*/
1338
1339 const double pi2 = M_PI * 2.0;
1340 const double ax = 0.007155352;
1341 const double omj = 99.99754*M_PI/180.0;
1342 const double jj = 1.30691*M_PI/180.0;
1343 const double ii = 3.10401*M_PI/180.0;
1344
1345 double l1, l2, l3, l4, phl, p1, p2, p3, p4;
1346 double pij, th2, th3, th4, psi, gg;
1347 double xi, nu, ze;
1348 double d; // : EXTENDED;
1349 Vec3 rs;
1350 Mat3 mt1;
1351
1352 // general jupiter moon data
1353 d = t*36525.0 + 8544.5;
1354 l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
1355 l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
1356 l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
1357 l4 = pi2 * frac((84.455823 + 21.57107087517961*d)/360.0);
1358 phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
1359 p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
1360 p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
1361 p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
1362 p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
1363 pij = 0.23510274429; // 13.470395
1364 th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
1365 th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
1366 th4 = pi2* frac((322.746564 - 0.00176018*d)/360.0);
1367 psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
1368 gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);
1369
1370 // Ganymede specific lines
1371 xi = -14691*cos(l3-p3) - 1758*cos(2*(l3-l4))
1372 +6333*cos(l2-l3);
1373 xi *= 1E-7;
1374 nu = -1488*sin(2*(pij+psi)) + 411*sin(-th3+psi)
1375 +346*sin(-th4+psi) - 2338*sin(gg)
1376 +6558*sin(p3-p4) + 523*sin(p1+p3-2*(pij+gg))
1377 +341*sin(phl) - 943*sin(l3-l4)
1378 +29387*sin(l3-p3) + 15800*sin(l3-p4)
1379 +3218*sin(2*(l3-l4)) + 226*sin(3*l3-2*l4)
1380 -12038*sin(l2-l3) - 662*sin(l1-2*l2+p4)
1381 -1246*sin(l1-2*l2+p3) + 699*sin(l1-2*l2+p2)
1382 +217*sin(l1-l3);
1383 nu *= 1E-7;
1384 ze = -2793*sin(l3-th2+nu) + 32387*sin(l3-th3+nu)
1385 +6871*sin(l3-th4+nu) - 16876*sin(l3-psi+nu);
1386 ze *= 1E-7;
1387
1388 // convert to rectangular coordinates in Jovian equatorial
1389
1390 rs.assign (ax * (1+xi) * cos(l3-psi+nu),
1391 ax * (1+xi) * sin(l3-psi+nu),
1392 ax * ze);
1393
1394 // convert into ecliptic of 1950.0
1395
1396 mt1 = xrot (-ii);
1397 rs = mxvct (mt1, rs);
1398 mt1 = zrot (-(psi-omj));
1399 rs = mxvct (mt1, rs);
1400 mt1 = xrot (-jj);
1401 rs = mxvct (mt1, rs);
1402 mt1 = zrot (-omj);
1403 rs = mxvct (mt1, rs);
1404
1405 // convert into ecliptic of date
1406 mt1 = pmatecl (-0.500002096, t);
1407 rs = mxvct (mt1, rs);
1408
1409 return rs;
1410
1411 }
1412
1413 /*--------------------- PosCallisto --------------------------------------*/
PosCallisto(double t)1414 Vec3 PosCallisto (double t)
1415 {
1416 /*
1417 Ecliptic coordinates (in A.U.)
1418 of Callisto with respect to Jupiter for Equinox of Date.
1419 t is the time in Julian centuries since J2000.
1420 =====================================================================*/
1421
1422 const double pi2 = M_PI * 2.0;
1423 const double ax = 0.012585436;
1424 const double omj = 99.99754*M_PI/180.0;
1425 const double jj = 1.30691*M_PI/180.0;
1426 const double ii = 3.10401*M_PI/180.0;
1427
1428 double l3, l4, p3, p4;
1429 double pij, th3, th4, psi, gp, gg, ph2;
1430 double xi, nu, ze;
1431 double d; // : EXTENDED;
1432 Vec3 rs;
1433 Mat3 mt1;
1434
1435 // general jupiter moon data
1436 d = t*36525.0 + 8544.5;
1437 l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
1438 l4 = pi2 * frac((84.455823 + 21.57107087517961*d)/360.0);
1439 // phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
1440 p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
1441 p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
1442 pij = 0.23510274429; // 13.470395
1443 th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
1444 th4 = pi2* frac((322.746564 - 0.00176018*d)/360.0);
1445 psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
1446 gp = pi2 * frac((31.9785280244 + 0.033459733896*d)/360.0);
1447 gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);
1448 ph2 = 0.91009489942; // 52.1445966929
1449
1450 // Callisto specific lines
1451 xi = 1656*cos(l4-p3) - 73328*cos(l4-p4)
1452 +182*cos(l4-pij) - 541*cos(l4+p4-2*(pij+gg))
1453 -269*cos(2*(l4-p4)) + 974*cos(l3-l4);
1454 xi *= 1E-7;
1455 nu = -407*sin(2*(psi-p4)) + 309*sin(-2*p4+th4+psi)
1456 -4840*sin(2*(psi-pij)) + 2074*sin(-th4+psi)
1457 -5605*sin(gg) - 204*sin(2*gg)
1458 -495*sin(5*gp-2*gg+ph2) + 234*sin(p4-pij)
1459 -6112*sin(p3-p4) - 3318*sin(l4-p3)
1460 +145573*sin(l4-p4) + 178*sin(l4-pij-gg)
1461 -363*sin(l4-pij) + 1085*sin(l4+p4-2*(pij+gg))
1462 +672*sin(2*(l4-p4)) + 218*sin(2*(l4-pij-gg))
1463 +167*sin(2*l4-th4-psi) - 142*sin(2*(l4-psi))
1464 +148*sin(l3-2*l4+p4) - 390*sin(l3-l4)
1465 -195*sin(2*(l3-l4)) + 185*sin(3*l3-7*l4+4*p4);
1466 nu *= 1E-7;
1467 ze = 773*sin(l4-2*(pij+gg)+psi+nu)
1468 -5075*sin(l4-th3+nu) + 44300*sin(l4-th4+nu)
1469 -76493*sin(l4-psi+nu);
1470 ze *= 1E-7;
1471
1472 rs.assign(ax * (1+xi) * cos(l4-psi+nu),
1473 ax * (1+xi) * sin(l4-psi+nu),
1474 ax * ze);
1475
1476 // convert into ecliptic of 1950.0
1477
1478 mt1 = xrot (-ii);
1479 rs = mxvct (mt1, rs);
1480 mt1 = zrot (-(psi-omj));
1481 rs = mxvct (mt1, rs);
1482 mt1 = xrot (-jj);
1483 rs = mxvct (mt1, rs);
1484 mt1 = zrot (-omj);
1485 rs = mxvct (mt1, rs);
1486
1487 // convert into ecliptic of date
1488 mt1 = pmatecl (-0.500002096, t);
1489 rs = mxvct (mt1, rs);
1490
1491 return rs;
1492
1493 }
1494
1495 /*-------------------- Mimas, Enceladus, Dione ---------------------*/
PosSMimas(double t)1496 Vec3 PosSMimas (double t)
1497 {
1498 /* Ecliptic coordinates (in A.U.) of Mimas with respect to Saturn
1499 for Equinox of Date.
1500 t is the time in Julian centuries since J2000.
1501 ===================================================================*/
1502 const double ie = 28.0771*M_PI/180.0;
1503 const double oe = 168.8263*M_PI/180.0;
1504
1505 double d, dt, tt, a, n, ecc, gam, th1, l1, p, m;
1506 Vec3 rs;
1507 Mat3 mt1;
1508
1509 d = 36525.0*t + 40452.0;
1510 dt = d / 365.25;
1511 tt = 5.0616*(((36525*t + 18262.577000)/365.25)+1950.0-1866.06);
1512 tt *= M_PI / 180.0;
1513 a = 0.00124171;
1514 n = 381.994516;
1515 ecc = 0.01986;
1516 gam = 0.0274017; // 1.570;
1517 th1 = 49.4 - 365.025*dt;
1518 l1 = 128.839 + n*d - 43.415*sin(tt) - 0.714*sin(3.0*tt) - 0.020*sin(5.0*tt);
1519 p = 107.0 + 365.560*dt;
1520
1521 // in-orbit-plane position vector
1522 m = l1 - p; // mean anomaly
1523 m = m * M_PI / 180.0;
1524 m = eccanom (m, ecc); // eccentric anomaly
1525 rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);
1526
1527 // convert to ecliptic coordinates
1528 m = p - th1;
1529 m = m * M_PI / 180.0;
1530 mt1 = zrot (-m);
1531 rs = mxvct (mt1, rs);
1532 mt1 = xrot (-gam);
1533 rs = mxvct (mt1, rs);
1534 m = th1*M_PI/180.0 - oe;
1535 mt1 = zrot (-m);
1536 rs = mxvct (mt1, rs);
1537 mt1 = xrot (-ie);
1538 rs = mxvct (mt1, rs);
1539 mt1 = zrot (-oe);
1540 rs = mxvct (mt1, rs);
1541
1542 // convert to ecliptic coordinates of date
1543 mt1 = pmatecl (-0.500002096, t);
1544 rs = mxvct (mt1, rs);
1545
1546 return rs;
1547 }
1548
PosSEnceladus(double t)1549 Vec3 PosSEnceladus (double t)
1550 {
1551 /* Ecliptic coordinates (in A.U.) of Enceladus with respect to Saturn
1552 for Equinox of Date.
1553 t is the time in Julian centuries since J2000.
1554 ===================================================================*/
1555 const double ie = 28.0771*M_PI/180.0;
1556 const double oe = 168.8263*M_PI/180.0;
1557
1558 double d, dt, /*tt,*/ a, n, ecc, gam, th1, l1, p, m;
1559 Vec3 rs;
1560 Mat3 mt1;
1561
1562 d = 36525.0*t + 40452.0;
1563 dt = d / 365.25;
1564 a = 0.00158935;
1565 n = 262.7319052;
1566 ecc = 0.00532;
1567 gam = 0.000628319; // 0.036;
1568 th1 = 145.0 - 152.7*dt;
1569 l1 = 200.155 + n*d + 0.256333 * sin((59.4+32.73*dt)*M_PI/180.0)
1570 + 0.217333 * sin((119.2+93.18*dt)*M_PI/180.0);
1571 p = 312.7 + 123.42*dt;
1572
1573 // in-orbit-plane position vector
1574 m = l1 - p; // mean anomaly
1575 m = m * M_PI / 180.0;
1576 m = eccanom (m, ecc); // eccentric anomaly
1577 rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);
1578
1579 // convert to ecliptic coordinates
1580 m = p - th1;
1581 m = m * M_PI / 180.0;
1582 mt1 = zrot (-m);
1583 rs = mxvct (mt1, rs);
1584 mt1 = xrot (-gam);
1585 rs = mxvct (mt1, rs);
1586 m = th1*M_PI/180.0 - oe;
1587 mt1 = zrot (-m);
1588 rs = mxvct (mt1, rs);
1589 mt1 = xrot (-ie);
1590 rs = mxvct (mt1, rs);
1591 mt1 = zrot (-oe);
1592 rs = mxvct (mt1, rs);
1593
1594 // convert to ecliptic coordinates of date
1595 mt1 = pmatecl (-0.500002096, t);
1596 rs = mxvct (mt1, rs);
1597
1598 return rs;
1599 }
1600
PosSDione(double t)1601 Vec3 PosSDione (double t)
1602 {
1603 /* Ecliptic coordinates (in A.U.) of Dione with respect to Saturn
1604 for Equinox of Date.
1605 t is the time in Julian centuries since J2000.
1606 ===================================================================*/
1607
1608 const double ie = 28.0771*M_PI/180.0;
1609 const double oe = 168.8263*M_PI/180.0;
1610
1611 double d, dt, /*tt,*/ a, n, ecc, gam, th1, l1, p, m;
1612 Vec3 rs;
1613 Mat3 mt1;
1614
1615 d = 36525.0*t + 40452.0;
1616 dt = d / 365.25;
1617 a = 0.00252413;
1618 n = 131.534920026;
1619 ecc = 0.001715;
1620 gam = 0.0005044; // 0.0289;
1621 th1 = 228.0 - 30.6197 * dt;
1622 l1 = 255.1183 + n*d - 0.0146667 * sin((59.4+32.73*dt)*M_PI/180.0)
1623 - 0.0125 * sin((119.2+93.18*dt)*M_PI/180.0);
1624 p = 173.6 + 30.8381 * dt;
1625
1626 // in-orbit-plane position vector
1627 m = l1 - p; // mean anomaly
1628 m = m * M_PI / 180.0;
1629 m = eccanom (m, ecc); // eccentric anomaly
1630 rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);
1631
1632 // convert to ecliptic coordinates
1633 m = p - th1;
1634 m = m * M_PI / 180.0;
1635 mt1 = zrot (-m);
1636 rs = mxvct (mt1, rs);
1637 mt1 = xrot (-gam);
1638 rs = mxvct (mt1, rs);
1639 m = th1*M_PI/180.0 - oe;
1640 mt1 = zrot (-m);
1641 rs = mxvct (mt1, rs);
1642 mt1 = xrot (-ie);
1643 rs = mxvct (mt1, rs);
1644 mt1 = zrot (-oe);
1645 rs = mxvct (mt1, rs);
1646
1647 // convert to ecliptic coordinates of date
1648 mt1 = pmatecl (-0.500002096, t);
1649 rs = mxvct (mt1, rs);
1650
1651 return rs;
1652 }
1653
1654 /*--------------------- JupIo --------------------------------------*/
JupIo(double t,Vec3 & rs,Vec3 & vs)1655 void JupIo (double t, Vec3& rs, Vec3& vs)
1656 {
1657 /*
1658 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1659 of Io with respect to Jupiter for Equinox of Date.
1660 t is the time in Julian centuries since J2000.
1661 ====================================================================*/
1662
1663 const double dlt = 3.472222E-4; // 30 sec
1664
1665 double td;
1666
1667 td = dlt / 36525;
1668 rs = PosJIo (t-td);
1669 vs = PosJIo (t+td);
1670 vs -= rs;
1671 td = 0.5/dlt;
1672 vs *= td;
1673 rs = PosJIo (t);
1674
1675 }
1676
1677 /*--------------------- JupEuropa --------------------------------------*/
JupEuropa(double t,Vec3 & rs,Vec3 & vs)1678 void JupEuropa (double t, Vec3& rs, Vec3& vs)
1679 {
1680 /*
1681 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1682 of Europa with respect to Jupiter for Equinox of Date.
1683 t is the time in Julian centuries since J2000.
1684 ====================================================================*/
1685
1686 const double dlt = 3.472222E-4; // 30 sec
1687
1688 double td;
1689
1690 td = dlt / 36525;
1691 rs = PosEuropa (t-td);
1692 vs = PosEuropa (t+td);
1693 vs -= rs;
1694 td = 0.5/dlt;
1695 vs *= td;
1696 rs = PosEuropa (t);
1697
1698 }
1699
1700 /*--------------------- JupGanymede ------------------------------------*/
JupGanymede(double t,Vec3 & rs,Vec3 & vs)1701 void JupGanymede (double t, Vec3& rs, Vec3& vs)
1702 {
1703 /*
1704 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1705 of Ganimede with respect to Jupiter for Equinox of Date.
1706 t is the time in Julian centuries since J2000.
1707 ====================================================================*/
1708
1709 const double dlt = 3.472222E-4; // 30 sec
1710
1711 double td;
1712
1713 td = dlt / 36525;
1714 rs = PosGanymede (t-td);
1715 vs = PosGanymede (t+td);
1716 vs -= rs;
1717 td = 0.5/dlt;
1718 vs *= td;
1719 rs = PosGanymede (t);
1720
1721 }
1722
1723 /*--------------------- JupCallisto ------------------------------------*/
JupCallisto(double t,Vec3 & rs,Vec3 & vs)1724 void JupCallisto (double t, Vec3& rs, Vec3& vs)
1725 {
1726 /*
1727 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1728 of Callisto with respect to Jupiter for Equinox of Date.
1729 t is the time in Julian centuries since J2000.
1730 ====================================================================*/
1731
1732 const double dlt = 3.472222E-4; // 30 sec
1733
1734 double td;
1735
1736 td = dlt / 36525;
1737 rs = PosCallisto (t-td);
1738 vs = PosCallisto (t+td);
1739 vs -= rs;
1740 td = 0.5/dlt;
1741 vs *= td;
1742 rs = PosCallisto (t);
1743
1744 }
1745
1746 /*----------------------- SatRhea ------------------------------------*/
SatRhea(double t,Vec3 & rs,Vec3 & vs)1747 void SatRhea (double t, Vec3& rs, Vec3& vs)
1748 {
1749 /*
1750 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1751 of Rhea with respect to Saturn for Equinox of Date.
1752 t is the time in Julian centuries since J2000.
1753 ====================================================================*/
1754
1755 const double pi2 = M_PI * 2.0;
1756 const double ie = 28.0771*M_PI/180.0;
1757 const double oe = 168.8263*M_PI/180.0;
1758 const double kap = 57.29578*M_PI/180.0;
1759 const double ax = 0.003524;
1760 const double sg0 = 5.7682811893E-3; // sin(0.3305)
1761
1762 double d, td, tt, pp, ot, nt;
1763 double lam, ii, omg, omb, ee;
1764 Mat3 mt1;
1765
1766 // various times
1767 tt = t*36525.0;
1768 d = tt + 40452.0;
1769 td = 0.5219*(tt+40177)/365.25;
1770 tt = d / 365.25;
1771
1772 // Titan perturbations
1773 pp = pi2 * frac((305.0 + 10.2077*tt)/360.0);
1774 ot = pi2 * frac((276.49 + td)/360.0);
1775 nt = pi2 * frac((44.5 - td)/360.0);
1776
1777 // osculating orbit elements
1778 ee = 0.00021*sin(pp) + 0.001*sin(ot);
1779 omg = 0.00021*cos(pp) + 0.001*cos(ot);
1780 omb = atan2(ee,omg);
1781 ii = sin(omb);
1782 if (fabs(ii) > 0.5) ee = fabs(ee/ii);
1783 else ee = fabs(omg/cos(omb));
1784 pp = pi2 * frac((356.87 - 10.2077*tt)/360.0);
1785 ot = sin(pp);
1786 lam = pi2 * frac((359.4727 + 79.69004007*d)/360.0);
1787 lam = lam + kap*sg0*tan(0.5*ie)*ot;
1788 ii = ie - 7.9412480966E-4 + kap*sg0*cos(pp)
1789 + 3.5081117965E-4*cos(nt);
1790 omg = oe - 1.3613568166E-4
1791 + (kap*sg0*ot + 3.5081117965E-4*sin(nt))/sin(ie);
1792
1793 // in-orbit-plane position vector
1794 pp = lam - omb; // mean anomaly
1795 pp = eccanom (pp, ee); // eccentric anomaly
1796 rs.assign (ax * (cos(pp) - ee), ax * sqrt(1.0 - ee*ee) * sin(pp), 0);
1797
1798 // in-orbit-plane velocity vector
1799 d = cos(pp);
1800 td = 1.0 - ee*d;
1801 td = 4.9000413575E-3 / td;
1802 vs.assign (-td*sin(pp), td*sqrt(1.0-ee*ee)*d, 0);
1803
1804 // convert to ecliptic coordinates
1805 pp = omb - omg;
1806 mt1 = zrot (-pp);
1807 rs = mxvct (mt1, rs);
1808 vs = mxvct (mt1, vs);
1809 mt1 = xrot (-ii);
1810 rs = mxvct (mt1, rs);
1811 vs = mxvct (mt1, vs);
1812 mt1 = zrot (-omg);
1813 rs = mxvct (mt1, rs);
1814 vs = mxvct (mt1, vs);
1815
1816 // convert to ecliptic coordinates of date
1817 mt1 = pmatecl (-0.500002096, t);
1818 rs = mxvct (mt1, rs);
1819 vs = mxvct (mt1, vs);
1820
1821 }
1822
1823 /*----------------------- SatTitan ------------------------------------*/
SatTitan(double t,Vec3 & rs,Vec3 & vs)1824 void SatTitan (double t, Vec3& rs, Vec3& vs)
1825 {
1826 /*
1827 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1828 of Titan with respect to Saturn for Equinox of Date.
1829 t is the time in Julian centuries since J2000.
1830 ===================================================================*/
1831
1832 const double pi2 = M_PI * 2.0;
1833 const double ie = 28.0771*M_PI/180.0;
1834 const double oe = 168.8263*M_PI/180.0;
1835 const double kap = 57.29578*M_PI/180.0;
1836 const double ax = 0.00816765;
1837 const double sg0 = 5.7682811893E-3; // sin(0.3305)
1838
1839 double d, tt, ts, gg, ls, is, os, lms, psi;
1840 double lam, ii, omg, omb, ee, s4;
1841 Mat3 mt1;
1842
1843 // various times
1844 ts = t + 1.0;
1845 tt = t*36525.0;
1846 d = tt+40177;
1847 tt = d / 365.25;
1848
1849 // solar perturbations
1850 ls = pi2 * frac((175.4762 + 1221.5515*ts)/360.0);
1851 is = pi2 * frac((2.4891 + 0.002435*ts)/360.0);
1852 os = pi2 * frac((113.35 - 0.2597*ts)/360.0);
1853 lms = pi2* frac((267.2635 + 1222.1136*ts)/360.0);
1854
1855 gg = pi2 * frac((41.28 - 0.5219*tt)/360.0);
1856 ii = ie - 1.0828022679E-2 + kap*5.2185107774E-3*cos(gg);
1857 s4 = sin(gg);
1858 omg = oe - 2.4748768793E-3 + kap*5.2185107774E-3*s4/sin(ie);
1859 omb = pi2 * frac((275.837 + 0.5219*tt)/360.0);
1860
1861 psi = sin(is)*sin(omg-os);
1862 gg = cos(omg-os);
1863 ts = cos(is)*sin(ii) - sin(is)*cos(ii)*gg;
1864 psi = atan2(psi,ts);
1865 ts = cos(is)*sin(ii)*gg - sin(is)*cos(ii);
1866 gg = sin(ii)*sin(omg-os);
1867 gg = atan2(gg,ts);
1868 lms = lms - gg - os;
1869 gg = omb - omg - psi;
1870
1871 // osculating orbit elements
1872 ts = 2*(lms-gg);
1873 ee = 0.028815 - 0.000184*cos(2*gg) + 0.000073*cos(ts);
1874 omb = omb + kap*(0.00630*sin(2*gg) + 0.0025*sin(ts));
1875 gg = 2*lms + psi;
1876 lam = pi2 * frac((261.3121 + 22.57697385*d)/360.0)
1877 +kap*(sg0*tan(0.5*ie)*s4 - 0.000176*sin(ls)
1878 -0.000215*sin(2*lms) + 0.000057*sin(gg));
1879 ii = ii + 0.000232*kap*cos(gg);
1880 omg = omg + 0.000503*kap*sin(gg);
1881
1882 // in-orbit-plane position vector
1883 gg = lam - omb; // mean anomaly
1884 if (gg < 0) gg += pi2;
1885 gg = eccanom (gg, ee); // eccentric anomaly
1886 rs.assign (ax * (cos(gg) - ee), ax * sqrt(1.0 - ee*ee) * sin(gg), 0);
1887
1888 // in-orbit-plane velocity vector
1889 d = cos(gg);
1890 ts = 1.0 - ee*d;
1891 ts = 3.2183196288E-3 / ts;
1892 vs.assign (-ts*sin(gg), ts*sqrt(1.0-ee*ee)*d, 0);
1893
1894 // convert to ecliptic coordinates
1895 gg = omb - omg;
1896 mt1 = zrot (-gg);
1897 rs = mxvct (mt1, rs);
1898 vs = mxvct (mt1, vs);
1899 mt1 = xrot (-ii);
1900 rs = mxvct (mt1, rs);
1901 vs = mxvct (mt1, vs);
1902 mt1 = zrot (-omg);
1903 rs = mxvct (mt1, rs);
1904 vs = mxvct (mt1, vs);
1905
1906 // convert to ecliptic coordinates of date
1907 mt1 = pmatecl (-0.500002096, t);
1908 rs = mxvct (mt1, rs);
1909 vs = mxvct (mt1, vs);
1910
1911 }
1912
1913
1914 /*----------------------- NepTriton ------------------------------------*/
NepTriton(double t,Vec3 & rs,Vec3 & vs)1915 void NepTriton (double t, Vec3& rs, Vec3& vs)
1916 {
1917 /*
1918 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1919 of Triton with respect to Neptune for Equinox of Date.
1920 t is the time in Julian centuries since J2000.
1921 ====================================================================*/
1922
1923 const double pi2 = M_PI * 2.0;
1924 const double ax = 0.00237142;
1925 const double gam = 158.996*M_PI/180.0;
1926
1927 double d, nn, u, jj, ap;
1928 Mat3 mt1;
1929
1930 // days since epoch
1931 d = t*36525 + 18262.5;
1932
1933 // pole of fixed plane in 1950.0
1934 nn = pi2 * frac((359.28 + 54.308*t)/360.0);
1935 rs.assign (1.0,
1936 pi2 * frac((298.72+2.58*sin(nn)-0.04*sin(2*nn))/360.0),
1937 pi2 * frac((42.63-1.9*cos(nn)+0.01*cos(2*nn))/360.0));
1938 rs = polcar (rs);
1939 mt1 = pmatequ (0.0, -0.500002096);
1940 rs = mxvct (mt1, rs);
1941 rs = carpol (rs);
1942 jj = M_PI/2.0;
1943 ap = rs[1] + jj; // R.A. of ascending node
1944 jj = jj - rs[2]; // inclination
1945
1946 // in-orbit-plane position vector
1947
1948 u = pi2 * frac((200.913 + 61.2588532*d)/360.0);
1949 rs.assign (ax * cos(u), ax * sin(u), 0);
1950
1951 // velocity is perpendicular to radius vector for e=0
1952 u = 1.0 / abs(rs);
1953 vs = u * rs;
1954 mt1 = zrot (-M_PI/2.0);
1955 vs = mxvct (mt1, vs);
1956 u = 2.53538612E-3; // mean orbital speed
1957 vs *= u;
1958
1959 // convert to equatorial coordinates
1960 mt1 = xrot (-gam);
1961 rs = mxvct (mt1, rs);
1962 vs = mxvct (mt1, vs);
1963 nn = pi2 * frac((151.401 + 0.57806*d/365.25)/360.0);
1964 mt1 = zrot (-nn);
1965 rs = mxvct (mt1, rs);
1966 vs = mxvct (mt1, vs);
1967 mt1 = xrot (-jj);
1968 rs = mxvct (mt1, rs);
1969 vs = mxvct (mt1, vs);
1970 mt1 = zrot (-ap);
1971 rs = mxvct (mt1, rs);
1972 vs = mxvct (mt1, vs);
1973
1974 // convert to ecliptic coordinates of date
1975 mt1 = pmatequ (-0.500002096, t);
1976 rs = mxvct (mt1, rs);
1977 vs = mxvct (mt1, vs);
1978 rs = equecl (t, rs);
1979 vs = equecl (t, vs);
1980
1981 }
1982
1983
1984 /*----------------------- PluCharon ------------------------------------*/
PluCharon(double t,Vec3 & rs,Vec3 & vs)1985 void PluCharon (double t, Vec3& rs, Vec3& vs)
1986 {
1987 /*
1988 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
1989 of Charon with respect to Pluto for Equinox of Date.
1990 t is the time in Julian centuries since J2000.
1991 ====================================================================*/
1992
1993 const double pi2 = M_PI * 2.0;
1994 const double ax = 0.00013102;
1995 const double jj = 94.3*M_PI/180.0;
1996 const double nn = 223.7*M_PI/180.0;
1997
1998 double d, u;
1999 Mat3 mt1;
2000
2001 // days since epoch
2002 d = t*36525 + 6544.5;
2003
2004 // in-orbit-plane position vector
2005 u = pi2 * frac((78.6 + 56.3625*d)/360.0);
2006 rs.assign (ax * cos(u), ax * sin(u), 0);
2007
2008 // velocity is perpendicular to radius vector for e=0
2009 d = 1.0 / abs(rs);
2010 vs = d * rs;
2011 mt1 = zrot (-M_PI/2.0);
2012 vs = mxvct (mt1, vs);
2013 d = 1.288837578E-4; // mean orbital speed
2014 vs *= d;
2015
2016 // convert to equatorial coordinates
2017 mt1 = xrot (-jj);
2018 rs = mxvct (mt1, rs);
2019 vs = mxvct (mt1, vs);
2020 mt1 = zrot (-nn);
2021 rs = mxvct (mt1, rs);
2022 vs = mxvct (mt1, vs);
2023
2024 // convert to ecliptic coordinates of date
2025 mt1 = pmatequ (-0.500002096, t);
2026 rs = mxvct (mt1, rs);
2027 vs = mxvct (mt1, vs);
2028 rs = equecl (t, rs);
2029 vs = equecl (t, vs);
2030
2031 }
2032
2033