BTA_lib/slalib/plantu.f

      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