path: root/src/slalib/pertue.f

SUBROUTINE sla_PERTUE (DATE, U, JSTAT)
*+
*     - - - - - - -
*      P E R T U E
*     - - - - - - -
*
*  Update the universal elements of an asteroid or comet by applying
*  planetary perturbations.
*
*  Given:
*     DATE     d       final epoch (TT MJD) for the updated elements
*
*  Given and returned:
*     U        d(13)   universal elements (updated in place)
*
*                (1)   combined mass (M+m)
*                (2)   total energy of the orbit (alpha)
*                (3)   reference (osculating) epoch (t0)
*              (4-6)   position at reference epoch (r0)
*              (7-9)   velocity at reference epoch (v0)
*               (10)   heliocentric distance at reference epoch
*               (11)   r0.v0
*               (12)   date (t)
*               (13)   universal eccentric anomaly (psi) of date, approx
*
*  Returned:
*     JSTAT    i       status:
*                          +102 = warning, distant epoch
*                          +101 = warning, large timespan ( > 100 years)
*                      +1 to +8 = coincident with major planet (Note 5)
*                             0 = OK
*                            -1 = numerical error
*
*  Called:  sla_PLANET, sla_UE2PV, sla_PV2UE
*
*  Notes:
*
*  1  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.
*
*  2  The universal elements are with respect to the J2000 equator and
*     equinox.
*
*  3  The epochs DATE, U(3) and U(12) are all Modified Julian Dates
*     (JD-2400000.5).
*
*  4  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 sla_PLANET.  Two
*     levels of interpolation are used, to enhance speed without
*     significantly degrading accuracy.  At a low frequency, the routine
*     sla_PLANET 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, relativistic effects
*     are not taken into account and the Earth-Moon system is treated as
*     a single body.
*
*     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.
*
*  5  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 (Mercury=1,...,Neptune=8) as a
*     warning.
*
*  References:
*
*     1  Sterne, Theodore E., "An Introduction to Celestial Mechanics",
*        Interscience Publishers Inc., 1960.  Section 6.7, p199.
*
*     2  Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
*
*  P.T.Wallace   Starlink   18 March 1999
*
*  Copyright (C) 1999 Rutherford Appleton Laboratory
*-

      IMPLICIT NONE
      DOUBLE PRECISION DATE,U(13)
      INTEGER JSTAT

*  Coefficient relating timestep to perturbing force
      DOUBLE PRECISION TSC
      PARAMETER (TSC=1D-4)

*  Minimum and maximum timestep (days)
      DOUBLE PRECISION TSMIN,TSMAX
      PARAMETER (TSMIN=0.01D0,TSMAX=10D0)

*  Age limit for major-planet state vector (days)
      DOUBLE PRECISION AGEPMO
      PARAMETER (AGEPMO=5D0)

*  Age limit for major-planet mean elements (days)
      DOUBLE PRECISION AGEPEL
      PARAMETER (AGEPEL=50D0)

*  Margin for error when deciding whether to renew the planetary data
      DOUBLE PRECISION TINY
      PARAMETER (TINY=1D-6)

*  Age limit for the body's osculating elements (before rectification)
      DOUBLE PRECISION AGEBEL
      PARAMETER (AGEBEL=100D0)

*  Gaussian gravitational constant (exact) and square
      DOUBLE PRECISION GCON,GCON2
      PARAMETER (GCON=0.01720209895D0,GCON2=GCON*GCON)

*  The final epoch
      DOUBLE PRECISION TFINAL

*  The body's current universal elements
      DOUBLE PRECISION UL(13)

*  Current reference epoch
      DOUBLE PRECISION T0

*  Timespan from latest orbit rectification to final epoch (days)
      DOUBLE PRECISION TSPAN

*  Time left to go before integration is complete
      DOUBLE PRECISION TLEFT

*  Time direction flag: +1=forwards, -1=backwards
      DOUBLE PRECISION FB

*  First-time flag
      LOGICAL FIRST

*
*  The current perturbations
*
*  Epoch (days relative to current reference epoch)
      DOUBLE PRECISION RTN
*  Position (AU)
      DOUBLE PRECISION PERP(3)
*  Velocity (AU/d)
      DOUBLE PRECISION PERV(3)
*  Acceleration (AU/d/d)
      DOUBLE PRECISION PERA(3)
*

*  Length of current timestep (days), and half that
      DOUBLE PRECISION TS,HTS

*  Epoch of middle of timestep
      DOUBLE PRECISION T

*  Epoch of planetary mean elements
      DOUBLE PRECISION TPEL

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

*  Planetary universal orbital elements
      DOUBLE PRECISION UP(13,8)

