source: trunk/MagicSoft/slalib/pertue.c@ 8723

Last change on this file since 8723 was 731, checked in by tbretz, 25 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 16.9 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void 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}
Note: See TracBrowser for help on using the repository browser.