Repatch (from linux) of OSX IRAF

1 files changed, 201 insertions, 0 deletions

diff --git a/math/slalib/rdplan.f b/math/slalib/rdplan.f
new file mode 100644
index 00000000..95ba3ef8
--- /dev/null
+++ b/math/slalib/rdplan.f
@@ -0,0 +1,201 @@
+      SUBROUTINE slRDPL (DATE, NP, ELONG, PHI, RA, DEC, DIAM)
+*+
+*     - - - - - - -
+*      R D P L
+*     - - - - - - -
+*
+*  Approximate topocentric apparent RA,Dec of a planet, and its
+*  angular diameter.
+*
+*  Given:
+*     DATE        d       MJD of observation (JD - 2400000.5)
+*     NP          i       planet: 1 = Mercury
+*                                 2 = Venus
+*                                 3 = Moon
+*                                 4 = Mars
+*                                 5 = Jupiter
+*                                 6 = Saturn
+*                                 7 = Uranus
+*                                 8 = Neptune
+*                                 9 = Pluto
+*                              else = Sun
+*     ELONG,PHI   d       observer's east longitude and geodetic
+*                                               latitude (radians)
+*
+*  Returned:
+*     RA,DEC      d        RA, Dec (topocentric apparent, radians)
+*     DIAM        d        angular diameter (equatorial, radians)
+*
+*  Notes:
+*
+*  1  The date is in a dynamical timescale (TDB, formerly ET) and is
+*     in the form of a Modified Julian Date (JD-2400000.5).  For all
+*     practical purposes, TT can be used instead of TDB, and for many
+*     applications UT will do (except for the Moon).
+*
+*  2  The longitude and latitude allow correction for geocentric
+*     parallax.  This is a major effect for the Moon, but in the
+*     context of the limited accuracy of the present routine its
+*     effect on planetary positions is small (negligible for the
+*     outer planets).  Geocentric positions can be generated by
+*     appropriate use of the routines slDMON and slPLNT.
+*
+*  3  The direction accuracy (arcsec, 1000-3000AD) is of order:
+*
+*            Sun              5
+*            Mercury          2
+*            Venus           10
+*            Moon            30
+*            Mars            50
+*            Jupiter         90
+*            Saturn          90
+*            Uranus          90
+*            Neptune         10
+*            Pluto            1   (1885-2099AD only)
+*
+*     The angular diameter accuracy is about 0.4% for the Moon,
+*     and 0.01% or better for the Sun and planets.
+*
+*  See the slPLNT routine for references.
+*
+*  Called: slGMST, slDT, slEPJ, slDMON, slPVOB, slPRNU,
+*          slPLNT, slDMXV, slDC2S, slDA2P
+*
+*  P.T.Wallace   Starlink   26 May 1997
+*
+*  Copyright (C) 1997 Rutherford Appleton Laboratory
+*
+*  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
+      INTEGER NP
+      DOUBLE PRECISION ELONG,PHI,RA,DEC,DIAM
+
+*  AU in km
+      DOUBLE PRECISION AUKM
+      PARAMETER (AUKM=1.49597870D8)
+
+*  Light time for unit distance (sec)
+      DOUBLE PRECISION TAU
+      PARAMETER (TAU=499.004782D0)
+
+      INTEGER IP,J,I
+      DOUBLE PRECISION EQRAU(0:9),STL,VGM(6),V(6),RMAT(3,3),
+     :                 VSE(6),VSG(6),VSP(6),VGO(6),DX,DY,DZ,R,TL
+      DOUBLE PRECISION slGMST,slDT,slEPJ,slDA2P
+
+*  Equatorial radii (km)
+      DATA EQRAU / 696000D0,2439.7D0,6051.9D0,1738D0,3397D0,71492D0,
+     :             60268D0,25559D0,24764D0,1151D0 /
+
+
+
+*  Classify NP
+      IP=NP
+      IF (IP.LT.0.OR.IP.GT.9) IP=0
+
+*  Approximate local ST
+      STL=slGMST(DATE-slDT(slEPJ(DATE))/86400D0)+ELONG
+
+*  Geocentre to Moon (mean of date)
+      CALL slDMON(DATE,V)
+
+*  Nutation to true of date
+      CALL slNUT(DATE,RMAT)
+      CALL slDMXV(RMAT,V,VGM)
+      CALL slDMXV(RMAT,V(4),VGM(4))
+
+*  Moon?
+      IF (IP.EQ.3) THEN
+
+*     Yes: geocentre to Moon (true of date)
+         DO I=1,6
+            V(I)=VGM(I)
+         END DO
+      ELSE
+
+*     No: precession/nutation matrix, J2000 to date
+         CALL slPRNU(2000D0,DATE,RMAT)
+
+*     Sun to Earth-Moon Barycentre (J2000)
+         CALL slPLNT(DATE,3,V,J)
+
+*     Precession and nutation to date
+         CALL slDMXV(RMAT,V,VSE)
+         CALL slDMXV(RMAT,V(4),VSE(4))
+
+*     Sun to geocentre (true of date)
+         DO I=1,6
+            VSG(I)=VSE(I)-0.012150581D0*VGM(I)
+         END DO
+
+*     Sun?
+         IF (IP.EQ.0) THEN
+
+*        Yes: geocentre to Sun
+            DO I=1,6
+               V(I)=-VSG(I)
+            END DO
+         ELSE
+
+*        No: Sun to Planet (J2000)
+            CALL slPLNT(DATE,IP,V,J)
+
+*        Precession and nutation to date
+            CALL slDMXV(RMAT,V,VSP)
+            CALL slDMXV(RMAT,V(4),VSP(4))
+
+*        Geocentre to planet
+            DO I=1,6
+               V(I)=VSP(I)-VSG(I)
+            END DO
+         END IF
+      END IF
+
+*  Refer to origin at the observer
+      CALL slPVOB(PHI,0D0,STL,VGO)
+      DO I=1,6
+         V(I)=V(I)-VGO(I)
+      END DO
+
+*  Geometric distance (AU)
+      DX=V(1)
+      DY=V(2)
+      DZ=V(3)
+      R=SQRT(DX*DX+DY*DY+DZ*DZ)
+
+*  Light time (sec)
+      TL=TAU*R
+
+*  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)
+
+*  Angular diameter (radians)
+      DIAM=2D0*ASIN(EQRAU(IP)/(R*AUKM))
+
+      END