1 | #include "slalib.h"
|
---|
2 | #include "slamac.h"
|
---|
3 | void slaPertue ( double date, double u[], int *jstat )
|
---|
4 | /*
|
---|
5 | ** - - - - - - - - - -
|
---|
6 | ** s l a P e r t u e
|
---|
7 | ** - - - - - - - - - -
|
---|
8 | **
|
---|
9 | ** Update the universal elements of an asteroid or comet by applying
|
---|
10 | ** planetary perturbations.
|
---|
11 | **
|
---|
12 | ** Given:
|
---|
13 | ** date double final epoch (TT MJD) for the updated elements
|
---|
14 | **
|
---|
15 | ** Given and returned:
|
---|
16 | **
|
---|
17 | ** u double[13] universal orbital elements (Note 1)
|
---|
18 | **
|
---|
19 | ** [0] combined mass (M+m)
|
---|
20 | ** [1] total energy of the orbit (alpha)
|
---|
21 | ** [2] reference (osculating) epoch (t0)
|
---|
22 | ** [3-5] position at reference epoch (r0)
|
---|
23 | ** [6-8] velocity at reference epoch (v0)
|
---|
24 | ** [9] heliocentric distance at reference epoch
|
---|
25 | ** [10] r0.v0
|
---|
26 | ** [11] date (t)
|
---|
27 | ** [12] universal eccentric anomaly (psi) of date, approx
|
---|
28 | **
|
---|
29 | ** Returned:
|
---|
30 | ** jstat int* status:
|
---|
31 | ** +102 = warning, distant epoch
|
---|
32 | ** +101 = warning, large timespan ( > 100 years)
|
---|
33 | ** +1 to +8 = coincident with major planet (Note 5)
|
---|
34 | ** 0 = OK
|
---|
35 | ** -1 = numerical error
|
---|
36 | **
|
---|
37 | ** Called: slaPlanet, slaUe2pv, slaPv2ue
|
---|
38 | **
|
---|
39 | ** Notes:
|
---|
40 | **
|
---|
41 | ** 1 The "universal" elements are those which define the orbit for the
|
---|
42 | ** purposes of the method of universal variables (see reference 2).
|
---|
43 | ** They consist of the combined mass of the two bodies, an epoch,
|
---|
44 | ** and the position and velocity vectors (arbitrary reference frame)
|
---|
45 | ** at that epoch. The parameter set used here includes also various
|
---|
46 | ** quantities that can, in fact, be derived from the other
|
---|
47 | ** information. This approach is taken to avoiding unnecessary
|
---|
48 | ** computation and loss of accuracy. The supplementary quantities
|
---|
49 | ** are (i) alpha, which is proportional to the total energy of the
|
---|
50 | ** orbit, (ii) the heliocentric distance at epoch, (iii) the
|
---|
51 | ** outwards component of the velocity at the given epoch, (iv) an
|
---|
52 | ** estimate of psi, the "universal eccentric anomaly" at a given
|
---|
53 | ** date and (v) that date.
|
---|
54 | **
|
---|
55 | ** 2 The universal elements are with respect to the J2000 equator and
|
---|
56 | ** equinox.
|
---|
57 | **
|
---|
58 | ** 3 The epochs date, u[2] and u[11] are all Modified Julian Dates
|
---|
59 | ** (JD-2400000.5).
|
---|
60 | **
|
---|
61 | ** 4 The algorithm is a simplified form of Encke's method. It takes as
|
---|
62 | ** a basis the unperturbed motion of the body, and numerically
|
---|
63 | ** integrates the perturbing accelerations from the major planets.
|
---|
64 | ** The expression used is essentially Sterne's 6.7-2 (reference 1).
|
---|
65 | ** Everhart and Pitkin (reference 2) suggest rectifying the orbit at
|
---|
66 | ** each integration step by propagating the new perturbed position
|
---|
67 | ** and velocity as the new universal variables. In the present
|
---|
68 | ** routine the orbit is rectified less frequently than this, in order
|
---|
69 | ** to gain a slight speed advantage. However, the rectification is
|
---|
70 | ** done directly in terms of position and velocity, as suggested by
|
---|
71 | ** Everhart and Pitkin, bypassing the use of conventional orbital
|
---|
72 | ** elements.
|
---|
73 | **
|
---|
74 | ** The f(q) part of the full Encke method is not used. The purpose
|
---|
75 | ** of this part is to avoid subtracting two nearly equal quantities
|
---|
76 | ** when calculating the "indirect member", which takes account of the
|
---|
77 | ** small change in the Sun's attraction due to the slightly displaced
|
---|
78 | ** position of the perturbed body. A simpler, direct calculation in
|
---|
79 | ** double precision proves to be faster and not significantly less
|
---|
80 | ** accurate.
|
---|
81 | **
|
---|
82 | ** Apart from employing a variable timestep, and occasionally
|
---|
83 | ** "rectifying the orbit" to keep the indirect member small, the
|
---|
84 | ** integration is done in a fairly straightforward way. The
|
---|
85 | ** acceleration estimated for the middle of the timestep is assumed
|
---|
86 | ** to apply throughout that timestep; it is also used in the
|
---|
87 | ** extrapolation of the perturbations to the middle of the next
|
---|
88 | ** timestep, to predict the new disturbed position. There is no
|
---|
89 | ** iteration within a timestep.
|
---|
90 | **
|
---|
91 | ** Measures are taken to reach a compromise between execution time
|
---|
92 | ** and accuracy. The starting-point is the goal of achieving
|
---|
93 | ** arcsecond accuracy for ordinary minor planets over a ten-year
|
---|
94 | ** timespan. This goal dictates how large the timesteps can be,
|
---|
95 | ** which in turn dictates how frequently the unperturbed motion has
|
---|
96 | ** to be recalculated from the osculating elements.
|
---|
97 | **
|
---|
98 | ** Within predetermined limits, the timestep for the numerical
|
---|
99 | ** integration is varied in length in inverse proportion to the
|
---|
100 | ** magnitude of the net acceleration on the body from the major
|
---|
101 | ** planets.
|
---|
102 | **
|
---|
103 | ** The numerical integration requires estimates of the major-planet
|
---|
104 | ** motions. Approximate positions for the major planets (Pluto
|
---|
105 | ** alone is omitted) are obtained from the routine slaPlanet. Two
|
---|
106 | ** levels of interpolation are used, to enhance speed without
|
---|
107 | ** significantly degrading accuracy. At a low frequency, the routine
|
---|
108 | ** slaPlanet is called to generate updated position+velocity "state
|
---|
109 | ** vectors". The only task remaining to be carried out at the full
|
---|
110 | ** frequency (i.e. at each integration step) is to use the state
|
---|
111 | ** vectors to extrapolate the planetary positions. In place of a
|
---|
112 | ** strictly linear extrapolation, some allowance is made for the
|
---|
113 | ** curvature of the orbit by scaling back the radius vector as the
|
---|
114 | ** linear extrapolation goes off at a tangent.
|
---|
115 | **
|
---|
116 | ** Various other approximations are made. For example, perturbations
|
---|
117 | ** by Pluto and the minor planets are neglected, relativistic effects
|
---|
118 | ** are not taken into account and the Earth-Moon system is treated as
|
---|
119 | ** a single body.
|
---|
120 | **
|
---|
121 | ** In the interests of simplicity, the background calculations for
|
---|
122 | ** the major planets are carried out en masse. The mean elements and
|
---|
123 | ** state vectors for all the planets are refreshed at the same time,
|
---|
124 | ** without regard for orbit curvature, mass or proximity.
|
---|
125 | **
|
---|
126 | ** 5 This routine is not intended to be used for major planets.
|
---|
127 | ** However, if major-planet elements are supplied, sensible results
|
---|
128 | ** will, in fact, be produced. This happens because the routine
|
---|
129 | ** checks the separation between the body and each of the planets and
|
---|
130 | ** interprets a suspiciously small value (0.001 AU) as an attempt to
|
---|
131 | ** apply the routine to the planet concerned. If this condition is
|
---|
132 | ** detected, the contribution from that planet is ignored, and the
|
---|
133 | ** status is set to the planet number (Mercury=1,...,Neptune=8) as a
|
---|
134 | ** warning.
|
---|
135 | **
|
---|
136 | ** References:
|
---|
137 | **
|
---|
138 | ** 1 Sterne, Theodore E., "An Introduction to Celestial Mechanics",
|
---|
139 | ** Interscience Publishers Inc., 1960. Section 6.7, p199.
|
---|
140 | **
|
---|
141 | ** 2 Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
|
---|
142 | **
|
---|
143 | ** Last revision: 18 March 1999
|
---|
144 | **
|
---|
145 | ** Copyright P.T.Wallace. All rights reserved.
|
---|
146 | */
|
---|
147 |
|
---|
148 | /* Coefficient relating timestep to perturbing force */
|
---|
149 | #define TSC 1e-4
|
---|
150 |
|
---|
151 | /* Minimum and maximum timestep (days) */
|
---|
152 | #define TSMIN 0.01
|
---|
153 | #define TSMAX 10.0
|
---|
154 |
|
---|
155 | /* Age limit for major-planet state vector (days) */
|
---|
156 | #define AGEPMO 5.0
|
---|
157 |
|
---|
158 | /* Age limit for major-planet mean elements (days) */
|
---|
159 | #define AGEPEL 50.0
|
---|
160 |
|
---|
161 | /* Margin for error when deciding whether to renew the planetary data */
|
---|
162 | #define TINY 1e-6
|
---|
163 |
|
---|
164 | /* Age limit for the body's osculating elements (before rectification) */
|
---|
165 | #define AGEBEL 100.0
|
---|
166 |
|
---|
167 | /* Gaussian gravitational constant (exact) and square */
|
---|
168 | #define GCON 0.01720209895
|
---|
169 | #define GCON2 (GCON*GCON)
|
---|
170 |
|
---|
171 | {
|
---|
172 |
|
---|
173 | /* The final epoch */
|
---|
174 | double tfinal;
|
---|
175 |
|
---|
176 | /* The body's current universal elements */
|
---|
177 | double ul[13];
|
---|
178 |
|
---|
179 | /* Current reference epoch */
|
---|
180 | double t0;
|
---|
181 |
|
---|
182 | /* Timespan from latest orbit rectification to final epoch (days) */
|
---|
183 | double tspan;
|
---|
184 |
|
---|
185 | /* Time left to go before integration is complete */
|
---|
186 | double tleft;
|
---|
187 |
|
---|
188 | /* Time direction flag: +1=forwards, -1=backwards */
|
---|
189 | double fb;
|
---|
190 |
|
---|
191 | /* First-time flag */
|
---|
192 | int first;
|
---|
193 |
|
---|
194 | /* The current perturbations */
|
---|
195 | double rtn, /* Epoch (days relative to current reference epoch) */
|
---|
196 | perp[3], /* Position (AU) */
|
---|
197 | perv[3], /* Velocity (AU/d) */
|
---|
198 | pera[3]; /* Acceleration (AU/d/d) */
|
---|
199 |
|
---|
200 | /* Length of current timestep (days), and half that */
|
---|
201 | double ts, hts;
|
---|
202 |
|
---|
203 | /* Epoch of middle of timestep */
|
---|
204 | double t;
|
---|
205 |
|
---|
206 | /* Epoch of planetary mean elements */
|
---|
207 | double tpel;
|
---|
208 |
|
---|
209 | /* Planet number (1=Mercury, 2=Venus, 3=EMB...8=Neptune) */
|
---|
210 | int np;
|
---|
211 |
|
---|
212 | /* Planetary universal orbital elements */
|
---|
213 | double up[8][13];
|
---|
214 |
|
---|
215 | /* Epoch of planetary state vectors */
|
---|
216 | double tpmo;
|
---|
217 |
|
---|
218 | /* State vectors for the major planets (AU,AU/s) */
|
---|
219 | double pvin[8][6];
|
---|
220 |
|
---|
221 | /* Correction terms for extrapolated major planet vectors */
|
---|
222 | double r2x3[8], /* Sun-to-planet distances squared multiplied by 3 */
|
---|
223 | gc[8], /* Sunward acceleration terms, G/2R^3 */
|
---|
224 | fc, /* Tangential-to-circular correction factor */
|
---|
225 | fg; /* Radial correction factor due to Sunwards acceleration */
|
---|
226 |
|
---|
227 | /* The body's unperturbed and perturbed state vectors (AU,AU/s) */
|
---|
228 | double pv0[6], pv[6];
|
---|
229 |
|
---|
230 | /* The body's perturbed and unperturbed heliocentric distances (AU) cubed */
|
---|
231 | double r03, r3;
|
---|
232 |
|
---|
233 | /* The perturbating accelerations, indirect and direct */
|
---|
234 | double fi[3], fd[3];
|
---|
235 |
|
---|
236 | /* Sun-to-planet vector, and distance cubed */
|
---|
237 | double rho[3], rho3;
|
---|
238 |
|
---|
239 | /* Body-to-planet vector, and distance cubed */
|
---|
240 | double delta[3], delta3;
|
---|
241 |
|
---|
242 | /* Miscellaneous */
|
---|
243 | int i, j, npm1;
|
---|
244 | double r2, w, dt, dt2, ft;
|
---|
245 |
|
---|
246 | /* Planetary inverse masses, Mercury through Neptune */
|
---|
247 | static double amas[] = {
|
---|
248 | 6023600.0,
|
---|
249 | 408523.5,
|
---|
250 | 328900.5,
|
---|
251 | 3098710.0,
|
---|
252 | 1047.355,
|
---|
253 | 3498.5,
|
---|
254 | 22869.0,
|
---|
255 | 19314.0
|
---|
256 | };
|
---|
257 |
|
---|
258 |
|
---|
259 |
|
---|
260 | /* Preset the status to OK. */
|
---|
261 | *jstat = 0;
|
---|
262 |
|
---|
263 | /* Copy the final epoch. */
|
---|
264 | tfinal = date;
|
---|
265 |
|
---|
266 | /* Copy the elements (which will be periodically updated). */
|
---|
267 | for ( i = 0; i < 13; i++ ) {
|
---|
268 | ul[i] = u[i];
|
---|
269 | }
|
---|
270 |
|
---|
271 | /* Initialize the working reference epoch. */
|
---|
272 | t0 = ul[2];
|
---|
273 |
|
---|
274 | /* Total timespan (days) and hence time left. */
|
---|
275 | tspan = tfinal - t0;
|
---|
276 | tleft = tspan;
|
---|
277 |
|
---|
278 | /* Warn if excessive. */
|
---|
279 | if ( fabs ( tspan ) > 36525.0 ) *jstat = 101;
|
---|
280 |
|
---|
281 | /* Time direction: +1 for forwards, -1 for backwards. */
|
---|
282 | fb = dsign ( 1.0, tspan );
|
---|
283 |
|
---|
284 | /* Initialize relative epoch for start of current timestep. */
|
---|
285 | rtn = 0.0;
|
---|
286 |
|
---|
287 | /* Reset the perturbations (position, velocity, acceleration). */
|
---|
288 | for ( i = 0; i < 3; i++ ) {
|
---|
289 | perp[i] = 0.0;
|
---|
290 | perv[i] = 0.0;
|
---|
291 | pera[i] = 0.0;
|
---|
292 | }
|
---|
293 |
|
---|
294 | /* Set "first iteration" flag. */
|
---|
295 | first = TRUE;
|
---|
296 |
|
---|
297 | /* Step through the time left. */
|
---|
298 | while ( fb * tleft > 0.0 ) {
|
---|
299 |
|
---|
300 | /* Magnitude of current acceleration due to planetary attractions. */
|
---|
301 | if ( first ) {
|
---|
302 | ts = TSMIN;
|
---|
303 | } else {
|
---|
304 | r2 = 0.0;
|
---|
305 | for ( i = 0; i < 3; i++ ) {
|
---|
306 | w = fd[i];
|
---|
307 | r2 += w * w;
|
---|
308 | }
|
---|
309 | w = sqrt ( r2 );
|
---|
310 |
|
---|
311 | /* Use the acceleration to decide how big a timestep can be tolerated. */
|
---|
312 | if ( w != 0.0 ) {
|
---|
313 | ts = TSC / w;
|
---|
314 | if ( ts > TSMAX ) {
|
---|
315 | ts = TSMAX;
|
---|
316 | } else if ( ts < TSMIN ) {
|
---|
317 | ts = TSMIN;
|
---|
318 | }
|
---|
319 | } else {
|
---|
320 | ts = TSMAX;
|
---|
321 | }
|
---|
322 | }
|
---|
323 | ts *= fb;
|
---|
324 |
|
---|
325 | /* Override if final epoch is imminent. */
|
---|
326 | tleft = tspan - rtn;
|
---|
327 | if ( fabs ( ts ) > fabs ( tleft ) ) ts = tleft;
|
---|
328 |
|
---|
329 | /* Epoch of middle of timestep. */
|
---|
330 | hts = ts / 2.0;
|
---|
331 | t = t0 + rtn + hts;
|
---|
332 |
|
---|
333 | /* Is it time to recompute the major-planet elements? */
|
---|
334 | if ( first || ( fabs ( t - tpel ) - AGEPEL ) >= TINY ) {
|
---|
335 |
|
---|
336 | /* Yes: go forward in time by just under the maximum allowed. */
|
---|
337 | tpel = t + fb * AGEPEL;
|
---|
338 |
|
---|
339 | /* Compute the state vector for the new epoch. */
|
---|
340 | for ( np = 1; np <= 8; np++ ) {
|
---|
341 | npm1 = np - 1;
|
---|
342 |
|
---|
343 | slaPlanet ( tpel, np, pv, &j );
|
---|
344 |
|
---|
345 | /* Warning if remote epoch, abort if error. */
|
---|
346 | if ( j == 1 ) {
|
---|
347 | *jstat = 102;
|
---|
348 | } else if ( j ) {
|
---|
349 | *jstat = -1;
|
---|
350 | return;
|
---|
351 | }
|
---|
352 |
|
---|
353 | /* Transform the vector into universal elements. */
|
---|
354 | slaPv2ue ( pv, tpel, 0.0, up[npm1], &j );
|
---|
355 | if ( j ) {
|
---|
356 | *jstat = -1;
|
---|
357 | return;
|
---|
358 | }
|
---|
359 | }
|
---|
360 | }
|
---|
361 |
|
---|
362 | /* Is it time to recompute the major-planet motions? */
|
---|
363 | if ( first || ( fabs ( t - tpmo ) - AGEPMO ) >= TINY ) {
|
---|
364 |
|
---|
365 | /* Yes: look ahead. */
|
---|
366 | tpmo = t + fb * AGEPMO;
|
---|
367 |
|
---|
368 | /* Compute the motions of each planet (AU,AU/d). */
|
---|
369 | for ( np = 1; np <= 8; np++ ) {
|
---|
370 | npm1 = np - 1;
|
---|
371 |
|
---|
372 | /* The planet's position and velocity (AU,AU/s). */
|
---|
373 | slaUe2pv ( tpmo, up[npm1], pvin[npm1], &j );
|
---|
374 | if ( j ) {
|
---|
375 | *jstat = -1;
|
---|
376 | return;
|
---|
377 | }
|
---|
378 |
|
---|
379 | /* Scale velocity to AU/d. */
|
---|
380 | for ( j = 3; j < 6; j++ ) {
|
---|
381 | pvin[npm1][j] *= 86400.0;
|
---|
382 | }
|
---|
383 |
|
---|
384 | /* Precompute also the extrapolation correction terms. */
|
---|
385 | r2 = 0.0;
|
---|
386 | for ( i = 0; i < 3; i++ ) {
|
---|
387 | w = pvin[npm1][i];
|
---|
388 | r2 += w * w;
|
---|
389 | }
|
---|
390 | r2x3[npm1] = r2 * 3.0;
|
---|
391 | gc[npm1] = GCON2 / ( 2.0 * r2 * sqrt ( r2 ) );
|
---|
392 | }
|
---|
393 | }
|
---|
394 |
|
---|
395 | /* Reset the first-time flag. */
|
---|
396 | first = FALSE;
|
---|
397 |
|
---|
398 | /* Unperturbed motion of the body at middle of timestep (AU,AU/s). */
|
---|
399 | slaUe2pv ( t, ul, pv0, &j );
|
---|
400 | if ( j ) {
|
---|
401 | *jstat = -1;
|
---|
402 | return;
|
---|
403 | }
|
---|
404 |
|
---|
405 | /* Perturbed position of the body (AU) and heliocentric distance cubed. */
|
---|
406 | r2 = 0.0;
|
---|
407 | for ( i = 0; i < 3; i++ ) {
|
---|
408 | w = pv0[i] + perp[i] + ( perv[i] + pera[i] * hts / 2.0 ) * hts;
|
---|
409 | pv[i] = w;
|
---|
410 | r2 += w * w;
|
---|
411 | }
|
---|
412 | r3 = r2 * sqrt ( r2 );
|
---|
413 |
|
---|
414 | /* The body's unperturbed heliocentric distance cubed. */
|
---|
415 | r2 = 0.0;
|
---|
416 | for ( i = 0; i < 3; i++ ) {
|
---|
417 | w = pv0[i];
|
---|
418 | r2 += w * w;
|
---|
419 | }
|
---|
420 | r03 = r2 * sqrt ( r2 );
|
---|
421 |
|
---|
422 | /* Compute indirect and initialize direct parts of the perturbation. */
|
---|
423 | for ( i = 0; i < 3; i++ ) {
|
---|
424 | fi[i] = pv0[i] / r03 - pv[i] / r3;
|
---|
425 | fd[i] = 0.0;
|
---|
426 | }
|
---|
427 |
|
---|
428 | /* Ready to compute the direct planetary effects. */
|
---|
429 |
|
---|
430 | /* Interval from state-vector epoch to middle of current timestep. */
|
---|
431 | dt = t - tpmo;
|
---|
432 | dt2 = dt * dt;
|
---|
433 |
|
---|
434 | /* Planet by planet. */
|
---|
435 | for ( np = 1; np <= 8; np++ ) {
|
---|
436 | npm1 = np - 1;
|
---|
437 |
|
---|
438 | /* First compute the extrapolation in longitude (squared). */
|
---|
439 | r2 = 0.0;
|
---|
440 | for ( j = 3; j < 6; j++ ) {
|
---|
441 | w = pvin[npm1][j] * dt;
|
---|
442 | r2 += w * w;
|
---|
443 | }
|
---|
444 |
|
---|
445 | /* Hence the tangential-to-circular correction factor. */
|
---|
446 | fc = 1.0 + r2 / r2x3[npm1];
|
---|
447 |
|
---|
448 | /* The radial correction factor due to the inwards acceleration. */
|
---|
449 | fg = 1.0 - gc[npm1] * dt2;
|
---|
450 |
|
---|
451 | /* Planet's position, and heliocentric distance cubed. */
|
---|
452 | r2 = 0.0;
|
---|
453 | for ( i = 0; i < 3; i++ ) {
|
---|
454 | w = fg * ( pvin[npm1][i] + fc * pvin[npm1][i+3] * dt );
|
---|
455 | rho[i] = w;
|
---|
456 | r2 += w * w;
|
---|
457 | }
|
---|
458 | rho3 = r2 * sqrt ( r2 );
|
---|
459 |
|
---|
460 | /* Body-to-planet vector, and distance cubed. */
|
---|
461 | r2 = 0.0;
|
---|
462 | for ( i = 0; i < 3; i++ ) {
|
---|
463 | w = rho[i] - pv[i];
|
---|
464 | delta[i] = w;
|
---|
465 | r2 += w * w;
|
---|
466 | }
|
---|
467 | delta3 = r2 * sqrt ( r2 );
|
---|
468 |
|
---|
469 | /* If too close, ignore this planet and set a warning. */
|
---|
470 | if ( r2 < 1e-6 ) {
|
---|
471 | *jstat = np;
|
---|
472 | } else {
|
---|
473 |
|
---|
474 | /* Accumulate "direct" part of perturbation acceleration. */
|
---|
475 | w = amas[npm1];
|
---|
476 | for ( i = 0; i < 3; i++ ) {
|
---|
477 | fd[i] += ( delta[i] / delta3 - rho[i] / rho3 ) / w;
|
---|
478 | }
|
---|
479 | }
|
---|
480 | }
|
---|
481 |
|
---|
482 | /* Update the perturbations to the end of the timestep. */
|
---|
483 | rtn = rtn + ts;
|
---|
484 | for ( i = 0; i < 3; i++ ) {
|
---|
485 | w = ( fi[i] + fd[i] ) * GCON2;
|
---|
486 | ft = w * ts;
|
---|
487 | perp[i] += ( perv[i] + ft / 2.0 ) * ts;
|
---|
488 | perv[i] += ft;
|
---|
489 | pera[i] = w;
|
---|
490 | }
|
---|
491 |
|
---|
492 | /* Time still to go. */
|
---|
493 | tleft = tspan - rtn;
|
---|
494 |
|
---|
495 | /* Is it either time to rectify the orbit or the last time through? */
|
---|
496 | if ( fabs ( rtn ) >= AGEBEL || ( fb * tleft ) <= 0.0 ) {
|
---|
497 |
|
---|
498 | /* Yes: update to the end of the current timestep. */
|
---|
499 | t0 += rtn;
|
---|
500 | rtn = 0.0;
|
---|
501 |
|
---|
502 | /* The body's unperturbed motion (AU,AU/s). */
|
---|
503 | slaUe2pv ( t0, ul, pv0, &j );
|
---|
504 | if ( j ) {
|
---|
505 | *jstat = -1;
|
---|
506 | return;
|
---|
507 | }
|
---|
508 |
|
---|
509 | /* Add and re-initialize the perturbations. */
|
---|
510 | for ( i = 0; i < 3; i++ ) {
|
---|
511 | j = i + 3;
|
---|
512 | pv[i] = pv0[i] + perp[i];
|
---|
513 | pv[j] = pv0[j] + perv[i] / 86400.0;
|
---|
514 | perp[i] = 0.0;
|
---|
515 | perv[i] = 0.0;
|
---|
516 | pera[i] = fd[i] * GCON2;
|
---|
517 | }
|
---|
518 |
|
---|
519 | /* Use the position and velocity to set up new universal elements. */
|
---|
520 | slaPv2ue ( pv, t0, 0.0, ul, &j );
|
---|
521 | if ( j ) {
|
---|
522 | *jstat = -1;
|
---|
523 | return;
|
---|
524 | }
|
---|
525 |
|
---|
526 | /* Adjust the timespan and time left. */
|
---|
527 | tspan = tfinal - t0;
|
---|
528 | tleft = tspan;
|
---|
529 | }
|
---|
530 |
|
---|
531 | /* Next timestep. */
|
---|
532 | }
|
---|
533 |
|
---|
534 | /* Return the updated universal-element set. */
|
---|
535 | for ( i = 0; i < 13; i++ ) {
|
---|
536 | u[i] = ul[i];
|
---|
537 | }
|
---|
538 |
|
---|
539 | }
|
---|