Permalink
Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
156 lines (138 sloc) 5.12 KB
SUBROUTINE sla_PLANTU (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: sla_GMST, sla_DT, sla_EPJ, sla_EPV, sla_UE2PV, sla_PRENUT,
* sla_DMXV, sla_PVOBS, sla_DCC2S, sla_DRANRM
*
* 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 sla_EPV (or sla_EVP) and
* sla_UE2PV.
*
* 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., 59 Temple Place, Suite 330,
* Boston, MA 02111-1307 USA
*
*-
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 sla_GMST,sla_DT,sla_EPJ,sla_DRANRM
* Sun to geocentre (J2000, velocity in AU/s).
CALL sla_EPV(DATE,VSG,VSG(4),DPB,DVB)
DO I=4,6
VSG(I)=VSG(I)/86400D0
END DO
* Sun to planet (J2000).
CALL sla_UE2PV(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 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