source: trunk/FACT++/pal/palPertue.c

Last change on this file was 18347, checked in by tbretz, 9 years ago
File size: 19.9 KB
Line 
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 <config.h>
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 <math.h>
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
212void 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 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 *jstat = 0;
348
349 /* Copy the final epoch. */
350 TFINAL = date;
351
352 /* Copy the elements (which will be periodically updated). */
353 for (I=0; I<13; I++) {
354 UL[I] = u[I];
355 }
356
357/* Initialize the working reference epoch. */
358 T0=UL[2];
359
360 /* Total timespan (days) and hence time left. */
361 TSPAN = TFINAL-T0;
362 TLEFT = TSPAN;
363
364 /* Warn if excessive. */
365 if (fabs(TSPAN) > 36525.0) *jstat=101;
366
367 /* Time direction: +1 for forwards, -1 for backwards. */
368 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 for (I=0; I<3; I++) {
375 PERP[I] = 0.0;
376 PERV[I] = 0.0;
377 PERA[I] = 0.0;
378 }
379
380 /* Set "first iteration" flag. */
381 FIRST = 1;
382
383 /* Step through the time left. */
384 while (FB*TLEFT > 0.0) {
385
386 /* Magnitude of current acceleration due to planetary attractions. */
387 if (FIRST) {
388 TS = TSMIN;
389 } else {
390 R2 = 0.0;
391 for (I=0; I<3; I++) {
392 W = FD[I];
393 R2 = R2+W*W;
394 }
395 W = sqrt(R2);
396
397 /* Use the acceleration to decide how big a timestep can be tolerated. */
398 if (W != 0.0) {
399 TS = DMIN(TSMAX,DMAX(TSMIN,TSC/W));
400 } else {
401 TS = TSMAX;
402 }
403 }
404 TS = TS*FB;
405
406 /* Override if final epoch is imminent. */
407 TLEFT = TSPAN-RTN;
408 if (fabs(TS) > fabs(TLEFT)) TS=TLEFT;
409
410 /* Epoch of middle of timestep. */
411 HTS = TS/2.0;
412 T = T0+RTN+HTS;
413
414 /* Is it time to recompute the major-planet elements? */
415 if (FIRST || fabs(T-TPEL)-AGEPEL >= TINY) {
416
417 /* Yes: go forward in time by just under the maximum allowed. */
418 TPEL = T+FB*AGEPEL;
419
420 /* Compute the state vector for the new epoch. */
421 for (NP=1; NP<=8; NP++) {
422 palPlanet(TPEL,NP,PV,&J);
423
424 /* Warning if remote epoch, abort if error. */
425 if (J == 1) {
426 *jstat = 102;
427 } else if (J != 0) {
428 goto ABORT;
429 }
430
431 /* Transform the vector into universal elements. */
432 palPv2ue(PV,TPEL,0.0,&(UP[NP-1][0]),&J);
433 if (J != 0) goto ABORT;
434 }
435 }
436
437 /* Is it time to recompute the major-planet motions? */
438 if (FIRST || fabs(T-TPMO)-AGEPMO >= TINY) {
439
440 /* Yes: look ahead. */
441 TPMO = T+FB*AGEPMO;
442
443 /* Compute the motions of each planet (AU,AU/d). */
444 for (NP=1; NP<=8; NP++) {
445
446 /* The planet's position and velocity (AU,AU/s). */
447 palUe2pv(TPMO,&(UP[NP-1][0]),&(PVIN[NP-1][0]),&J);
448 if (J != 0) goto ABORT;
449
450 /* Scale velocity to AU/d. */
451 for (J=3; J<6; J++) {
452 PVIN[NP-1][J] = PVIN[NP-1][J]*PAL__SPD;
453 }
454
455 /* Precompute also the extrapolation correction terms. */
456 R2 = 0.0;
457 for (I=0; I<3; I++) {
458 W = PVIN[NP-1][I];
459 R2 = R2+W*W;
460 }
461 R2X3[NP-1] = R2*3.0;
462 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 palUe2pv(T,UL,PV0,&J);
471 if (J != 0) goto ABORT;
472
473 /* Perturbed position of the body (AU) and heliocentric distance cubed. */
474 R2 = 0.0;
475 for (I=0; I<3; I++) {
476 W = PV0[I]+PERP[I]+(PERV[I]+PERA[I]*HTS/2.0)*HTS;
477 PV[I] = W;
478 R2 = R2+W*W;
479 }
480 R3 = R2*sqrt(R2);
481
482 /* The body's unperturbed heliocentric distance cubed. */
483 R2 = 0.0;
484 for (I=0; I<3; I++) {
485 W = PV0[I];
486 R2 = R2+W*W;
487 }
488 R03 = R2*sqrt(R2);
489
490 /* Compute indirect and initialize direct parts of the perturbation. */
491 for (I=0; I<3; I++) {
492 FI[I] = PV0[I]/R03-PV[I]/R3;
493 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 DT = T-TPMO;
503 DT2 = DT*DT;
504
505 /* Planet by planet, including separate Earth and Moon. */
506 for (NP=1; NP<10; NP++) {
507
508 /* Which perturbing body? */
509 if (NP <= 8) {
510
511 /* Planet: compute the extrapolation in longitude (squared). */
512 R2 = 0.0;
513 for (J=3; J<6; J++) {
514 W = PVIN[NP-1][J]*DT;
515 R2 = R2+W*W;
516 }
517
518 /* Hence the tangential-to-circular correction factor. */
519 FC = 1.0+R2/R2X3[NP-1];
520
521 /* The radial correction factor due to the inwards acceleration. */
522 FG = 1.0-GC[NP-1]*DT2;
523
524 /* Planet's position. */
525 for (I=0; I<3; I++) {
526 RHO[I] = FG*(PVIN[NP-1][I]+FC*PVIN[NP-1][I+3]*DT);
527 }
528
529 } else if (NE) {
530
531 /* Near-Earth and either Earth or Moon. */
532
533 if (NP == 9) {
534
535 /* Earth: position. */
536 palEpv(T,PE,VH,PB,VB);
537 for (I=0; I<3; I++) {
538 RHO[I] = PE[I];
539 }
540
541 } else {
542
543 /* Moon: position. */
544 palPrec(palEpj(T),2000.0,PMAT);
545 palDmoon(T,PVM);
546 eraRxp(PMAT,PVM,PM);
547 for (I=0; I<3; I++) {
548 RHO[I] = PM[I]+PE[I];
549 }
550 }
551 }
552
553 /* Proceed unless Earth or Moon and not the near-Earth case. */
554 if (NP <= 8 || NE) {
555
556 /* Heliocentric distance cubed. */
557 R2 = 0.0;
558 for (I=0; I<3; I++) {
559 W = RHO[I];
560 R2 = R2+W*W;
561 }
562 R = sqrt(R2);
563 RHO3 = R2*R;
564
565 /* Body-to-planet vector, and distance. */
566 R2 = 0.0;
567 for (I=0; I<3; I++) {
568 W = RHO[I]-PV[I];
569 DELTA[I] = W;
570 R2 = R2+W*W;
571 }
572 R = sqrt(R2);
573
574 /* If this is the EMB, set the near-Earth flag appropriately. */
575 if (NP == 3 && R < RNE) NE = 1;
576
577 /* Proceed unless EMB and this is the near-Earth case. */
578 if ( ! (NE && NP == 3) ) {
579
580 /* If too close, ignore this planet and set a warning. */
581 if (R < COINC) {
582 *jstat = NP;
583
584 } else {
585
586 /* Accumulate "direct" part of perturbation acceleration. */
587 DELTA3 = R2*R;
588 W = AMAS[NP-1];
589 for (I=0; I<3; I++) {
590 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 RTN += TS;
599 for (I=0; I<3; I++) {
600 W = (FI[I]+FD[I])*GCON2;
601 FT = W*TS;
602 PERP[I] = PERP[I]+(PERV[I]+FT/2.0)*TS;
603 PERV[I] = PERV[I]+FT;
604 PERA[I] = W;
605 }
606
607 /* Time still to go. */
608 TLEFT = TSPAN-RTN;
609
610 /* Is it either time to rectify the orbit or the last time through? */
611 if (fabs(RTN) >= AGEBEL || FB*TLEFT <= 0.0) {
612
613 /* Yes: update to the end of the current timestep. */
614 T0 += RTN;
615 RTN = 0.0;
616
617 /* The body's unperturbed motion (AU,AU/s). */
618 palUe2pv(T0,UL,PV0,&J);
619 if (J != 0) goto ABORT;
620
621 /* Add and re-initialize the perturbations. */
622 for (I=0; I<3; I++) {
623 J = I+3;
624 PV[I] = PV0[I]+PERP[I];
625 PV[J] = PV0[J]+PERV[I]/PAL__SPD;
626 PERP[I] = 0.0;
627 PERV[I] = 0.0;
628 PERA[I] = FD[I]*GCON2;
629 }
630
631 /* Use the position and velocity to set up new universal elements. */
632 palPv2ue(PV,T0,0.0,UL,&J);
633 if (J != 0) goto ABORT;
634
635 /* Adjust the timespan and time left. */
636 TSPAN = TFINAL-T0;
637 TLEFT = TSPAN;
638 }
639
640 /* Next timestep. */
641 }
642
643 /* Return the updated universal-element set. */
644 for (I=0; I<13; I++) {
645 u[I] = UL[I];
646 }
647
648 /* Finished. */
649 return;
650
651 /* Miscellaneous numerical error. */
652 ABORT:
653 *jstat = -1;
654 return;
655}
Note: See TracBrowser for help on using the repository browser.