palPertue.c

/*
*+
*  Name:
*     palPertue

*  Purpose:
*     Update the universal elements by applying planetary perturbations

*  Language:
*     Starlink ANSI C

*  Type of Module:
*     Library routine

*  Invocation:
*     void palPertue( double date, double u[13], int *jstat );

*  Arguments:
*     date = double (Given)
*        Final epoch (TT MJD) for the update elements.
*     u = const double [13] (Given & Returned)
*        Universal orbital elements (Note 1)
*            (0)  combined mass (M+m)
*            (1)  total energy of the orbit (alpha)
*            (2)  reference (osculating) epoch (t0)
*          (3-5)  position at reference epoch (r0)
*          (6-8)  velocity at reference epoch (v0)
*            (9)  heliocentric distance at reference epoch
*           (10)  r0.v0
*           (11)  date (t)
*           (12)  universal eccentric anomaly (psi) of date, approx
*     jstat = int * (Returned)
*        status:
*                   +102 = warning, distant epoch
*                   +101 = warning, large timespan ( > 100 years)
*              +1 to +10 = coincident with major planet (Note 5)
*                      0 = OK
*                     -1 = numerical error

*  Description:
*     Update the universal elements of an asteroid or comet by applying
*     planetary perturbations.

*  Authors:
*     PTW: Pat Wallace (STFC)
*     TIMJ: Tim Jenness (JAC, Hawaii)
*     {enter_new_authors_here}

*  Notes:
*     - The "universal" elements are those which define the orbit for the
*       purposes of the method of universal variables (see reference 2).
*       They consist of the combined mass of the two bodies, an epoch,
*       and the position and velocity vectors (arbitrary reference frame)
*       at that epoch.  The parameter set used here includes also various
*       quantities that can, in fact, be derived from the other
*       information.  This approach is taken to avoiding unnecessary
*       computation and loss of accuracy.  The supplementary quantities
*       are (i) alpha, which is proportional to the total energy of the
*       orbit, (ii) the heliocentric distance at epoch, (iii) the
*       outwards component of the velocity at the given epoch, (iv) an
*       estimate of psi, the "universal eccentric anomaly" at a given
*       date and (v) that date.
*     - The universal elements are with respect to the J2000 equator and
*       equinox.
*     - The epochs DATE, U(3) and U(12) are all Modified Julian Dates
*       (JD-2400000.5).
*     - The algorithm is a simplified form of Encke's method.  It takes as
*       a basis the unperturbed motion of the body, and numerically
*       integrates the perturbing accelerations from the major planets.
*       The expression used is essentially Sterne's 6.7-2 (reference 1).
*       Everhart and Pitkin (reference 2) suggest rectifying the orbit at
*       each integration step by propagating the new perturbed position
*       and velocity as the new universal variables.  In the present
*       routine the orbit is rectified less frequently than this, in order
*       to gain a slight speed advantage.  However, the rectification is
*       done directly in terms of position and velocity, as suggested by
*       Everhart and Pitkin, bypassing the use of conventional orbital
*       elements.
*
*       The f(q) part of the full Encke method is not used.  The purpose
*       of this part is to avoid subtracting two nearly equal quantities
*       when calculating the "indirect member", which takes account of the
*       small change in the Sun's attraction due to the slightly displaced
*       position of the perturbed body.  A simpler, direct calculation in
*       double precision proves to be faster and not significantly less
*       accurate.
*
*       Apart from employing a variable timestep, and occasionally
*       "rectifying the orbit" to keep the indirect member small, the
*       integration is done in a fairly straightforward way.  The
*       acceleration estimated for the middle of the timestep is assumed
*       to apply throughout that timestep;  it is also used in the
*       extrapolation of the perturbations to the middle of the next
*       timestep, to predict the new disturbed position.  There is no
*       iteration within a timestep.
*
*       Measures are taken to reach a compromise between execution time
*       and accuracy.  The starting-point is the goal of achieving
*       arcsecond accuracy for ordinary minor planets over a ten-year
*       timespan.  This goal dictates how large the timesteps can be,
*       which in turn dictates how frequently the unperturbed motion has
*       to be recalculated from the osculating elements.
*
*       Within predetermined limits, the timestep for the numerical
*       integration is varied in length in inverse proportion to the
*       magnitude of the net acceleration on the body from the major
*       planets.
*
*       The numerical integration requires estimates of the major-planet
*       motions.  Approximate positions for the major planets (Pluto
*       alone is omitted) are obtained from the routine palPlanet.  Two
*       levels of interpolation are used, to enhance speed without
*       significantly degrading accuracy.  At a low frequency, the routine
*       palPlanet is called to generate updated position+velocity "state
*       vectors".  The only task remaining to be carried out at the full
*       frequency (i.e. at each integration step) is to use the state
*       vectors to extrapolate the planetary positions.  In place of a
*       strictly linear extrapolation, some allowance is made for the
*       curvature of the orbit by scaling back the radius vector as the
*       linear extrapolation goes off at a tangent.
*
*       Various other approximations are made.  For example, perturbations
*       by Pluto and the minor planets are neglected and relativistic
*       effects are not taken into account.
*
*       In the interests of simplicity, the background calculations for
*       the major planets are carried out en masse.  The mean elements and
*       state vectors for all the planets are refreshed at the same time,
*       without regard for orbit curvature, mass or proximity.
*
*       The Earth-Moon system is treated as a single body when the body is
*       distant but as separate bodies when closer to the EMB than the
*       parameter RNE, which incurs a time penalty but improves accuracy
*       for near-Earth objects.
*
*     - This routine is not intended to be used for major planets.
*       However, if major-planet elements are supplied, sensible results
*       will, in fact, be produced.  This happens because the routine
*       checks the separation between the body and each of the planets and
*       interprets a suspiciously small value (0.001 AU) as an attempt to
*       apply the routine to the planet concerned.  If this condition is
*       detected, the contribution from that planet is ignored, and the
*       status is set to the planet number (1-10 = Mercury, Venus, EMB,
*       Mars, Jupiter, Saturn, Uranus, Neptune, Earth, Moon) as a warning.

*  See Also:
*     - Sterne, Theodore E., "An Introduction to Celestial Mechanics",
*       Interscience Publishers Inc., 1960.  Section 6.7, p199.
*     - Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.

*  History:
*     2012-03-12 (TIMJ):
*        Initial version direct conversion of SLA/F.
*        Adapted with permission from the Fortran SLALIB library.
*     2012-06-21 (TIMJ):
*        Support a lack of copysign() function.
*     2012-06-22 (TIMJ):
*        Check __STDC_VERSION__
*     {enter_further_changes_here}

*  Copyright:
*     Copyright (C) 2004 Patrick T. Wallace
*     Copyright (C) 2012 Science and Technology Facilities Council.
*     All Rights Reserved.

*  Licence:
*     This program is free software; you can redistribute it and/or
*     modify it under the terms of the GNU General Public License as
*     published by the Free Software Foundation; either version 3 of
*     the License, or (at your option) any later version.
*
*     This program is distributed in the hope that it will be
*     useful, but WITHOUT ANY WARRANTY; without even the implied
*     warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
*     PURPOSE. See the GNU General Public License for more details.
*
*     You should have received a copy of the GNU General Public License
*     along with this program; if not, write to the Free Software
*     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
*     MA 02110-1301, USA.

*  Bugs:
*     {note_any_bugs_here}
*-
*/

