line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
/* |
2
|
|
|
|
|
|
|
*+ |
3
|
|
|
|
|
|
|
* Name: |
4
|
|
|
|
|
|
|
* palPertue |
5
|
|
|
|
|
|
|
|
6
|
|
|
|
|
|
|
* Purpose: |
7
|
|
|
|
|
|
|
* Update the universal elements by applying planetary perturbations |
8
|
|
|
|
|
|
|
|
9
|
|
|
|
|
|
|
* Language: |
10
|
|
|
|
|
|
|
* Starlink ANSI C |
11
|
|
|
|
|
|
|
|
12
|
|
|
|
|
|
|
* Type of Module: |
13
|
|
|
|
|
|
|
* Library routine |
14
|
|
|
|
|
|
|
|
15
|
|
|
|
|
|
|
* Invocation: |
16
|
|
|
|
|
|
|
* void palPertue( double date, double u[13], int *jstat ); |
17
|
|
|
|
|
|
|
|
18
|
|
|
|
|
|
|
* Arguments: |
19
|
|
|
|
|
|
|
* date = double (Given) |
20
|
|
|
|
|
|
|
* Final epoch (TT MJD) for the update elements. |
21
|
|
|
|
|
|
|
* u = const double [13] (Given & Returned) |
22
|
|
|
|
|
|
|
* Universal orbital elements (Note 1) |
23
|
|
|
|
|
|
|
* (0) combined mass (M+m) |
24
|
|
|
|
|
|
|
* (1) total energy of the orbit (alpha) |
25
|
|
|
|
|
|
|
* (2) reference (osculating) epoch (t0) |
26
|
|
|
|
|
|
|
* (3-5) position at reference epoch (r0) |
27
|
|
|
|
|
|
|
* (6-8) velocity at reference epoch (v0) |
28
|
|
|
|
|
|
|
* (9) heliocentric distance at reference epoch |
29
|
|
|
|
|
|
|
* (10) r0.v0 |
30
|
|
|
|
|
|
|
* (11) date (t) |
31
|
|
|
|
|
|
|
* (12) universal eccentric anomaly (psi) of date, approx |
32
|
|
|
|
|
|
|
* jstat = int * (Returned) |
33
|
|
|
|
|
|
|
* status: |
34
|
|
|
|
|
|
|
* +102 = warning, distant epoch |
35
|
|
|
|
|
|
|
* +101 = warning, large timespan ( > 100 years) |
36
|
|
|
|
|
|
|
* +1 to +10 = coincident with major planet (Note 5) |
37
|
|
|
|
|
|
|
* 0 = OK |
38
|
|
|
|
|
|
|
* -1 = numerical error |
39
|
|
|
|
|
|
|
|
40
|
|
|
|
|
|
|
* Description: |
41
|
|
|
|
|
|
|
* Update the universal elements of an asteroid or comet by applying |
42
|
|
|
|
|
|
|
* planetary perturbations. |
43
|
|
|
|
|
|
|
|
44
|
|
|
|
|
|
|
* Authors: |
45
|
|
|
|
|
|
|
* PTW: Pat Wallace (STFC) |
46
|
|
|
|
|
|
|
* TIMJ: Tim Jenness (JAC, Hawaii) |
47
|
|
|
|
|
|
|
* {enter_new_authors_here} |
48
|
|
|
|
|
|
|
|
49
|
|
|
|
|
|
|
* Notes: |
50
|
|
|
|
|
|
|
* - The "universal" elements are those which define the orbit for the |
51
|
|
|
|
|
|
|
* purposes of the method of universal variables (see reference 2). |
52
|
|
|
|
|
|
|
* They consist of the combined mass of the two bodies, an epoch, |
53
|
|
|
|
|
|
|
* and the position and velocity vectors (arbitrary reference frame) |
54
|
|
|
|
|
|
|
* at that epoch. The parameter set used here includes also various |
55
|
|
|
|
|
|
|
* quantities that can, in fact, be derived from the other |
56
|
|
|
|
|
|
|
* information. This approach is taken to avoiding unnecessary |
57
|
|
|
|
|
|
|
* computation and loss of accuracy. The supplementary quantities |
58
|
|
|
|
|
|
|
* are (i) alpha, which is proportional to the total energy of the |
59
|
|
|
|
|
|
|
* orbit, (ii) the heliocentric distance at epoch, (iii) the |
60
|
|
|
|
|
|
|
* outwards component of the velocity at the given epoch, (iv) an |
61
|
|
|
|
|
|
|
* estimate of psi, the "universal eccentric anomaly" at a given |
62
|
|
|
|
|
|
|
* date and (v) that date. |
63
|
|
|
|
|
|
|
* - The universal elements are with respect to the J2000 equator and |
64
|
|
|
|
|
|
|
* equinox. |
65
|
|
|
|
|
|
|
* - The epochs DATE, U(3) and U(12) are all Modified Julian Dates |
66
|
|
|
|
|
|
|
* (JD-2400000.5). |
67
|
|
|
|
|
|
|
* - The algorithm is a simplified form of Encke's method. It takes as |
68
|
|
|
|
|
|
|
* a basis the unperturbed motion of the body, and numerically |
69
|
|
|
|
|
|
|
* integrates the perturbing accelerations from the major planets. |
70
|
|
|
|
|
|
|
* The expression used is essentially Sterne's 6.7-2 (reference 1). |
71
|
|
|
|
|
|
|
* Everhart and Pitkin (reference 2) suggest rectifying the orbit at |
72
|
|
|
|
|
|
|
* each integration step by propagating the new perturbed position |
73
|
|
|
|
|
|
|
* and velocity as the new universal variables. In the present |
74
|
|
|
|
|
|
|
* routine the orbit is rectified less frequently than this, in order |
75
|
|
|
|
|
|
|
* to gain a slight speed advantage. However, the rectification is |
76
|
|
|
|
|
|
|
* done directly in terms of position and velocity, as suggested by |
77
|
|
|
|
|
|
|
* Everhart and Pitkin, bypassing the use of conventional orbital |
78
|
|
|
|
|
|
|
* elements. |
79
|
|
|
|
|
|
|
* |
80
|
|
|
|
|
|
|
* The f(q) part of the full Encke method is not used. The purpose |
81
|
|
|
|
|
|
|
* of this part is to avoid subtracting two nearly equal quantities |
82
|
|
|
|
|
|
|
* when calculating the "indirect member", which takes account of the |
83
|
|
|
|
|
|
|
* small change in the Sun's attraction due to the slightly displaced |
84
|
|
|
|
|
|
|
* position of the perturbed body. A simpler, direct calculation in |
85
|
|
|
|
|
|
|
* double precision proves to be faster and not significantly less |
86
|
|
|
|
|
|
|
* accurate. |
87
|
|
|
|
|
|
|
* |
88
|
|
|
|
|
|
|
* Apart from employing a variable timestep, and occasionally |
89
|
|
|
|
|
|
|
* "rectifying the orbit" to keep the indirect member small, the |
90
|
|
|
|
|
|
|
* integration is done in a fairly straightforward way. The |
91
|
|
|
|
|
|
|
* acceleration estimated for the middle of the timestep is assumed |
92
|
|
|
|
|
|
|
* to apply throughout that timestep; it is also used in the |
93
|
|
|
|
|
|
|
* extrapolation of the perturbations to the middle of the next |
94
|
|
|
|
|
|
|
* timestep, to predict the new disturbed position. There is no |
95
|
|
|
|
|
|
|
* iteration within a timestep. |
96
|
|
|
|
|
|
|
* |
97
|
|
|
|
|
|
|
* Measures are taken to reach a compromise between execution time |
98
|
|
|
|
|
|
|
* and accuracy. The starting-point is the goal of achieving |
99
|
|
|
|
|
|
|
* arcsecond accuracy for ordinary minor planets over a ten-year |
100
|
|
|
|
|
|
|
* timespan. This goal dictates how large the timesteps can be, |
101
|
|
|
|
|
|
|
* which in turn dictates how frequently the unperturbed motion has |
102
|
|
|
|
|
|
|
* to be recalculated from the osculating elements. |
103
|
|
|
|
|
|
|
* |
104
|
|
|
|
|
|
|
* Within predetermined limits, the timestep for the numerical |
105
|
|
|
|
|
|
|
* integration is varied in length in inverse proportion to the |
106
|
|
|
|
|
|
|
* magnitude of the net acceleration on the body from the major |
107
|
|
|
|
|
|
|
* planets. |
108
|
|
|
|
|
|
|
* |
109
|
|
|
|
|
|
|
* The numerical integration requires estimates of the major-planet |
110
|
|
|
|
|
|
|
* motions. Approximate positions for the major planets (Pluto |
111
|
|
|
|
|
|
|
* alone is omitted) are obtained from the routine palPlanet. Two |
112
|
|
|
|
|
|
|
* levels of interpolation are used, to enhance speed without |
113
|
|
|
|
|
|
|
* significantly degrading accuracy. At a low frequency, the routine |
114
|
|
|
|
|
|
|
* palPlanet is called to generate updated position+velocity "state |
115
|
|
|
|
|
|
|
* vectors". The only task remaining to be carried out at the full |
116
|
|
|
|
|
|
|
* frequency (i.e. at each integration step) is to use the state |
117
|
|
|
|
|
|
|
* vectors to extrapolate the planetary positions. In place of a |
118
|
|
|
|
|
|
|
* strictly linear extrapolation, some allowance is made for the |
119
|
|
|
|
|
|
|
* curvature of the orbit by scaling back the radius vector as the |
120
|
|
|
|
|
|
|
* linear extrapolation goes off at a tangent. |
121
|
|
|
|
|
|
|
* |
122
|
|
|
|
|
|
|
* Various other approximations are made. For example, perturbations |
123
|
|
|
|
|
|
|
* by Pluto and the minor planets are neglected and relativistic |
124
|
|
|
|
|
|
|
* effects are not taken into account. |
125
|
|
|
|
|
|
|
* |
126
|
|
|
|
|
|
|
* In the interests of simplicity, the background calculations for |
127
|
|
|
|
|
|
|
* the major planets are carried out en masse. The mean elements and |
128
|
|
|
|
|
|
|
* state vectors for all the planets are refreshed at the same time, |
129
|
|
|
|
|
|
|
* without regard for orbit curvature, mass or proximity. |
130
|
|
|
|
|
|
|
* |
131
|
|
|
|
|
|
|
* The Earth-Moon system is treated as a single body when the body is |
132
|
|
|
|
|
|
|
* distant but as separate bodies when closer to the EMB than the |
133
|
|
|
|
|
|
|
* parameter RNE, which incurs a time penalty but improves accuracy |
134
|
|
|
|
|
|
|
* for near-Earth objects. |
135
|
|
|
|
|
|
|
* |
136
|
|
|
|
|
|
|
* - This routine is not intended to be used for major planets. |
137
|
|
|
|
|
|
|
* However, if major-planet elements are supplied, sensible results |
138
|
|
|
|
|
|
|
* will, in fact, be produced. This happens because the routine |
139
|
|
|
|
|
|
|
* checks the separation between the body and each of the planets and |
140
|
|
|
|
|
|
|
* interprets a suspiciously small value (0.001 AU) as an attempt to |
141
|
|
|
|
|
|
|
* apply the routine to the planet concerned. If this condition is |
142
|
|
|
|
|
|
|
* detected, the contribution from that planet is ignored, and the |
143
|
|
|
|
|
|
|
* status is set to the planet number (1-10 = Mercury, Venus, EMB, |
144
|
|
|
|
|
|
|
* Mars, Jupiter, Saturn, Uranus, Neptune, Earth, Moon) as a warning. |
145
|
|
|
|
|
|
|
|
146
|
|
|
|
|
|
|
* See Also: |
147
|
|
|
|
|
|
|
* - Sterne, Theodore E., "An Introduction to Celestial Mechanics", |
148
|
|
|
|
|
|
|
* Interscience Publishers Inc., 1960. Section 6.7, p199. |
149
|
|
|
|
|
|
|
* - Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983. |
150
|
|
|
|
|
|
|
|
151
|
|
|
|
|
|
|
* History: |
152
|
|
|
|
|
|
|
* 2012-03-12 (TIMJ): |
153
|
|
|
|
|
|
|
* Initial version direct conversion of SLA/F. |
154
|
|
|
|
|
|
|
* Adapted with permission from the Fortran SLALIB library. |
155
|
|
|
|
|
|
|
* 2012-06-21 (TIMJ): |
156
|
|
|
|
|
|
|
* Support a lack of copysign() function. |
157
|
|
|
|
|
|
|
* 2012-06-22 (TIMJ): |
158
|
|
|
|
|
|
|
* Check __STDC_VERSION__ |
159
|
|
|
|
|
|
|
* {enter_further_changes_here} |
160
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
* Copyright: |
162
|
|
|
|
|
|
|
* Copyright (C) 2004 Patrick T. Wallace |
163
|
|
|
|
|
|
|
* Copyright (C) 2012 Science and Technology Facilities Council. |
164
|
|
|
|
|
|
|
* All Rights Reserved. |
165
|
|
|
|
|
|
|
|
166
|
|
|
|
|
|
|
* Licence: |
167
|
|
|
|
|
|
|
* This program is free software; you can redistribute it and/or |
168
|
|
|
|
|
|
|
* modify it under the terms of the GNU General Public License as |
169
|
|
|
|
|
|
|
* published by the Free Software Foundation; either version 3 of |
170
|
|
|
|
|
|
|
* the License, or (at your option) any later version. |
171
|
|
|
|
|
|
|
* |
172
|
|
|
|
|
|
|
* This program is distributed in the hope that it will be |
173
|
|
|
|
|
|
|
* useful, but WITHOUT ANY WARRANTY; without even the implied |
174
|
|
|
|
|
|
|
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR |
175
|
|
|
|
|
|
|
* PURPOSE. See the GNU General Public License for more details. |
176
|
|
|
|
|
|
|
* |
177
|
|
|
|
|
|
|
* You should have received a copy of the GNU General Public License |
178
|
|
|
|
|
|
|
* along with this program; if not, write to the Free Software |
179
|
|
|
|
|
|
|
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, |
180
|
|
|
|
|
|
|
* MA 02110-1301, USA. |
181
|
|
|
|
|
|
|
|
182
|
|
|
|
|
|
|
* Bugs: |
183
|
|
|
|
|
|
|
* {note_any_bugs_here} |
184
|
|
|
|
|
|
|
*- |
185
|
|
|
|
|
|
|
*/ |
186
|
|
|
|
|
|
|
|
187
|
|
|
|
|
|
|
/* Use the config file if we have one, else look at |
188
|
|
|
|
|
|
|
compiler defines to see if we have C99 */ |
189
|
|
|
|
|
|
|
#if HAVE_CONFIG_H |
190
|
|
|
|
|
|
|
#include |
191
|
|
|
|
|
|
|
#else |
192
|
|
|
|
|
|
|
#ifdef __STDC_VERSION__ |
193
|
|
|
|
|
|
|
# if (__STDC_VERSION__ >= 199901L) |
194
|
|
|
|
|
|
|
# define HAVE_COPYSIGN 1 |
195
|
|
|
|
|
|
|
# endif |
196
|
|
|
|
|
|
|
#endif |
197
|
|
|
|
|
|
|
#endif |
198
|
|
|
|
|
|
|
|
199
|
|
|
|
|
|
|
#include |
200
|
|
|
|
|
|
|
|
201
|
|
|
|
|
|
|
#include "pal.h" |
202
|
|
|
|
|
|
|
#include "palmac.h" |
203
|
|
|
|
|
|
|
#include "pal1sofa.h" |
204
|
|
|
|
|
|
|
|
205
|
|
|
|
|
|
|
/* copysign is C99 */ |
206
|
|
|
|
|
|
|
#if HAVE_COPYSIGN |
207
|
|
|
|
|
|
|
# define COPYSIGN copysign |
208
|
|
|
|
|
|
|
#else |
209
|
|
|
|
|
|
|
# define COPYSIGN(a,b) DSIGN(a,b) |
210
|
|
|
|
|
|
|
#endif |
211
|
|
|
|
|
|
|
|
212
|
2
|
|
|
|
|
|
void palPertue( double date, double u[13], int *jstat ) { |
213
|
|
|
|
|
|
|
|
214
|
|
|
|
|
|
|
/* Distance from EMB at which Earth and Moon are treated separately */ |
215
|
|
|
|
|
|
|
const double RNE=1.0; |
216
|
|
|
|
|
|
|
|
217
|
|
|
|
|
|
|
/* Coincidence with major planet distance */ |
218
|
|
|
|
|
|
|
const double COINC=0.0001; |
219
|
|
|
|
|
|
|
|
220
|
|
|
|
|
|
|
/* Coefficient relating timestep to perturbing force */ |
221
|
|
|
|
|
|
|
const double TSC=1e-4; |
222
|
|
|
|
|
|
|
|
223
|
|
|
|
|
|
|
/* Minimum and maximum timestep (days) */ |
224
|
|
|
|
|
|
|
const double TSMIN = 0.01; |
225
|
|
|
|
|
|
|
const double TSMAX = 10.0; |
226
|
|
|
|
|
|
|
|
227
|
|
|
|
|
|
|
/* Age limit for major-planet state vector (days) */ |
228
|
|
|
|
|
|
|
const double AGEPMO=5.0; |
229
|
|
|
|
|
|
|
|
230
|
|
|
|
|
|
|
/* Age limit for major-planet mean elements (days) */ |
231
|
|
|
|
|
|
|
const double AGEPEL=50.0; |
232
|
|
|
|
|
|
|
|
233
|
|
|
|
|
|
|
/* Margin for error when deciding whether to renew the planetary data */ |
234
|
|
|
|
|
|
|
const double TINY=1e-6; |
235
|
|
|
|
|
|
|
|
236
|
|
|
|
|
|
|
/* Age limit for the body's osculating elements (before rectification) */ |
237
|
|
|
|
|
|
|
const double AGEBEL=100.0; |
238
|
|
|
|
|
|
|
|
239
|
|
|
|
|
|
|
/* Gaussian gravitational constant squared */ |
240
|
|
|
|
|
|
|
const double GCON2 = PAL__GCON * PAL__GCON; |
241
|
|
|
|
|
|
|
|
242
|
|
|
|
|
|
|
/* The final epoch */ |
243
|
|
|
|
|
|
|
double TFINAL; |
244
|
|
|
|
|
|
|
|
245
|
|
|
|
|
|
|
/* The body's current universal elements */ |
246
|
|
|
|
|
|
|
double UL[13]; |
247
|
|
|
|
|
|
|
|
248
|
|
|
|
|
|
|
/* Current reference epoch */ |
249
|
|
|
|
|
|
|
double T0; |
250
|
|
|
|
|
|
|
|
251
|
|
|
|
|
|
|
/* Timespan from latest orbit rectification to final epoch (days) */ |
252
|
|
|
|
|
|
|
double TSPAN; |
253
|
|
|
|
|
|
|
|
254
|
|
|
|
|
|
|
/* Time left to go before integration is complete */ |
255
|
|
|
|
|
|
|
double TLEFT; |
256
|
|
|
|
|
|
|
|
257
|
|
|
|
|
|
|
/* Time direction flag: +1=forwards, -1=backwards */ |
258
|
|
|
|
|
|
|
double FB; |
259
|
|
|
|
|
|
|
|
260
|
|
|
|
|
|
|
/* First-time flag */ |
261
|
|
|
|
|
|
|
int FIRST = 0; |
262
|
|
|
|
|
|
|
|
263
|
|
|
|
|
|
|
/* |
264
|
|
|
|
|
|
|
* The current perturbations |
265
|
|
|
|
|
|
|
*/ |
266
|
|
|
|
|
|
|
|
267
|
|
|
|
|
|
|
/* Epoch (days relative to current reference epoch) */ |
268
|
|
|
|
|
|
|
double RTN; |
269
|
|
|
|
|
|
|
/* Position (AU) */ |
270
|
|
|
|
|
|
|
double PERP[3]; |
271
|
|
|
|
|
|
|
/* Velocity (AU/d) */ |
272
|
|
|
|
|
|
|
double PERV[3]; |
273
|
|
|
|
|
|
|
/* Acceleration (AU/d/d) */ |
274
|
|
|
|
|
|
|
double PERA[3]; |
275
|
|
|
|
|
|
|
|
276
|
|
|
|
|
|
|
/* Length of current timestep (days), and half that */ |
277
|
|
|
|
|
|
|
double TS,HTS; |
278
|
|
|
|
|
|
|
|
279
|
|
|
|
|
|
|
/* Epoch of middle of timestep */ |
280
|
|
|
|
|
|
|
double T; |
281
|
|
|
|
|
|
|
|
282
|
|
|
|
|
|
|
/* Epoch of planetary mean elements */ |
283
|
|
|
|
|
|
|
double TPEL = 0.0; |
284
|
|
|
|
|
|
|
|
285
|
|
|
|
|
|
|
/* Planet number (1=Mercury, 2=Venus, 3=EMB...8=Neptune) */ |
286
|
|
|
|
|
|
|
int NP; |
287
|
|
|
|
|
|
|
|
288
|
|
|
|
|
|
|
/* Planetary universal orbital elements */ |
289
|
|
|
|
|
|
|
double UP[8][13]; |
290
|
|
|
|
|
|
|
|
291
|
|
|
|
|
|
|
/* Epoch of planetary state vectors */ |
292
|
|
|
|
|
|
|
double TPMO = 0.0; |
293
|
|
|
|
|
|
|
|
294
|
|
|
|
|
|
|
/* State vectors for the major planets (AU,AU/s) */ |
295
|
|
|
|
|
|
|
double PVIN[8][6]; |
296
|
|
|
|
|
|
|
|
297
|
|
|
|
|
|
|
/* Earth velocity and position vectors (AU,AU/s) */ |
298
|
|
|
|
|
|
|
double VB[3],PB[3],VH[3],PE[3]; |
299
|
|
|
|
|
|
|
|
300
|
|
|
|
|
|
|
/* Moon geocentric state vector (AU,AU/s) and position part */ |
301
|
|
|
|
|
|
|
double PVM[6],PM[3]; |
302
|
|
|
|
|
|
|
|
303
|
|
|
|
|
|
|
/* Date to J2000 de-precession matrix */ |
304
|
|
|
|
|
|
|
double PMAT[3][3]; |
305
|
|
|
|
|
|
|
|
306
|
|
|
|
|
|
|
/* |
307
|
|
|
|
|
|
|
* Correction terms for extrapolated major planet vectors |
308
|
|
|
|
|
|
|
*/ |
309
|
|
|
|
|
|
|
|
310
|
|
|
|
|
|
|
/* Sun-to-planet distances squared multiplied by 3 */ |
311
|
|
|
|
|
|
|
double R2X3[8]; |
312
|
|
|
|
|
|
|
/* Sunward acceleration terms, G/2R^3 */ |
313
|
|
|
|
|
|
|
double GC[8]; |
314
|
|
|
|
|
|
|
/* Tangential-to-circular correction factor */ |
315
|
|
|
|
|
|
|
double FC; |
316
|
|
|
|
|
|
|
/* Radial correction factor due to Sunwards acceleration */ |
317
|
|
|
|
|
|
|
double FG; |
318
|
|
|
|
|
|
|
|
319
|
|
|
|
|
|
|
/* The body's unperturbed and perturbed state vectors (AU,AU/s) */ |
320
|
|
|
|
|
|
|
double PV0[6],PV[6]; |
321
|
|
|
|
|
|
|
|
322
|
|
|
|
|
|
|
/* The body's perturbed and unperturbed heliocentric distances (AU) cubed */ |
323
|
|
|
|
|
|
|
double R03,R3; |
324
|
|
|
|
|
|
|
|
325
|
|
|
|
|
|
|
/* The perturbating accelerations, indirect and direct */ |
326
|
|
|
|
|
|
|
double FI[3],FD[3]; |
327
|
|
|
|
|
|
|
|
328
|
|
|
|
|
|
|
/* Sun-to-planet vector, and distance cubed */ |
329
|
|
|
|
|
|
|
double RHO[3],RHO3; |
330
|
|
|
|
|
|
|
|
331
|
|
|
|
|
|
|
/* Body-to-planet vector, and distance cubed */ |
332
|
|
|
|
|
|
|
double DELTA[3],DELTA3; |
333
|
|
|
|
|
|
|
|
334
|
|
|
|
|
|
|
/* Miscellaneous */ |
335
|
|
|
|
|
|
|
int I,J; |
336
|
|
|
|
|
|
|
double R2,W,DT,DT2,R,FT; |
337
|
|
|
|
|
|
|
int NE; |
338
|
|
|
|
|
|
|
|
339
|
|
|
|
|
|
|
/* Planetary inverse masses, Mercury through Neptune then Earth and Moon */ |
340
|
2
|
|
|
|
|
|
const double AMAS[10] = { |
341
|
|
|
|
|
|
|
6023600., 408523.5, 328900.5, 3098710., |
342
|
|
|
|
|
|
|
1047.355, 3498.5, 22869., 19314., |
343
|
|
|
|
|
|
|
332946.038, 27068709. |
344
|
|
|
|
|
|
|
}; |
345
|
|
|
|
|
|
|
|
346
|
|
|
|
|
|
|
/* Preset the status to OK. */ |
347
|
2
|
|
|
|
|
|
*jstat = 0; |
348
|
|
|
|
|
|
|
|
349
|
|
|
|
|
|
|
/* Copy the final epoch. */ |
350
|
|
|
|
|
|
|
TFINAL = date; |
351
|
|
|
|
|
|
|
|
352
|
|
|
|
|
|
|
/* Copy the elements (which will be periodically updated). */ |
353
|
28
|
100
|
|
|
|
|
for (I=0; I<13; I++) { |
354
|
26
|
|
|
|
|
|
UL[I] = u[I]; |
355
|
|
|
|
|
|
|
} |
356
|
|
|
|
|
|
|
|
357
|
|
|
|
|
|
|
/* Initialize the working reference epoch. */ |
358
|
2
|
|
|
|
|
|
T0=UL[2]; |
359
|
|
|
|
|
|
|
|
360
|
|
|
|
|
|
|
/* Total timespan (days) and hence time left. */ |
361
|
2
|
|
|
|
|
|
TSPAN = TFINAL-T0; |
362
|
|
|
|
|
|
|
TLEFT = TSPAN; |
363
|
|
|
|
|
|
|
|
364
|
|
|
|
|
|
|
/* Warn if excessive. */ |
365
|
2
|
50
|
|
|
|
|
if (fabs(TSPAN) > 36525.0) *jstat=101; |
366
|
|
|
|
|
|
|
|
367
|
|
|
|
|
|
|
/* Time direction: +1 for forwards, -1 for backwards. */ |
368
|
2
|
|
|
|
|
|
FB = COPYSIGN(1.0,TSPAN); |
369
|
|
|
|
|
|
|
|
370
|
|
|
|
|
|
|
/* Initialize relative epoch for start of current timestep. */ |
371
|
|
|
|
|
|
|
RTN = 0.0; |
372
|
|
|
|
|
|
|
|
373
|
|
|
|
|
|
|
/* Reset the perturbations (position, velocity, acceleration). */ |
374
|
8
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
375
|
6
|
|
|
|
|
|
PERP[I] = 0.0; |
376
|
6
|
|
|
|
|
|
PERV[I] = 0.0; |
377
|
6
|
|
|
|
|
|
PERA[I] = 0.0; |
378
|
|
|
|
|
|
|
} |
379
|
|
|
|
|
|
|
|
380
|
|
|
|
|
|
|
/* Set "first iteration" flag. */ |
381
|
|
|
|
|
|
|
FIRST = 1; |
382
|
|
|
|
|
|
|
|
383
|
|
|
|
|
|
|
/* Step through the time left. */ |
384
|
670
|
100
|
|
|
|
|
while (FB*TLEFT > 0.0) { |
385
|
|
|
|
|
|
|
|
386
|
|
|
|
|
|
|
/* Magnitude of current acceleration due to planetary attractions. */ |
387
|
668
|
100
|
|
|
|
|
if (FIRST) { |
388
|
|
|
|
|
|
|
TS = TSMIN; |
389
|
|
|
|
|
|
|
} else { |
390
|
|
|
|
|
|
|
R2 = 0.0; |
391
|
2664
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
392
|
1998
|
|
|
|
|
|
W = FD[I]; |
393
|
1998
|
|
|
|
|
|
R2 = R2+W*W; |
394
|
|
|
|
|
|
|
} |
395
|
666
|
|
|
|
|
|
W = sqrt(R2); |
396
|
|
|
|
|
|
|
|
397
|
|
|
|
|
|
|
/* Use the acceleration to decide how big a timestep can be tolerated. */ |
398
|
666
|
50
|
|
|
|
|
if (W != 0.0) { |
399
|
666
|
50
|
|
|
|
|
TS = DMIN(TSMAX,DMAX(TSMIN,TSC/W)); |
|
|
100
|
|
|
|
|
|
|
|
50
|
|
|
|
|
|
400
|
|
|
|
|
|
|
} else { |
401
|
|
|
|
|
|
|
TS = TSMAX; |
402
|
|
|
|
|
|
|
} |
403
|
|
|
|
|
|
|
} |
404
|
668
|
|
|
|
|
|
TS = TS*FB; |
405
|
|
|
|
|
|
|
|
406
|
|
|
|
|
|
|
/* Override if final epoch is imminent. */ |
407
|
668
|
|
|
|
|
|
TLEFT = TSPAN-RTN; |
408
|
668
|
100
|
|
|
|
|
if (fabs(TS) > fabs(TLEFT)) TS=TLEFT; |
409
|
|
|
|
|
|
|
|
410
|
|
|
|
|
|
|
/* Epoch of middle of timestep. */ |
411
|
668
|
|
|
|
|
|
HTS = TS/2.0; |
412
|
668
|
|
|
|
|
|
T = T0+RTN+HTS; |
413
|
|
|
|
|
|
|
|
414
|
|
|
|
|
|
|
/* Is it time to recompute the major-planet elements? */ |
415
|
668
|
100
|
|
|
|
|
if (FIRST || fabs(T-TPEL)-AGEPEL >= TINY) { |
|
|
100
|
|
|
|
|
|
416
|
|
|
|
|
|
|
|
417
|
|
|
|
|
|
|
/* Yes: go forward in time by just under the maximum allowed. */ |
418
|
22
|
|
|
|
|
|
TPEL = T+FB*AGEPEL; |
419
|
|
|
|
|
|
|
|
420
|
|
|
|
|
|
|
/* Compute the state vector for the new epoch. */ |
421
|
198
|
100
|
|
|
|
|
for (NP=1; NP<=8; NP++) { |
422
|
176
|
|
|
|
|
|
palPlanet(TPEL,NP,PV,&J); |
423
|
|
|
|
|
|
|
|
424
|
|
|
|
|
|
|
/* Warning if remote epoch, abort if error. */ |
425
|
176
|
50
|
|
|
|
|
if (J == 1) { |
426
|
0
|
|
|
|
|
|
*jstat = 102; |
427
|
176
|
50
|
|
|
|
|
} else if (J != 0) { |
428
|
|
|
|
|
|
|
goto ABORT; |
429
|
|
|
|
|
|
|
} |
430
|
|
|
|
|
|
|
|
431
|
|
|
|
|
|
|
/* Transform the vector into universal elements. */ |
432
|
176
|
|
|
|
|
|
palPv2ue(PV,TPEL,0.0,&(UP[NP-1][0]),&J); |
433
|
176
|
50
|
|
|
|
|
if (J != 0) goto ABORT; |
434
|
|
|
|
|
|
|
} |
435
|
|
|
|
|
|
|
} |
436
|
|
|
|
|
|
|
|
437
|
|
|
|
|
|
|
/* Is it time to recompute the major-planet motions? */ |
438
|
668
|
100
|
|
|
|
|
if (FIRST || fabs(T-TPMO)-AGEPMO >= TINY) { |
|
|
100
|
|
|
|
|
|
439
|
|
|
|
|
|
|
|
440
|
|
|
|
|
|
|
/* Yes: look ahead. */ |
441
|
176
|
|
|
|
|
|
TPMO = T+FB*AGEPMO; |
442
|
|
|
|
|
|
|
|
443
|
|
|
|
|
|
|
/* Compute the motions of each planet (AU,AU/d). */ |
444
|
1584
|
100
|
|
|
|
|
for (NP=1; NP<=8; NP++) { |
445
|
|
|
|
|
|
|
|
446
|
|
|
|
|
|
|
/* The planet's position and velocity (AU,AU/s). */ |
447
|
1408
|
|
|
|
|
|
palUe2pv(TPMO,&(UP[NP-1][0]),&(PVIN[NP-1][0]),&J); |
448
|
1408
|
50
|
|
|
|
|
if (J != 0) goto ABORT; |
449
|
|
|
|
|
|
|
|
450
|
|
|
|
|
|
|
/* Scale velocity to AU/d. */ |
451
|
5632
|
100
|
|
|
|
|
for (J=3; J<6; J++) { |
452
|
4224
|
|
|
|
|
|
PVIN[NP-1][J] = PVIN[NP-1][J]*PAL__SPD; |
453
|
|
|
|
|
|
|
} |
454
|
|
|
|
|
|
|
|
455
|
|
|
|
|
|
|
/* Precompute also the extrapolation correction terms. */ |
456
|
|
|
|
|
|
|
R2 = 0.0; |
457
|
5632
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
458
|
4224
|
|
|
|
|
|
W = PVIN[NP-1][I]; |
459
|
4224
|
|
|
|
|
|
R2 = R2+W*W; |
460
|
|
|
|
|
|
|
} |
461
|
1408
|
|
|
|
|
|
R2X3[NP-1] = R2*3.0; |
462
|
1408
|
|
|
|
|
|
GC[NP-1] = GCON2/(2.0*R2*sqrt(R2)); |
463
|
|
|
|
|
|
|
} |
464
|
|
|
|
|
|
|
} |
465
|
|
|
|
|
|
|
|
466
|
|
|
|
|
|
|
/* Reset the first-time flag. */ |
467
|
|
|
|
|
|
|
FIRST = 0; |
468
|
|
|
|
|
|
|
|
469
|
|
|
|
|
|
|
/* Unperturbed motion of the body at middle of timestep (AU,AU/s). */ |
470
|
668
|
|
|
|
|
|
palUe2pv(T,UL,PV0,&J); |
471
|
668
|
50
|
|
|
|
|
if (J != 0) goto ABORT; |
472
|
|
|
|
|
|
|
|
473
|
|
|
|
|
|
|
/* Perturbed position of the body (AU) and heliocentric distance cubed. */ |
474
|
|
|
|
|
|
|
R2 = 0.0; |
475
|
2672
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
476
|
2004
|
|
|
|
|
|
W = PV0[I]+PERP[I]+(PERV[I]+PERA[I]*HTS/2.0)*HTS; |
477
|
2004
|
|
|
|
|
|
PV[I] = W; |
478
|
2004
|
|
|
|
|
|
R2 = R2+W*W; |
479
|
|
|
|
|
|
|
} |
480
|
668
|
|
|
|
|
|
R3 = R2*sqrt(R2); |
481
|
|
|
|
|
|
|
|
482
|
|
|
|
|
|
|
/* The body's unperturbed heliocentric distance cubed. */ |
483
|
|
|
|
|
|
|
R2 = 0.0; |
484
|
2672
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
485
|
2004
|
|
|
|
|
|
W = PV0[I]; |
486
|
2004
|
|
|
|
|
|
R2 = R2+W*W; |
487
|
|
|
|
|
|
|
} |
488
|
668
|
|
|
|
|
|
R03 = R2*sqrt(R2); |
489
|
|
|
|
|
|
|
|
490
|
|
|
|
|
|
|
/* Compute indirect and initialize direct parts of the perturbation. */ |
491
|
2672
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
492
|
2004
|
|
|
|
|
|
FI[I] = PV0[I]/R03-PV[I]/R3; |
493
|
2004
|
|
|
|
|
|
FD[I] = 0.0; |
494
|
|
|
|
|
|
|
} |
495
|
|
|
|
|
|
|
|
496
|
|
|
|
|
|
|
/* Ready to compute the direct planetary effects. */ |
497
|
|
|
|
|
|
|
|
498
|
|
|
|
|
|
|
/* Reset the "near-Earth" flag. */ |
499
|
|
|
|
|
|
|
NE = 0; |
500
|
|
|
|
|
|
|
|
501
|
|
|
|
|
|
|
/* Interval from state-vector epoch to middle of current timestep. */ |
502
|
668
|
|
|
|
|
|
DT = T-TPMO; |
503
|
668
|
|
|
|
|
|
DT2 = DT*DT; |
504
|
|
|
|
|
|
|
|
505
|
|
|
|
|
|
|
/* Planet by planet, including separate Earth and Moon. */ |
506
|
6680
|
100
|
|
|
|
|
for (NP=1; NP<10; NP++) { |
507
|
|
|
|
|
|
|
|
508
|
|
|
|
|
|
|
/* Which perturbing body? */ |
509
|
6012
|
100
|
|
|
|
|
if (NP <= 8) { |
510
|
|
|
|
|
|
|
|
511
|
|
|
|
|
|
|
/* Planet: compute the extrapolation in longitude (squared). */ |
512
|
|
|
|
|
|
|
R2 = 0.0; |
513
|
21376
|
100
|
|
|
|
|
for (J=3; J<6; J++) { |
514
|
16032
|
|
|
|
|
|
W = PVIN[NP-1][J]*DT; |
515
|
16032
|
|
|
|
|
|
R2 = R2+W*W; |
516
|
|
|
|
|
|
|
} |
517
|
|
|
|
|
|
|
|
518
|
|
|
|
|
|
|
/* Hence the tangential-to-circular correction factor. */ |
519
|
5344
|
|
|
|
|
|
FC = 1.0+R2/R2X3[NP-1]; |
520
|
|
|
|
|
|
|
|
521
|
|
|
|
|
|
|
/* The radial correction factor due to the inwards acceleration. */ |
522
|
5344
|
|
|
|
|
|
FG = 1.0-GC[NP-1]*DT2; |
523
|
|
|
|
|
|
|
|
524
|
|
|
|
|
|
|
/* Planet's position. */ |
525
|
21376
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
526
|
16032
|
|
|
|
|
|
RHO[I] = FG*(PVIN[NP-1][I]+FC*PVIN[NP-1][I+3]*DT); |
527
|
|
|
|
|
|
|
} |
528
|
|
|
|
|
|
|
|
529
|
668
|
100
|
|
|
|
|
} else if (NE) { |
530
|
|
|
|
|
|
|
|
531
|
|
|
|
|
|
|
/* Near-Earth and either Earth or Moon. */ |
532
|
|
|
|
|
|
|
|
533
|
15
|
50
|
|
|
|
|
if (NP == 9) { |
534
|
|
|
|
|
|
|
|
535
|
|
|
|
|
|
|
/* Earth: position. */ |
536
|
15
|
|
|
|
|
|
palEpv(T,PE,VH,PB,VB); |
537
|
60
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
538
|
45
|
|
|
|
|
|
RHO[I] = PE[I]; |
539
|
|
|
|
|
|
|
} |
540
|
|
|
|
|
|
|
|
541
|
|
|
|
|
|
|
} else { |
542
|
|
|
|
|
|
|
|
543
|
|
|
|
|
|
|
/* Moon: position. */ |
544
|
0
|
|
|
|
|
|
palPrec(palEpj(T),2000.0,PMAT); |
545
|
0
|
|
|
|
|
|
palDmoon(T,PVM); |
546
|
0
|
|
|
|
|
|
eraRxp(PMAT,PVM,PM); |
547
|
0
|
0
|
|
|
|
|
for (I=0; I<3; I++) { |
548
|
0
|
|
|
|
|
|
RHO[I] = PM[I]+PE[I]; |
549
|
|
|
|
|
|
|
} |
550
|
|
|
|
|
|
|
} |
551
|
|
|
|
|
|
|
} |
552
|
|
|
|
|
|
|
|
553
|
|
|
|
|
|
|
/* Proceed unless Earth or Moon and not the near-Earth case. */ |
554
|
6012
|
100
|
|
|
|
|
if (NP <= 8 || NE) { |
555
|
|
|
|
|
|
|
|
556
|
|
|
|
|
|
|
/* Heliocentric distance cubed. */ |
557
|
|
|
|
|
|
|
R2 = 0.0; |
558
|
21436
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
559
|
16077
|
|
|
|
|
|
W = RHO[I]; |
560
|
16077
|
|
|
|
|
|
R2 = R2+W*W; |
561
|
|
|
|
|
|
|
} |
562
|
5359
|
|
|
|
|
|
R = sqrt(R2); |
563
|
5359
|
|
|
|
|
|
RHO3 = R2*R; |
564
|
|
|
|
|
|
|
|
565
|
|
|
|
|
|
|
/* Body-to-planet vector, and distance. */ |
566
|
|
|
|
|
|
|
R2 = 0.0; |
567
|
21436
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
568
|
16077
|
|
|
|
|
|
W = RHO[I]-PV[I]; |
569
|
16077
|
|
|
|
|
|
DELTA[I] = W; |
570
|
16077
|
|
|
|
|
|
R2 = R2+W*W; |
571
|
|
|
|
|
|
|
} |
572
|
5359
|
|
|
|
|
|
R = sqrt(R2); |
573
|
|
|
|
|
|
|
|
574
|
|
|
|
|
|
|
/* If this is the EMB, set the near-Earth flag appropriately. */ |
575
|
5359
|
100
|
|
|
|
|
if (NP == 3 && R < RNE) NE = 1; |
|
|
100
|
|
|
|
|
|
576
|
|
|
|
|
|
|
|
577
|
|
|
|
|
|
|
/* Proceed unless EMB and this is the near-Earth case. */ |
578
|
5359
|
100
|
|
|
|
|
if ( ! (NE && NP == 3) ) { |
579
|
|
|
|
|
|
|
|
580
|
|
|
|
|
|
|
/* If too close, ignore this planet and set a warning. */ |
581
|
5344
|
50
|
|
|
|
|
if (R < COINC) { |
582
|
0
|
|
|
|
|
|
*jstat = NP; |
583
|
|
|
|
|
|
|
|
584
|
|
|
|
|
|
|
} else { |
585
|
|
|
|
|
|
|
|
586
|
|
|
|
|
|
|
/* Accumulate "direct" part of perturbation acceleration. */ |
587
|
5344
|
|
|
|
|
|
DELTA3 = R2*R; |
588
|
5344
|
|
|
|
|
|
W = AMAS[NP-1]; |
589
|
21376
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
590
|
16032
|
|
|
|
|
|
FD[I] = FD[I]+(DELTA[I]/DELTA3-RHO[I]/RHO3)/W; |
591
|
|
|
|
|
|
|
} |
592
|
|
|
|
|
|
|
} |
593
|
|
|
|
|
|
|
} |
594
|
|
|
|
|
|
|
} |
595
|
|
|
|
|
|
|
} |
596
|
|
|
|
|
|
|
|
597
|
|
|
|
|
|
|
/* Update the perturbations to the end of the timestep. */ |
598
|
668
|
|
|
|
|
|
RTN += TS; |
599
|
2672
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
600
|
2004
|
|
|
|
|
|
W = (FI[I]+FD[I])*GCON2; |
601
|
2004
|
|
|
|
|
|
FT = W*TS; |
602
|
2004
|
|
|
|
|
|
PERP[I] = PERP[I]+(PERV[I]+FT/2.0)*TS; |
603
|
2004
|
|
|
|
|
|
PERV[I] = PERV[I]+FT; |
604
|
2004
|
|
|
|
|
|
PERA[I] = W; |
605
|
|
|
|
|
|
|
} |
606
|
|
|
|
|
|
|
|
607
|
|
|
|
|
|
|
/* Time still to go. */ |
608
|
668
|
|
|
|
|
|
TLEFT = TSPAN-RTN; |
609
|
|
|
|
|
|
|
|
610
|
|
|
|
|
|
|
/* Is it either time to rectify the orbit or the last time through? */ |
611
|
668
|
100
|
|
|
|
|
if (fabs(RTN) >= AGEBEL || FB*TLEFT <= 0.0) { |
|
|
100
|
|
|
|
|
|
612
|
|
|
|
|
|
|
|
613
|
|
|
|
|
|
|
/* Yes: update to the end of the current timestep. */ |
614
|
22
|
|
|
|
|
|
T0 += RTN; |
615
|
|
|
|
|
|
|
RTN = 0.0; |
616
|
|
|
|
|
|
|
|
617
|
|
|
|
|
|
|
/* The body's unperturbed motion (AU,AU/s). */ |
618
|
22
|
|
|
|
|
|
palUe2pv(T0,UL,PV0,&J); |
619
|
22
|
50
|
|
|
|
|
if (J != 0) goto ABORT; |
620
|
|
|
|
|
|
|
|
621
|
|
|
|
|
|
|
/* Add and re-initialize the perturbations. */ |
622
|
88
|
100
|
|
|
|
|
for (I=0; I<3; I++) { |
623
|
66
|
|
|
|
|
|
J = I+3; |
624
|
66
|
|
|
|
|
|
PV[I] = PV0[I]+PERP[I]; |
625
|
66
|
|
|
|
|
|
PV[J] = PV0[J]+PERV[I]/PAL__SPD; |
626
|
66
|
|
|
|
|
|
PERP[I] = 0.0; |
627
|
66
|
|
|
|
|
|
PERV[I] = 0.0; |
628
|
66
|
|
|
|
|
|
PERA[I] = FD[I]*GCON2; |
629
|
|
|
|
|
|
|
} |
630
|
|
|
|
|
|
|
|
631
|
|
|
|
|
|
|
/* Use the position and velocity to set up new universal elements. */ |
632
|
22
|
|
|
|
|
|
palPv2ue(PV,T0,0.0,UL,&J); |
633
|
22
|
50
|
|
|
|
|
if (J != 0) goto ABORT; |
634
|
|
|
|
|
|
|
|
635
|
|
|
|
|
|
|
/* Adjust the timespan and time left. */ |
636
|
668
|
|
|
|
|
|
TSPAN = TFINAL-T0; |
637
|
|
|
|
|
|
|
TLEFT = TSPAN; |
638
|
|
|
|
|
|
|
} |
639
|
|
|
|
|
|
|
|
640
|
|
|
|
|
|
|
/* Next timestep. */ |
641
|
|
|
|
|
|
|
} |
642
|
|
|
|
|
|
|
|
643
|
|
|
|
|
|
|
/* Return the updated universal-element set. */ |
644
|
28
|
100
|
|
|
|
|
for (I=0; I<13; I++) { |
645
|
26
|
|
|
|
|
|
u[I] = UL[I]; |
646
|
|
|
|
|
|
|
} |
647
|
|
|
|
|
|
|
|
648
|
|
|
|
|
|
|
/* Finished. */ |
649
|
|
|
|
|
|
|
return; |
650
|
|
|
|
|
|
|
|
651
|
|
|
|
|
|
|
/* Miscellaneous numerical error. */ |
652
|
|
|
|
|
|
|
ABORT: |
653
|
0
|
|
|
|
|
|
*jstat = -1; |
654
|
0
|
|
|
|
|
|
return; |
655
|
|
|
|
|
|
|
} |