diff options
Diffstat (limited to 'math/slalib/rdplan.f')
| -rw-r--r-- | math/slalib/rdplan.f | 201 |
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 |