diff options
Diffstat (limited to 'src/slalib/rdplan.f')
| -rw-r--r-- | src/slalib/rdplan.f | 183 |
1 files changed, 183 insertions, 0 deletions
diff --git a/src/slalib/rdplan.f b/src/slalib/rdplan.f new file mode 100644 index 0000000..e00ec1d --- /dev/null +++ b/src/slalib/rdplan.f @@ -0,0 +1,183 @@ + SUBROUTINE sla_RDPLAN (DATE, NP, ELONG, PHI, RA, DEC, DIAM) +*+ +* - - - - - - - +* R D P L A N +* - - - - - - - +* +* 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 sla_DMOON and sla_PLANET. +* +* 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 sla_PLANET routine for references. +* +* Called: sla_GMST, sla_DT, sla_EPJ, sla_DMOON, sla_PVOBS, sla_PRENUT, +* sla_PLANET, sla_DMXV, sla_DCC2S, sla_DRANRM +* +* P.T.Wallace Starlink 26 May 1997 +* +* Copyright (C) 1997 Rutherford Appleton Laboratory +*- + + 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 sla_GMST,sla_DT,sla_EPJ,sla_DRANRM + +* 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=sla_GMST(DATE-sla_DT(sla_EPJ(DATE))/86400D0)+ELONG + +* Geocentre to Moon (mean of date) + CALL sla_DMOON(DATE,V) + +* Nutation to true of date + CALL sla_NUT(DATE,RMAT) + CALL sla_DMXV(RMAT,V,VGM) + CALL sla_DMXV(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 sla_PRENUT(2000D0,DATE,RMAT) + +* Sun to Earth-Moon Barycentre (J2000) + CALL sla_PLANET(DATE,3,V,J) + +* Precession and nutation to date + CALL sla_DMXV(RMAT,V,VSE) + CALL sla_DMXV(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 sla_PLANET(DATE,IP,V,J) + +* Precession and nutation to date + CALL sla_DMXV(RMAT,V,VSP) + CALL sla_DMXV(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 sla_PVOBS(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 sla_DCC2S(V,RA,DEC) + RA=sla_DRANRM(RA) + +* Angular diameter (radians) + DIAM=2D0*ASIN(EQRAU(IP)/(R*AUKM)) + + END |