FACT++  1.0
palPertue.c
Go to the documentation of this file.
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 
212 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  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 }
void palPertue(double date, double u[13], int *jstat)
Definition: palPertue.c:212
void palPrec(double ep0, double ep1, double rmatp[3][3])
Definition: palPrec.c:76
double palEpj(double date)
Definition: palOne2One.c:1169
void palPv2ue(const double pv[6], double date, double pmass, double u[13], int *jstat)
Definition: palPv2ue.c:114
#define DMIN(A, B)
Definition: palmac.h:129
static const double PAL__GCON
Definition: palmac.h:117
void palDmoon(double date, double pv[6])
Definition: palDmoon.c:95
#define COPYSIGN(a, b)
Definition: palPertue.c:209
void palEpv(double date, double ph[3], double vh[3], double pb[3], double vb[3])
Definition: palEpv.c:79
void eraRxp(double r[3][3], double p[3], double rp[3])
Definition: rxp.c:3
#define DMAX(A, B)
Definition: palmac.h:126
void palPlanet(double date, int np, double pv[6], int *j)
Definition: palPlanet.c:82
#define NP
static const double PAL__SPD
Definition: palmac.h:102
void palUe2pv(double date, double u[13], double pv[], int *jstat)