*  Epoch of planetary state vectors
      DOUBLE PRECISION TPMO

*  State vectors for the major planets (AU,AU/s)
      DOUBLE PRECISION PVIN(6,8)

*
*  Correction terms for extrapolated major planet vectors
*
*  Sun-to-planet distances squared multiplied by 3
      DOUBLE PRECISION R2X3(8)
*  Sunward acceleration terms, G/2R^3
      DOUBLE PRECISION GC(8)
*  Tangential-to-circular correction factor
      DOUBLE PRECISION FC
*  Radial correction factor due to Sunwards acceleration
      DOUBLE PRECISION FG
*

*  The body's unperturbed and perturbed state vectors (AU,AU/s)
      DOUBLE PRECISION PV0(6),PV(6)

*  The body's perturbed and unperturbed heliocentric distances (AU) cubed
      DOUBLE PRECISION R03,R3

*  The perturbating accelerations, indirect and direct
      DOUBLE PRECISION FI(3),FD(3)

*  Sun-to-planet vector, and distance cubed
      DOUBLE PRECISION RHO(3),RHO3

*  Body-to-planet vector, and distance cubed
      DOUBLE PRECISION DELTA(3),DELTA3

*  Miscellaneous
      INTEGER I,J
      DOUBLE PRECISION R2,W,DT,DT2,FT

*  Planetary inverse masses, Mercury through Neptune
      DOUBLE PRECISION AMAS(8)
      DATA AMAS / 6023600D0,408523.5D0,328900.5D0,3098710D0,
     :            1047.355D0,3498.5D0,22869D0,19314D0 /

*---------------------------------------------------------------------*


*  Preset the status to OK.
      JSTAT = 0

*  Copy the final epoch.
      TFINAL = DATE

*  Copy the elements (which will be periodically updated).
      DO I=1,13
         UL(I) = U(I)
      END DO

*  Initialize the working reference epoch.
      T0=UL(3)

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

*  Warn if excessive.
      IF (ABS(TSPAN).GT.36525D0) JSTAT=101

*  Time direction: +1 for forwards, -1 for backwards.
      FB = SIGN(1D0,TSPAN)

*  Initialize relative epoch for start of current timestep.
      RTN = 0D0

*  Reset the perturbations (position, velocity, acceleration).
      DO I=1,3
         PERP(I) = 0D0
         PERV(I) = 0D0
         PERA(I) = 0D0
      END DO

*  Set "first iteration" flag.
      FIRST = .TRUE.

*  Step through the time left.
      DO WHILE (FB*TLEFT.GT.0D0)

*     Magnitude of current acceleration due to planetary attractions.
         IF (FIRST) THEN
            TS = TSMIN
         ELSE
            R2 = 0D0
            DO I=1,3
               W = FD(I)
               R2 = R2+W*W
            END DO
            W = SQRT(R2)

*        Use the acceleration to decide how big a timestep can be tolerated.
            IF (W.NE.0D0) THEN
               TS = MIN(TSMAX,MAX(TSMIN,TSC/W))
            ELSE
               TS = TSMAX
            END IF
         END IF
         TS = TS*FB

*     Override if final epoch is imminent.
         TLEFT = TSPAN-RTN
         IF (ABS(TS).GT.ABS(TLEFT)) TS=TLEFT

*     Epoch of middle of timestep.
         HTS = TS/2D0
         T = T0+RTN+HTS

*     Is it time to recompute the major-planet elements?
         IF (FIRST.OR.ABS(T-TPEL)-AGEPEL.GE.TINY) THEN

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

*        Compute the state vector for the new epoch.
            DO NP=1,8
               CALL sla_PLANET(TPEL,NP,PV,J)

*           Warning if remote epoch, abort if error.
               IF (J.EQ.1) THEN
                  JSTAT = 102
               ELSE IF (J.NE.0) THEN
                  GO TO 9010
               END IF

*           Transform the vector into universal elements.
               CALL sla_PV2UE(PV,TPEL,0D0,UP(1,NP),J)
               IF (J.NE.0) GO TO 9010
            END DO
         END IF

*     Is it time to recompute the major-planet motions?
         IF (FIRST.OR.ABS(T-TPMO)-AGEPMO.GE.TINY) THEN

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

*        Compute the motions of each planet (AU,AU/d).
            DO NP=1,8

*           The planet's position and velocity (AU,AU/s).
               CALL sla_UE2PV(TPMO,UP(1,NP),PVIN(1,NP),J)
               IF (J.NE.0) GO TO 9010

*           Scale velocity to AU/d.
               DO J=4,6
                  PVIN(J,NP) = PVIN(J,NP)*86400D0
               END DO

