Tim:


Give the appended version of sla_UE2PV a go.  Changes:

* Newton's method for solving Everhart & Pitkin expression (17) is
  now damped 50% first time through.

* Number of iterations increased by 1 (for neatness as it was in fact
  stopping one iteration short of the stated number).

The code correctly implements the Everhart & Pitkin procedure, but there
is clearly some lingering scope for incipient instability.  The damping
fixes your example, but I wouldn't like to guarantee that the solution is
universal.


Patrick Wallace
____________________________________________________________________________
SSTD/HMNAO                                   Internet:  [log in to unmask]
Rutherford Appleton Laboratory                    Tel:       +44-1235-445372
Chilton, Didcot,                                  Fax:       +44-1235-446362
Oxon OX11 0QX, UK
____________________________________________________________________________

      SUBROUTINE sla_UE2PV (DATE, U, PV, JSTAT)
*+
*     - - - - - -
*      U E 2 P V
*     - - - - - -
*
*  Heliocentric position and velocity of a planet, asteroid or comet,
*  starting from orbital elements in the "universal variables" form.
*
*  Given:
*     DATE     d       date, Modified Julian Date (JD-2400000.5)
*
*  Given and returned:
*     U        d(13)   universal orbital elements (updated; Note 1)
*
*       given    (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
*     returned  (12)   date (t)
*         "     (13)   universal eccentric anomaly (psi) of date
*
*  Returned:
*     PV       d(6)    position (AU) and velocity (AU/s)
*     JSTAT    i       status:  0 = OK
*                              -1 = radius vector zero
*                              -2 = failed to converge
*
*  Notes
*
*  1  The "universal" elements are those which define the orbit for the
*     purposes of the method of universal variables (see reference).
*     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 companion routine is sla_EL2UE.  This takes the conventional
*     orbital elements and transforms them into the set of numbers
*     needed by the present routine.  A single prediction requires one
*     one call to sla_EL2UE followed by one call to the present routine;
*     for convenience, the two calls are packaged as the routine
*     sla_PLANEL.  Multiple predictions may be made by again
*     calling sla_EL2UE once, but then calling the present routine
*     multiple times, which is faster than multiple calls to sla_PLANEL.
*
*     It is not obligatory to use sla_EL2UE to obtain the parameters.
*     However, it should be noted that because sla_EL2UE performs its
*     own validation, no checks on the contents of the array U are made
*     by the present routine.
*
*  3  DATE is the instant for which the prediction is required.  It is
*     in the TT timescale (formerly Ephemeris Time, ET) and is a
*     Modified Julian Date (JD-2400000.5).
*
*  4  The universal elements supplied in the array U are in canonical
*     units (solar masses, AU and canonical days).  The position and
*     velocity are not sensitive to the choice of reference frame.  The
*     sla_EL2UE routine in fact produces coordinates with respect to the
*     J2000 equator and equinox.
*
*  5  The algorithm was originally adapted from the EPHSLA program of
*     D.H.P.Jones (private communication, 1996).  The method is based
*     on Stumpff's Universal Variables.
*
*  Reference:  Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
*
*  P.T.Wallace   Starlink   19 March 1999
*
*  Copyright (C) 1999 Rutherford Appleton Laboratory
*-

      IMPLICIT NONE

      DOUBLE PRECISION DATE,U(13),PV(6)
      INTEGER JSTAT

*  Gaussian gravitational constant (exact)
      DOUBLE PRECISION GCON
      PARAMETER (GCON=0.01720209895D0)

*  Canonical days to seconds
      DOUBLE PRECISION CD2S
      PARAMETER (CD2S=GCON/86400D0)

*  Test value for solution and maximum number of iterations
      DOUBLE PRECISION TEST
      INTEGER NITMAX
      PARAMETER (TEST=1D-13,NITMAX=20)

      INTEGER I,NIT,N

      DOUBLE PRECISION CM,ALPHA,T0,P0(3),V0(3),R0,SIGMA0,T,PSI,DT,W,
     :                 TOL,PSJ,PSJ2,BETA,S0,S1,S2,S3,FF,R,F,G,FD,GD



*  Unpack the parameters.
      CM = U(1)
      ALPHA = U(2)
      T0 = U(3)
      DO I=1,3
         P0(I) = U(I+3)
         V0(I) = U(I+6)
      END DO
      R0 = U(10)
      SIGMA0 = U(11)
      T = U(12)
      PSI = U(13)

*  Approximately update the universal eccentric anomaly.
      PSI = PSI+(DATE-T)*GCON/R0

*  Time from reference epoch to date (in Canonical Days: a canonical
*  day is 58.1324409... days, defined as 1/GCON).
      DT = (DATE-T0)*GCON

*  Refine the universal eccentric anomaly.
      NIT = 1
      W = 1D0
      TOL = 0D0
      DO WHILE (ABS(W).GE.TOL)

*     Form half angles until BETA small enough.
         N = 0
         PSJ = PSI
         PSJ2 = PSJ*PSJ
         BETA = ALPHA*PSJ2
         DO WHILE (ABS(BETA).GT.0.7D0)
            N = N+1
            BETA = BETA/4D0
            PSJ = PSJ/2D0
            PSJ2 = PSJ2/4D0
         END DO

*     Calculate Universal Variables S0,S1,S2,S3 by nested series.
         S3 = PSJ*PSJ2*((((((BETA/210D0+1D0)
     :                      *BETA/156D0+1D0)
     :                      *BETA/110D0+1D0)
     :                      *BETA/72D0+1D0)
     :                      *BETA/42D0+1D0)
     :                      *BETA/20D0+1D0)/6D0
         S2 = PSJ2*((((((BETA/182D0+1D0)
     :                  *BETA/132D0+1D0)
     :                  *BETA/90D0+1D0)
     :                  *BETA/56D0+1D0)
     :                  *BETA/30D0+1D0)
     :                  *BETA/12D0+1D0)/2D0
         S1 = PSJ+ALPHA*S3
         S0 = 1D0+ALPHA*S2

*     Undo the angle-halving.
         TOL = TEST
         DO WHILE (N.GT.0)
            S3 = 2D0*(S0*S3+PSJ*S2)
            S2 = 2D0*S1*S1
            S1 = 2D0*S0*S1
            S0 = 2D0*S0*S0-1D0
            PSJ = PSJ+PSJ
            TOL = TOL+TOL
            N = N-1
         END DO

*     Improve the approximation to PSI (damped 1st time through).
         FF = R0*S1+SIGMA0*S2+CM*S3-DT
         R = R0*S0+SIGMA0*S1+CM*S2
         IF (R.EQ.0D0) GO TO 9010
         W = FF/R
         IF (NIT.LT.2) W=W/2d0
         PSI = PSI-W

*     Next iteration, unless too many already.
         IF (NIT.GT.NITMAX) GO TO 9020
         NIT = NIT+1
      END DO

*  Project the position and velocity vectors (scaling velocity to AU/s).
      W = CM*S2
      F = 1D0-W/R0
      G = DT-CM*S3
      FD = -CM*S1/(R0*R)
      GD = 1D0-W/R
      DO I=1,3
         PV(I) = P0(I)*F+V0(I)*G
         PV(I+3) = CD2S*(P0(I)*FD+V0(I)*GD)
      END DO

*  Update the parameters to allow speedy prediction of PSI next time.
      U(12) = DATE
      U(13) = PSI

*  OK exit.
      JSTAT = 0
      GO TO 9999

*  Null radius vector.
 9010 CONTINUE
      JSTAT = -1
      GO TO 9999

*  Failed to converge.
 9020 CONTINUE
      JSTAT = -2

 9999 CONTINUE
      END