aboutsummaryrefslogtreecommitdiff
path: root/src/slalib/pv2el.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/slalib/pv2el.f')
-rw-r--r--src/slalib/pv2el.f351
1 files changed, 351 insertions, 0 deletions
diff --git a/src/slalib/pv2el.f b/src/slalib/pv2el.f
new file mode 100644
index 0000000..e45ad64
--- /dev/null
+++ b/src/slalib/pv2el.f
@@ -0,0 +1,351 @@
+ SUBROUTINE sla_PV2EL (PV, DATE, PMASS, JFORMR,
+ : JFORM, EPOCH, ORBINC, ANODE, PERIH,
+ : AORQ, E, AORL, DM, JSTAT)
+*+
+* - - - - - -
+* P V 2 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.
+*
+* Reference: Sterne, Theodore E., "An Introduction to Celestial
+* Mechanics", Interscience Publishers, 1960
+*
+* Called: sla_DRANRM
+*
+* P.T.Wallace Starlink 13 February 1999
+*
+* Copyright (C) 1999 Rutherford Appleton Laboratory
+*-
+
+ 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 sla_DRANRM
+
+
+* 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.AND.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)
+
+* 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 = sla_DRANRM(BIGOM)
+ E = ECC
+ IF (JF.EQ.1) THEN
+ PERIH = sla_DRANRM(PL)
+ AORL = sla_DRANRM(EL)
+ DM = DN
+ ELSE
+ PERIH = sla_DRANRM(OM)
+ IF (JF.EQ.2) AORL = sla_DRANRM(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
generated by cgit (git 2.34.1) at 2025-09-20 05:13:41 -0400