path: root/math/slalib/pv2el.f

SUBROUTINE slPVEL (PV, DATE, PMASS, JFORMR,
     :                      JFORM, EPOCH, ORBINC, ANODE, PERIH,
     :                      AORQ, E, AORL, DM, JSTAT)
*+
*     - - - - - -
*      P V E L
*     - - - - - -
*
*  Heliocentric osculating elements obtained from instantaneous position
*  and velocity.
*
*  Given:
*     PV        d(6)   heliocentric x,y,z,xdot,ydot,zdot of date,
*                      J2000 equatorial triad (AU,AU/s; Note 1)
*     DATE      d      date (TT Modified Julian Date = JD-2400000.5)
*     PMASS     d      mass of the planet (Sun=1; Note 2)
*     JFORMR    i      requested element set (1-3; Note 3)
*
*  Returned:
*     JFORM     d      element set actually returned (1-3; Note 4)
*     EPOCH     d      epoch of elements (TT MJD)
*     ORBINC    d      inclination (radians)
*     ANODE     d      longitude of the ascending node (radians)
*     PERIH     d      longitude or argument of perihelion (radians)
*     AORQ      d      mean distance or perihelion distance (AU)
*     E         d      eccentricity
*     AORL      d      mean anomaly or longitude (radians, JFORM=1,2 only)
*     DM        d      daily motion (radians, JFORM=1 only)
*     JSTAT     i      status:  0 = OK
*                              -1 = illegal PMASS
*                              -2 = illegal JFORMR
*                              -3 = position/velocity out of range
*
*  Notes
*
*  1  The PV 6-vector is with respect to the mean equator and equinox of
*     epoch J2000.  The orbital elements produced are with respect to
*     the J2000 ecliptic and mean equinox.
*
*  2  The mass, PMASS, is important only for the larger planets.  For
*     most purposes (e.g. asteroids) use 0D0.  Values less than zero
*     are illegal.
*
*  3  Three different element-format options are supported:
*
*     Option JFORM=1, suitable for the major planets:
*
*     EPOCH  = epoch of elements (TT MJD)
*     ORBINC = inclination i (radians)
*     ANODE  = longitude of the ascending node, big omega (radians)
*     PERIH  = longitude of perihelion, curly pi (radians)
*     AORQ   = mean distance, a (AU)
*     E      = eccentricity, e
*     AORL   = mean longitude L (radians)
*     DM     = daily motion (radians)
*
*     Option JFORM=2, suitable for minor planets:
*
*     EPOCH  = epoch of elements (TT MJD)
*     ORBINC = inclination i (radians)
*     ANODE  = longitude of the ascending node, big omega (radians)
*     PERIH  = argument of perihelion, little omega (radians)
*     AORQ   = mean distance, a (AU)
*     E      = eccentricity, e
*     AORL   = mean anomaly M (radians)
*
*     Option JFORM=3, suitable for comets:
*
*     EPOCH  = epoch of perihelion (TT MJD)
*     ORBINC = inclination i (radians)
*     ANODE  = longitude of the ascending node, big omega (radians)
*     PERIH  = argument of perihelion, little omega (radians)
*     AORQ   = perihelion distance, q (AU)
*     E      = eccentricity, e
*
*  4  It may not be possible to generate elements in the form
*     requested through JFORMR.  The caller is notified of the form
*     of elements actually returned by means of the JFORM argument:
*
*      JFORMR   JFORM     meaning
*
*        1        1       OK - elements are in the requested format
*        1        2       never happens
*        1        3       orbit not elliptical
*
*        2        1       never happens
*        2        2       OK - elements are in the requested format
*        2        3       orbit not elliptical
*
*        3        1       never happens
*        3        2       never happens
*        3        3       OK - elements are in the requested format
*
*  5  The arguments returned for each value of JFORM (cf Note 5: JFORM
*     may not be the same as JFORMR) are as follows:
*
*         JFORM         1              2              3
*         EPOCH         t0             t0             T
*         ORBINC        i              i              i
*         ANODE         Omega          Omega          Omega
*         PERIH         curly pi       omega          omega
*         AORQ          a              a              q
*         E             e              e              e
*         AORL          L              M              -
*         DM            n              -              -
*
*     where:
*
*         t0           is the epoch of the elements (MJD, TT)
*         T              "    epoch of perihelion (MJD, TT)
*         i              "    inclination (radians)
*         Omega          "    longitude of the ascending node (radians)
*         curly pi       "    longitude of perihelion (radians)
*         omega          "    argument of perihelion (radians)
*         a              "    mean distance (AU)
*         q              "    perihelion distance (AU)
*         e              "    eccentricity
*         L              "    longitude (radians, 0-2pi)
*         M              "    mean anomaly (radians, 0-2pi)
*         n              "    daily motion (radians)
*         -             means no value is set
*
*  6  At very small inclinations, the longitude of the ascending node
*     ANODE becomes indeterminate and under some circumstances may be
*     set arbitrarily to zero.  Similarly, if the orbit is close to
*     circular, the true anomaly becomes indeterminate and under some
*     circumstances may be set arbitrarily to zero.  In such cases,
*     the other elements are automatically adjusted to compensate,
*     and so the elements remain a valid description of the orbit.
*
*  7  The osculating epoch for the returned elements is the argument
*     DATE.
*
*  Reference:  Sterne, Theodore E., "An Introduction to Celestial
*              Mechanics", Interscience Publishers, 1960
*
*  Called:  slDA2P
*
*  Last revision:   8 September 2005
*
*  Copyright P.T.Wallace.  All rights reserved.
*
*  License:
*    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 2 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 (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      DOUBLE PRECISION PV(6),DATE,PMASS
      INTEGER JFORMR,JFORM
      DOUBLE PRECISION EPOCH,ORBINC,ANODE,PERIH,AORQ,E,AORL,DM
      INTEGER JSTAT

*  Seconds to days
      DOUBLE PRECISION DAY
      PARAMETER (DAY=86400D0)

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

*  Sin and cos of J2000 mean obliquity (IAU 1976)
      DOUBLE PRECISION SE,CE
      PARAMETER (SE=0.3977771559319137D0,
     :           CE=0.9174820620691818D0)

*  Minimum allowed distance (AU) and speed (AU/day)
      DOUBLE PRECISION RMIN,VMIN
      PARAMETER (RMIN=1D-3,VMIN=1D-8)

*  How close to unity the eccentricity has to be to call it a parabola
      DOUBLE PRECISION PARAB
      PARAMETER (PARAB=1D-8)

      DOUBLE PRECISION X,Y,Z,XD,YD,ZD,R,V2,V,RDV,GMU,HX,HY,HZ,
     :                 HX2PY2,H2,H,OI,BIGOM,AR,ECC,S,C,AT,U,OM,
     :                 GAR3,EM1,EP1,HAT,SHAT,CHAT,AE,AM,DN,PL,
     :                 EL,Q,TP,THAT,THHF,F

      INTEGER JF

      DOUBLE PRECISION slDA2P


*  Validate arguments PMASS and JFORMR.
      IF (PMASS.LT.0D0) THEN
         JSTAT = -1
         GO TO 999
      END IF
      IF (JFORMR.LT.1.OR.JFORMR.GT.3) THEN
         JSTAT = -2
         GO TO 999
      END IF

*  Provisionally assume the elements will be in the chosen form.
      JF = JFORMR

*  Rotate the position from equatorial to ecliptic coordinates.
      X = PV(1)
      Y = PV(2)*CE+PV(3)*SE
      Z = -PV(2)*SE+PV(3)*CE

*  Rotate the velocity similarly, scaling to AU/day.
      XD = DAY*PV(4)
      YD = DAY*(PV(5)*CE+PV(6)*SE)
      ZD = DAY*(-PV(5)*SE+PV(6)*CE)

*  Distance and speed.
      R = SQRT(X*X+Y*Y+Z*Z)
      V2 = XD*XD+YD*YD+ZD*ZD
      V = SQRT(V2)

*  Reject unreasonably small values.
      IF (R.LT.RMIN.OR.V.LT.VMIN) THEN
         JSTAT = -3
         GO TO 999
      END IF

*  R dot V.
      RDV = X*XD+Y*YD+Z*ZD

*  Mu.
      GMU = (1D0+PMASS)*GCON*GCON

*  Vector angular momentum per unit reduced mass.
      HX = Y*ZD-Z*YD
      HY = Z*XD-X*ZD
      HZ = X*YD-Y*XD

*  Areal constant.
      HX2PY2 = HX*HX+HY*HY
      H2 = HX2PY2+HZ*HZ
      H = SQRT(H2)

*  Inclination.
      OI = ATAN2(SQRT(HX2PY2),HZ)

*  Longitude of ascending node.
      IF (HX.NE.0D0.OR.HY.NE.0D0) THEN
         BIGOM = ATAN2(HX,-HY)
      ELSE
         BIGOM=0D0
      END IF

*  Reciprocal of mean distance etc.
      AR = 2D0/R-V2/GMU

*  Eccentricity.
      ECC = SQRT(MAX(1D0-AR*H2/GMU,0D0))

*  True anomaly.
      S = H*RDV
      C = H2-R*GMU
      IF (S.NE.0D0.OR.C.NE.0D0) THEN
         AT = ATAN2(S,C)
      ELSE
         AT = 0D0
      END IF

*  Argument of the latitude.
      S = SIN(BIGOM)
      C = COS(BIGOM)
      U = ATAN2((-X*S+Y*C)*COS(OI)+Z*SIN(OI),X*C+Y*S)

*  Argument of perihelion.
      OM = U-AT

*  Capture near-parabolic cases.
      IF (ABS(ECC-1D0).LT.PARAB) ECC=1D0

*  Comply with JFORMR = 1 or 2 only if orbit is elliptical.
      IF (ECC.GE.1D0) JF=3

*  Functions.
      GAR3 = GMU*AR*AR*AR
      EM1 = ECC-1D0
      EP1 = ECC+1D0
      HAT = AT/2D0
      SHAT = SIN(HAT)
      CHAT = COS(HAT)

*  Variable initializations to avoid compiler warnings.
      AM = 0D0
      DN = 0D0
      PL = 0D0
      EL = 0D0
      Q = 0D0
      TP = 0D0

*  Ellipse?
      IF (ECC.LT.1D0 ) THEN

*     Eccentric anomaly.
         AE = 2D0*ATAN2(SQRT(-EM1)*SHAT,SQRT(EP1)*CHAT)

*     Mean anomaly.
         AM = AE-ECC*SIN(AE)

*     Daily motion.
         DN = SQRT(GAR3)
      END IF

*  "Major planet" element set?
      IF (JF.EQ.1) THEN

*     Longitude of perihelion.
         PL = BIGOM+OM

*     Longitude at epoch.
         EL = PL+AM
      END IF

*  "Comet" element set?
      IF (JF.EQ.3) THEN

*     Perihelion distance.
         Q = H2/(GMU*EP1)

*     Ellipse, parabola, hyperbola?
         IF (ECC.LT.1D0) THEN

*        Ellipse: epoch of perihelion.
            TP = DATE-AM/DN
         ELSE

*        Parabola or hyperbola: evaluate tan ( ( true anomaly ) / 2 )
            THAT = SHAT/CHAT
            IF (ECC.EQ.1D0) THEN

*           Parabola: epoch of perihelion.
               TP = DATE-THAT*(1D0+THAT*THAT/3D0)*H*H2/(2D0*GMU*GMU)
            ELSE

*           Hyperbola: epoch of perihelion.
               THHF = SQRT(EM1/EP1)*THAT
               F = LOG(1D0+THHF)-LOG(1D0-THHF)
               TP = DATE-(ECC*SINH(F)-F)/SQRT(-GAR3)
            END IF
         END IF
      END IF

*  Return the appropriate set of elements.
      JFORM = JF
      ORBINC = OI
      ANODE = slDA2P(BIGOM)
      E = ECC
      IF (JF.EQ.1) THEN
         PERIH = slDA2P(PL)
         AORL = slDA2P(EL)
         DM = DN
      ELSE
         PERIH = slDA2P(OM)
         IF (JF.EQ.2) AORL = slDA2P(AM)
      END IF
      IF (JF.NE.3) THEN
         EPOCH = DATE
         AORQ = 1D0/AR
      ELSE
         EPOCH = TP
         AORQ = Q
      END IF
      JSTAT = 0

 999  CONTINUE
      END