Repatch (from linux) of OSX IRAF

1 files changed, 156 insertions, 0 deletions

diff --git a/math/slalib/plantu.f b/math/slalib/plantu.f
new file mode 100644
index 00000000..81e65148
--- /dev/null
+++ b/math/slalib/plantu.f
@@ -0,0 +1,156 @@
+      SUBROUTINE slPLTU (DATE, ELONG, PHI, U, RA, DEC, R, JSTAT)
+*+
+*     - - - - - - -
+*      P L A N T U
+*     - - - - - - -
+*
+*  Topocentric apparent RA,Dec of a Solar-System object whose
+*  heliocentric universal elements are known.
+*
+*  Given:
+*     DATE      d     TT MJD of observation (JD - 2400000.5)
+*     ELONG     d     observer's east longitude (radians)
+*     PHI       d     observer's geodetic latitude (radians)
+*     U       d(13)   universal elements
+*
+*               (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:
+*     RA,DEC    d     RA, Dec (topocentric apparent, radians)
+*     R         d     distance from observer (AU)
+*     JSTAT     i     status:  0 = OK
+*                             -1 = radius vector zero
+*                             -2 = failed to converge
+*
+*  Called: slGMST, slDT, slEPJ, slEPV, slUEPV, slPRNU,
+*          slDMXV, slPVOB, slDC2S, slDA2P
+*
+*  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 near-Earth asteroids.  Geocentric positions can be
+*     generated by appropriate use of routines slEPV (or slEVP) and
+*     slUEPV.
+*
+*  3  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.
+*
+*  4  The universal elements are with respect to the J2000 equator and
+*     equinox.
+*
+*     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.
+*
+*  Last revision:   19 February 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 DATE,ELONG,PHI,U(13),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 slGMST,slDT,slEPJ,slDA2P
+
+
+
+*  Sun to geocentre (J2000, velocity in AU/s).
+      CALL slEPV(DATE,VSG,VSG(4),DPB,DVB)
+      DO I=4,6
+         VSG(I)=VSG(I)/86400D0
+      END DO
+
+*  Sun to planet (J2000).
+      CALL slUEPV(DATE,U,VSP,JSTAT)
+
+*  Geocentre to planet (J2000).
+      DO I=1,6
+         V(I)=VSP(I)-VSG(I)
+      END DO
+
+*  Precession and nutation to date.
+      CALL slPRNU(2000D0,DATE,RMAT)
+      CALL slDMXV(RMAT,V,VGP)
+      CALL slDMXV(RMAT,V(4),VGP(4))
+
+*  Geocentre to observer (date).
+      STL=slGMST(DATE-slDT(slEPJ(DATE))/86400D0)+ELONG
+      CALL slPVOB(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 slDC2S(V,RA,DEC)
+      RA=slDA2P(RA)
+      R=D
+
+      END