/* Use the config file if we have one, else look at
   compiler defines to see if we have C99 */
#if HAVE_CONFIG_H
#include <config.h>
#else
#ifdef __STDC_VERSION__
#  if (__STDC_VERSION__ >= 199901L)
#    define HAVE_COPYSIGN 1
#  endif
#endif
#endif

#include <math.h>

#include "pal.h"
#include "palmac.h"
#include "pal1sofa.h"

/* copysign is C99 */
#if HAVE_COPYSIGN
# define COPYSIGN copysign
#else
# define COPYSIGN(a,b) DSIGN(a,b)
#endif

void palPertue( double date, double u[13], int *jstat ) {

  /*  Distance from EMB at which Earth and Moon are treated separately */
  const double RNE=1.0;

  /*  Coincidence with major planet distance */
  const double COINC=0.0001;

  /*  Coefficient relating timestep to perturbing force */
  const double TSC=1e-4;

  /*  Minimum and maximum timestep (days) */
  const double TSMIN = 0.01;
  const double TSMAX = 10.0;

  /*  Age limit for major-planet state vector (days) */
  const double AGEPMO=5.0;

  /*  Age limit for major-planet mean elements (days) */
  const double AGEPEL=50.0;

  /*  Margin for error when deciding whether to renew the planetary data */
  const double TINY=1e-6;

  /*  Age limit for the body's osculating elements (before rectification) */
  const double AGEBEL=100.0;

  /*  Gaussian gravitational constant squared */
  const double GCON2 = PAL__GCON * PAL__GCON;

  /*  The final epoch */
  double TFINAL;

  /*  The body's current universal elements */
  double UL[13];

  /*  Current reference epoch */
  double T0;

  /*  Timespan from latest orbit rectification to final epoch (days) */
  double TSPAN;

  /*  Time left to go before integration is complete */
  double TLEFT;

  /*  Time direction flag: +1=forwards, -1=backwards */
  double FB;

  /*  First-time flag */
  int FIRST = 0;

  /*
   *  The current perturbations
   */

  /*  Epoch (days relative to current reference epoch) */
  double RTN;
  /*  Position (AU) */
  double PERP[3];
  /*  Velocity (AU/d) */
  double PERV[3];
  /*  Acceleration (AU/d/d) */
  double PERA[3];

  /*  Length of current timestep (days), and half that */
  double TS,HTS;

  /*  Epoch of middle of timestep */
  double T;

  /*  Epoch of planetary mean elements */
  double TPEL = 0.0;

  /*  Planet number (1=Mercury, 2=Venus, 3=EMB...8=Neptune) */
  int NP;

  /*  Planetary universal orbital elements */
  double UP[8][13];

  /*  Epoch of planetary state vectors */
  double TPMO = 0.0;

  /*  State vectors for the major planets (AU,AU/s) */
  double PVIN[8][6];

  /*  Earth velocity and position vectors (AU,AU/s) */
  double VB[3],PB[3],VH[3],PE[3];

  /*  Moon geocentric state vector (AU,AU/s) and position part */
  double PVM[6],PM[3];

  /*  Date to J2000 de-precession matrix */
  double PMAT[3][3];

  /*
   *  Correction terms for extrapolated major planet vectors
   */

  /*  Sun-to-planet distances squared multiplied by 3 */
  double R2X3[8];
  /*  Sunward acceleration terms, G/2R^3 */
  double GC[8];
  /*  Tangential-to-circular correction factor */
  double FC;
  /*  Radial correction factor due to Sunwards acceleration */
  double FG;

  /*  The body's unperturbed and perturbed state vectors (AU,AU/s) */
  double PV0[6],PV[6];

  /*  The body's perturbed and unperturbed heliocentric distances (AU) cubed */
  double R03,R3;

  /*  The perturbating accelerations, indirect and direct */
  double FI[3],FD[3];

  /*  Sun-to-planet vector, and distance cubed */
  double RHO[3],RHO3;

  /*  Body-to-planet vector, and distance cubed */
  double DELTA[3],DELTA3;

  /*  Miscellaneous */
  int I,J;
  double R2,W,DT,DT2,R,FT;
  int NE;

  /*  Planetary inverse masses, Mercury through Neptune then Earth and Moon */
  const double AMAS[10] = {
    6023600., 408523.5, 328900.5, 3098710.,
    1047.355, 3498.5, 22869., 19314.,
    332946.038, 27068709.
  };

  /*  Preset the status to OK. */
  *jstat = 0;

  /*  Copy the final epoch. */
  TFINAL = date;

  /*  Copy the elements (which will be periodically updated). */
  for (I=0; I<13; I++) {
    UL[I] = u[I];
  }

/*  Initialize the working reference epoch. */
  T0=UL[2];

  /*  Total timespan (days) and hence time left. */
  TSPAN = TFINAL-T0;
  TLEFT = TSPAN;

  /*  Warn if excessive. */
  if (fabs(TSPAN) > 36525.0) *jstat=101;

  /*  Time direction: +1 for forwards, -1 for backwards. */
  FB = COPYSIGN(1.0,TSPAN);

  /*  Initialize relative epoch for start of current timestep. */
  RTN = 0.0;

  /*  Reset the perturbations (position, velocity, acceleration). */
  for (I=0; I<3; I++) {
    PERP[I] = 0.0;
    PERV[I] = 0.0;
    PERA[I] = 0.0;
  }

  /*  Set "first iteration" flag. */
  FIRST = 1;

  /*  Step through the time left. */
  while (FB*TLEFT > 0.0) {

    /*     Magnitude of current acceleration due to planetary attractions. */
    if (FIRST) {
      TS = TSMIN;
    } else {
      R2 = 0.0;
      for (I=0; I<3; I++) {
        W = FD[I];
        R2 = R2+W*W;
      }
      W = sqrt(R2);

      /*        Use the acceleration to decide how big a timestep can be tolerated. */
      if (W != 0.0) {
        TS = DMIN(TSMAX,DMAX(TSMIN,TSC/W));
      } else {
        TS = TSMAX;
      }
    }
    TS = TS*FB;

    /*     Override if final epoch is imminent. */
    TLEFT = TSPAN-RTN;
    if (fabs(TS) > fabs(TLEFT)) TS=TLEFT;

    /*     Epoch of middle of timestep. */
    HTS = TS/2.0;
    T = T0+RTN+HTS;

    /*     Is it time to recompute the major-planet elements? */
    if (FIRST || fabs(T-TPEL)-AGEPEL >= TINY) {

      /*        Yes: go forward in time by just under the maximum allowed. */
      TPEL = T+FB*AGEPEL;

      /*        Compute the state vector for the new epoch. */
      for (NP=1; NP<=8; NP++) {
        palPlanet(TPEL,NP,PV,&J);

        /*           Warning if remote epoch, abort if error. */
        if (J == 1) {
          *jstat = 102;
        } else if (J != 0) {
          goto ABORT;
        }

        /*           Transform the vector into universal elements. */
        palPv2ue(PV,TPEL,0.0,&(UP[NP-1][0]),&J);
        if (J != 0) goto ABORT;
      }
    }

    /*     Is it time to recompute the major-planet motions? */
    if (FIRST || fabs(T-TPMO)-AGEPMO >= TINY) {

      /*        Yes: look ahead. */
      TPMO = T+FB*AGEPMO;

      /*        Compute the motions of each planet (AU,AU/d). */
      for (NP=1; NP<=8; NP++) {

        /*           The planet's position and velocity (AU,AU/s). */
        palUe2pv(TPMO,&(UP[NP-1][0]),&(PVIN[NP-1][0]),&J);
        if (J != 0) goto ABORT;

        /*           Scale velocity to AU/d. */
        for (J=3; J<6; J++) {
          PVIN[NP-1][J] = PVIN[NP-1][J]*PAL__SPD;
        }

        /*           Precompute also the extrapolation correction terms. */
        R2 = 0.0;
        for (I=0; I<3; I++) {
          W = PVIN[NP-1][I];
          R2 = R2+W*W;
        }
        R2X3[NP-1] = R2*3.0;
        GC[NP-1] = GCON2/(2.0*R2*sqrt(R2));
      }
    }

    /*     Reset the first-time flag. */
    FIRST = 0;

    /*     Unperturbed motion of the body at middle of timestep (AU,AU/s). */
    palUe2pv(T,UL,PV0,&J);
    if (J != 0) goto ABORT;

    /*     Perturbed position of the body (AU) and heliocentric distance cubed. */
    R2 = 0.0;
    for (I=0; I<3; I++) {
      W = PV0[I]+PERP[I]+(PERV[I]+PERA[I]*HTS/2.0)*HTS;
      PV[I] = W;
      R2 = R2+W*W;
    }
    R3 = R2*sqrt(R2);

    /*     The body's unperturbed heliocentric distance cubed. */
    R2 = 0.0;
    for (I=0; I<3; I++) {
      W = PV0[I];
      R2 = R2+W*W;
    }
    R03 = R2*sqrt(R2);

    /*     Compute indirect and initialize direct parts of the perturbation. */
    for (I=0; I<3; I++) {
      FI[I] = PV0[I]/R03-PV[I]/R3;
      FD[I] = 0.0;
    }

    /*     Ready to compute the direct planetary effects. */

    /*     Reset the "near-Earth" flag. */
    NE = 0;

    /*     Interval from state-vector epoch to middle of current timestep. */
    DT = T-TPMO;
    DT2 = DT*DT;

    /*     Planet by planet, including separate Earth and Moon. */
    for (NP=1; NP<10; NP++) {

      /*        Which perturbing body? */
      if (NP <= 8) {

        /*           Planet: compute the extrapolation in longitude (squared). */
        R2 = 0.0;
        for (J=3; J<6; J++) {
          W = PVIN[NP-1][J]*DT;
          R2 = R2+W*W;
        }

        /*           Hence the tangential-to-circular correction factor. */
        FC = 1.0+R2/R2X3[NP-1];

        /*           The radial correction factor due to the inwards acceleration. */
        FG = 1.0-GC[NP-1]*DT2;

        /*           Planet's position. */
        for (I=0; I<3; I++) {
          RHO[I] = FG*(PVIN[NP-1][I]+FC*PVIN[NP-1][I+3]*DT);
        }

      } else if (NE) {

        /*           Near-Earth and either Earth or Moon. */

        if (NP == 9) {

          /*              Earth: position. */
          palEpv(T,PE,VH,PB,VB);
          for (I=0; I<3; I++) {
            RHO[I] = PE[I];
          }

        } else {

          /*              Moon: position. */
          palPrec(palEpj(T),2000.0,PMAT);
          palDmoon(T,PVM);
          eraRxp(PMAT,PVM,PM);
          for (I=0; I<3; I++) {
            RHO[I] = PM[I]+PE[I];
          }
        }
      }

      /*        Proceed unless Earth or Moon and not the near-Earth case. */
      if (NP <= 8 || NE) {

        /*           Heliocentric distance cubed. */
        R2 = 0.0;
        for (I=0; I<3; I++) {
          W = RHO[I];
          R2 = R2+W*W;
        }
        R = sqrt(R2);
        RHO3 = R2*R;

        /*           Body-to-planet vector, and distance. */
        R2 = 0.0;
        for (I=0; I<3; I++) {
          W = RHO[I]-PV[I];
          DELTA[I] = W;
          R2 = R2+W*W;
        }
        R = sqrt(R2);

        /*           If this is the EMB, set the near-Earth flag appropriately. */
        if (NP == 3 && R < RNE) NE = 1;

        /*           Proceed unless EMB and this is the near-Earth case. */
        if ( ! (NE && NP == 3) ) {

          /*              If too close, ignore this planet and set a warning. */
          if (R < COINC) {
            *jstat = NP;

          } else {

            /*                 Accumulate "direct" part of perturbation acceleration. */
            DELTA3 = R2*R;
            W = AMAS[NP-1];
            for (I=0; I<3; I++) {
              FD[I] = FD[I]+(DELTA[I]/DELTA3-RHO[I]/RHO3)/W;
            }
          }
        }
      }
    }

    /*     Update the perturbations to the end of the timestep. */
    RTN += TS;
    for (I=0; I<3; I++) {
      W = (FI[I]+FD[I])*GCON2;
      FT = W*TS;
      PERP[I] = PERP[I]+(PERV[I]+FT/2.0)*TS;
      PERV[I] = PERV[I]+FT;
      PERA[I] = W;
    }

    /*     Time still to go. */
    TLEFT = TSPAN-RTN;

    /*     Is it either time to rectify the orbit or the last time through? */
    if (fabs(RTN) >= AGEBEL || FB*TLEFT <= 0.0) {

      /*        Yes: update to the end of the current timestep. */
      T0 += RTN;
      RTN = 0.0;

      /*        The body's unperturbed motion (AU,AU/s). */
      palUe2pv(T0,UL,PV0,&J);
      if (J != 0) goto ABORT;

      /*        Add and re-initialize the perturbations. */
      for (I=0; I<3; I++) {
        J = I+3;
        PV[I] = PV0[I]+PERP[I];
        PV[J] = PV0[J]+PERV[I]/PAL__SPD;
        PERP[I] = 0.0;
        PERV[I] = 0.0;
        PERA[I] = FD[I]*GCON2;
      }

      /*        Use the position and velocity to set up new universal elements. */
      palPv2ue(PV,T0,0.0,UL,&J);
      if (J != 0) goto ABORT;

      /*        Adjust the timespan and time left. */
      TSPAN = TFINAL-T0;
      TLEFT = TSPAN;
    }

    /*     Next timestep. */
  }

  /*  Return the updated universal-element set. */
  for (I=0; I<13; I++) {
    u[I] = UL[I];
  }

  /*  Finished. */
  return;

  /*  Miscellaneous numerical error. */
 ABORT:
  *jstat = -1;
  return;
}