File Coverage

palsrc/palPertue.c
Criterion Covered Total %
statement 118 127 92.9
branch 85 100 85.0
condition n/a
subroutine n/a
pod n/a
total 203 227 89.4


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             }