*           Precompute also the extrapolation correction terms.
               R2 = 0D0
               DO I=1,3
                  W = PVIN(I,NP)
                  R2 = R2+W*W
               END DO
               R2X3(NP) = R2*3D0
               GC(NP) = GCON2/(2D0*R2*SQRT(R2))
            END DO
         END IF

*     Reset the first-time flag.
         FIRST = .FALSE.

*     Unperturbed motion of the body at middle of timestep (AU,AU/s).
         CALL sla_UE2PV(T,UL,PV0,J)
         IF (J.NE.0) GO TO 9010

*     Perturbed position of the body (AU) and heliocentric distance cubed.
         R2 = 0D0
         DO I=1,3
            W = PV0(I)+PERP(I)+(PERV(I)+PERA(I)*HTS/2D0)*HTS
            PV(I) = W
            R2 = R2+W*W
         END DO
         R3 = R2*SQRT(R2)

*     The body's unperturbed heliocentric distance cubed.
         R2 = 0D0
         DO I=1,3
            W = PV0(I)
            R2 = R2+W*W
         END DO
         R03 = R2*SQRT(R2)

*     Compute indirect and initialize direct parts of the perturbation.
         DO I=1,3
            FI(I) = PV0(I)/R03-PV(I)/R3
            FD(I) = 0D0
         END DO

*     Ready to compute the direct planetary effects.

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

*     Planet by planet.
         DO NP=1,8

*        First compute the extrapolation in longitude (squared).
            R2 = 0D0
            DO J=4,6
               W = PVIN(J,NP)*DT
               R2 = R2+W*W
            END DO

*        Hence the tangential-to-circular correction factor.
            FC = 1D0+R2/R2X3(NP)

*        The radial correction factor due to the inwards acceleration.
            FG = 1D0-GC(NP)*DT2

*        Planet's position, and heliocentric distance cubed.
            R2 = 0D0
            DO I=1,3
               W = FG*(PVIN(I,NP)+FC*PVIN(I+3,NP)*DT)
               RHO(I) = W
               R2 = R2+W*W
            END DO
            RHO3 = R2*SQRT(R2)

*        Body-to-planet vector, light-time, and distance cubed.
            R2 = 0D0
            DO I=1,3
               W = RHO(I)-PV(I)
               DELTA(I) = W
               R2 = R2+W*W
            END DO
            W = SQRT(R2)
            DELTA3 = R2*SQRT(R2)

*        If too close, ignore this planet and set a warning.
            IF (R2.LT.1D-6) THEN
               JSTAT = NP
            ELSE

*           Accumulate "direct" part of perturbation acceleration.
               W = AMAS(NP)
               DO I=1,3
                  FD(I) = FD(I)+(DELTA(I)/DELTA3-RHO(I)/RHO3)/W
               END DO
            END IF
         END DO

*     Update the perturbations to the end of the timestep.
         RTN = RTN+TS
         DO I=1,3
            W = (FI(I)+FD(I))*GCON2
            FT = W*TS
            PERP(I) = PERP(I)+(PERV(I)+FT/2D0)*TS
            PERV(I) = PERV(I)+FT
            PERA(I) = W
         END DO

*     Time still to go.
         TLEFT = TSPAN-RTN

*     Is it either time to rectify the orbit or the last time through?
         IF (ABS(RTN).GE.AGEBEL.OR.FB*TLEFT.LE.0D0) THEN

*        Yes: update to the end of the current timestep.
            T0 = T0+RTN
            RTN = 0D0

*        The body's unperturbed motion (AU,AU/s).
            CALL sla_UE2PV(T0,UL,PV0,J)
            IF (J.NE.0) GO TO 9010

*        Add and re-initialize the perturbations.
            DO I=1,3
               J = I+3
               PV(I) = PV0(I)+PERP(I)
               PV(J) = PV0(J)+PERV(I)/86400D0
               PERP(I) = 0D0
               PERV(I) = 0D0
               PERA(I) = FD(I)*GCON2
            END DO

*        Use the position and velocity to set up new universal elements.
            CALL sla_PV2UE(PV,T0,0D0,UL,J)
            IF (J.NE.0) GO TO 9010

*        Adjust the timespan and time left.
            TSPAN = TFINAL-T0
            TLEFT = TSPAN
         END IF

*     Next timestep.
      END DO

*  Return the updated universal-element set.
      DO I=1,13
         U(I) = UL(I)
      END DO

*  Finished.
      GO TO 9999

*  Miscellaneous numerical error.
 9010 CONTINUE
      JSTAT = -1

 9999 CONTINUE
      END