Initial commit

1 files changed, 157 insertions, 0 deletions

diff --git a/src/slalib/plante.f b/src/slalib/plante.f
new file mode 100644
index 0000000..e8f89d9
--- /dev/null
+++ b/src/slalib/plante.f
@@ -0,0 +1,157 @@
+      SUBROUTINE sla_PLANTE (DATE, ELONG, PHI, JFORM, EPOCH,
+     :                       ORBINC, ANODE, PERIH, AORQ, E,
+     :                       AORL, DM, RA, DEC, R, JSTAT)
+*+
+*     - - - - - - -
+*      P L A N T E
+*     - - - - - - -
+*
+*  Topocentric apparent RA,Dec of a Solar-System object whose
+*  heliocentric orbital elements are known.
+*
+*  Given:
+*     DATE      d     MJD of observation (JD - 2400000.5)
+*     ELONG     d     observer's east longitude (radians)
+*     PHI       d     observer's geodetic latitude (radians)
+*     JFORM     i     choice of element set (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 )
+*
+*  Returned:
+*     RA,DEC    d     RA, Dec (topocentric apparent, radians)
+*     R         d     distance from observer (AU)
+*     JSTAT     i     status:  0 = OK
+*                             -1 = illegal JFORM
+*                             -2 = illegal E
+*                             -3 = illegal AORQ
+*                             -4 = illegal DM
+*                             -5 = numerical error
+*
+*  Notes:
+*
+*  1  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).
+*
+*  2  The longitude and latitude allow correction for geocentric
+*     parallax.  This is usually a small effect, but can become
+*     important for Earth-crossing asteroids.  Geocentric positions
+*     can be generated by appropriate use of routines sla_EVP and
+*     sla_PLANEL.
+*
+*  3  The elements are with respect to the J2000 ecliptic and equinox.
+*
+*  4  Three different element-format options are available:
+*
+*     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
+*
+*  5  Unused elements (DM for JFORM=2, AORL and DM for JFORM=3) are
+*     not accessed.
+*
+*  Called: sla_GMST, sla_DT, sla_EPJ, sla_PVOBS, sla_PRENUT,
+*          sla_PLANEL, sla_DMXV, sla_DCC2S, sla_DRANRM
+*
+*  P.T.Wallace   Starlink   17 March 1999
+*
+*  Copyright (C) 1999 Rutherford Appleton Laboratory
+*-
+
+      IMPLICIT NONE
+
+      DOUBLE PRECISION DATE,ELONG,PHI
+      INTEGER JFORM
+      DOUBLE PRECISION EPOCH,ORBINC,ANODE,PERIH,AORQ,E,
+     :                 AORL,DM,RA,DEC,R
+      INTEGER JSTAT
+
+*  Light time for unit distance (sec)
+      DOUBLE PRECISION TAU
+      PARAMETER (TAU=499.004782D0)
+
+      INTEGER I
+      DOUBLE PRECISION DVB(3),DPB(3),VSG(6),VSP(6),V(6),RMAT(3,3),
+     :                 VGP(6),STL,VGO(6),DX,DY,DZ,D,TL
+      DOUBLE PRECISION sla_GMST,sla_DT,sla_EPJ,sla_DRANRM
+
+
+
+*  Sun to geocentre (J2000).
+      CALL sla_EVP(DATE,2000D0,DVB,DPB,VSG(4),VSG)
+
+*  Sun to planet (J2000).
+      CALL sla_PLANEL(DATE,JFORM,EPOCH,ORBINC,ANODE,PERIH,AORQ,
+     :                E,AORL,DM,VSP,JSTAT)
+
+*  Geocentre to planet (J2000).
+      DO I=1,6
+         V(I)=VSP(I)-VSG(I)
+      END DO
+
+*  Precession and nutation to date.
+      CALL sla_PRENUT(2000D0,DATE,RMAT)
+      CALL sla_DMXV(RMAT,V,VGP)
+      CALL sla_DMXV(RMAT,V(4),VGP(4))
+
+*  Geocentre to observer (date).
+      STL=sla_GMST(DATE-sla_DT(sla_EPJ(DATE))/86400D0)+ELONG
+      CALL sla_PVOBS(PHI,0D0,STL,VGO)
+
+*  Observer to planet (date).
+      DO I=1,6
+         V(I)=VGP(I)-VGO(I)
+      END DO
+
+*  Geometric distance (AU).
+      DX=V(1)
+      DY=V(2)
+      DZ=V(3)
+      D=SQRT(DX*DX+DY*DY+DZ*DZ)
+
+*  Light time (sec).
+      TL=TAU*D
+
+*  Correct position for planetary aberration
+      DO I=1,3
+         V(I)=V(I)-TL*V(I+3)
+      END DO
+
+*  To RA,Dec.
+      CALL sla_DCC2S(V,RA,DEC)
+      RA=sla_DRANRM(RA)
+      R=D
+
+      END