;+
; NAME:
;   HPRSTATN
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Compute high precision earth station positions in inertial coordinates
;
; MAJOR TOPICS:
;   Geometry
;
; CALLING SEQUENCE:
;   HPRSTATN, JDTT, R_ITRF, R_ECI, V_ECI, $
;             [ /JPL, /USE_EOP, /NO_UT1, TBASE= ]
;
; DESCRIPTION:
;
;  The procedure HPRSTATN computes the coordinates and velocities of
;  an earth station in J2000 equatorial earth-centered inertial
;  coordinates (ECI).  This may be useful in any application where an
;  earthbound observatory is used to collect data on a non-terrestrial
;  phenomenon.
;
;  The user must have the following routines involved: HPRNUTANG;
;  Markwardt Quaternion Library; JPLEPHREAD and JPLEPHINTERP (if JPL
;  keyword is used); EOPDATA (if USE_EOP keyword is set); and TAI_UTC.
;  Also, the appropriate data files for TAI-UTC and Earth Orientation
;  Parameters must be installed.
;
;  The user must specify the position of the earth station in
;  earth-centered, earth-fixed, cartesian coordinates of the ITRF.
;  The Z-axis points to terrestrial north, the X-axis lies in the
;  terrestrial equator pointing towards the Greenwich meridian, and
;  the Y-axis forms the right handed coordinate system.  Any
;  positional units may be specified.
;
;  For the highest precision, the preferred method is to know the
;  coordinates from a direct IRTF reduction, i.e., via VLBI.  If you
;  have geodetic longitude/latitude, a procedure is given below which
;  can covert geodetic coordinates referred to WGS84 to ITRF cartesian
;  coordinates.
;
;  The values returned are the earth-centered inertial J2000
;  coordinates and velocities of the station.  All the effects of
;  earth rotation, precession, nutation, and polar motion can be
;  included.  The user has a choice of the kinds of transformations
;  that are included (see JPL, USE_EOP and NO_UT1 keywords).
;
;  The returned positional units are the same as the input units.  The
;  returned velocity units are (units of input) PER SECOND.
;
;  It is possible specify more than one time, or more than one station
;  position, or both.  If both more than one time and position are
;  specified, then there must be an equal quantity of both.
;
;
; CARTESIAN COORDINATES FROM GEODETIC COORDINATES
;
;  For a station whose geodetic latitude LAT, longitude LON (where
;  positive is east) and elevation H specified relative to the
;  ellipsoid, the cartesian coordinates are:
;
;         R_ITRF = [ (A*C + H)*COS(LAT)*COS(LON), $
;                    (A*C + H)*COS(LAT)*SIN(LON), $
;                    (A*S + H)*SIN(LAT) ]
;
;  where for the WGS84 reference ellipsoid, the equatorial radius is
;  set to A = 6378.137 km, and the flattening factor F =
;  1/298.257223563, and
;
;      C = 1.0/SQRT(COS(LAT)^2 + (1 - F)^2*SIN(LAT)^2)
;      S = (1 - F)^2 * C
;
;  Reference: Explanatory Supplement to the Astronomical Almanac,
;  eqns. 3.351-1 and 3.352-2.
;
;
; INPUTS:
;
;   JDTT - a scalar or N-vector, the TT epoch time(s) for which
;          station coordinates are to be computed.
;
;          For reference, JDTT = JDTAI + 32.184/86400d, where JDTAI is
;          the international atomic time measured in days.  The value
;          of the keyword TBASE is added to JDTT to arrive at the
;          actual Julian date.
;
;   R_ITRF - cartesian coordinates of earth station.  Either a
;            3-vector, or a 3xN array.  Units can be any positional
;            units.
;
;
; OUTPUTS:
;
;   R_ECI - upon output, the coordinates of the station at the given
;           time(s), referred to the earth-centered J2000 coordinate
;           system.  Either a 3-vector or 3xN array depending on the
;           input.  Units are the same as for R_ITRF.
;
;   V_ECI - upon output, the velocities of the station at the given
;           time(s), referred to the earth-centered J2000 coordinate
;           system.  Either a 3-vector or 3xN array depending on the
;           input.  Units are (units of R_ITRF) PER SECOND.
;
;
; KEYWORD PARAMETERS:
;
;   TBASE - a fixed epoch time (Julian days) to be added to each value
;           of JDTT.  Since subtraction of large numbers occurs with
;           TBASE first, the greatest precision is achieved when TBASE
;           is expressed as a nearby julian epoch, JDTT is expressed
;           as a small offset from the fixed epoch.  
;           Default: 0
;
;   JPL - if set, then the JPL ephemeris is used to compute nutation
;         angles.  Otherwise the series representation of HPRNUTANG is
;         used.
;
;   USE_EOP - if set, then use earth orientation parameters, returned
;             by EOPDATA, to further refine the station coordinates.
;             Otherwise, only precession and nutation are used.
;
;   NO_UT1 - if set, then do not use the UT1-UTC conversion.
;   NO_PRECESSION - disable precession calculation.
;   NO_NUTATION - disable nutation calculation.
;   NO_POLAR_MOTION - disable polar motion calculation.
;
;
; EXAMPLE:
;
;   ;; ITRF coordinates of Deep Space Network Antenna 63 (METERS)
;   R_DSN63 = [+4849092.647d, -0360180.569d, +4115109.113d]
;
;   ;; Time: 2000/01/01 01:30 TT
;   JDTT = jday(2000d,1,1) + 1.5/24
;
;   ;; Compute position of antenna in J2000 coordinate system using
;   ;; full Earth Orientation Parameters.
;   HPRSTATN, JDTT, R_DSN63, R_ECI, V_ECI, /USE_EOP
;
;
; SEE ALSO:
;
;   HPRNUTANG, TAI_UTC (Markwardt Library)
;   EOPDATA, JPLEPHREAD, JPLEPHINTERP
;
;
; REFERENCES:
;
;   McCarthy, D. D. (ed.) 1996: IERS Conventions, IERS T.N. 21.
;     http://maia.usno.navy.mil/conventions.html
;
;   Seidelmann, P.K. 1992, *Explanatory Supplement to the Astronomical
;     Almanac*, ISBN 0-935702-68-7
;
;
; MODIFICATION HISTORY:
;   Written, 6 May 2002, CM
;   Documented, 12 May 2002, CM
;   Corrected discussion of geodetic coordinates, 26 May 2002, CM
;   Add NO_POLAR_MOTION keyword; only compute V_ECI if the variable
;    is passed, 07 Mar 2007, CM
;   Save some memory by deleting variables that are no longer used,
;    19 Dec 2008, CM
;   Some small documentation improvements, 16 Jan 2010, CM
;   Correct error in calculation of conversion from geodetic lat/lon
;     to cartesian, and add references, 2010-05-12, CM
;
;  $Id: hprstatn.pro,v 1.8 2010/05/12 22:17:32 craigm Exp $
;
;-
; Copyright (C) 2002, 2007, 2008, 2010, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

pro hprstatn, jdtt, r_iers, r_eci, v_eci, tbase=tbase0, $
              jpl=jpl, use_eop=use_eop, no_ut1=no_ut1, $
              no_precession=noprec, no_nutation=nonut, $
              no_polar_motion=nopm

  ;; Compute the earth-orientation angles
  hprnutang, jdtt, tbase=tbase0, zeta, theta, z, dpsi, deps, $
    mean_obliq=eps0, true_obliq=eps, $
    gas_time=gast, $
    polar_x=pmx, polar_y=pmy, use_eop=use_eop, no_ut1=no_ut1, $
    no_nutation=nonut, jpl=jpl

  ;; Standard value of rotational angular velocity of the earth (Aoki
  ;; et al 1982)
  omega = 7.2921151467d-5 ;; RADIAN/SEC

  if keyword_set(nopm) then begin
      pmx = 0
      pmy = 0
  endif

  ;; Rotation from terrestrial reference pole to the true celestial
  ;; pole of date, then application of earth spin using apparent
  ;; sidereal time.
  ;;              GAST   PMY   PMX
  qter = qteuler(['z',   'x',  'y'], $
                 +gast,  -pmy, -pmx)
  gast = 0 & pmx = 0 & pmy = 0

  ;; Position referred to true celestial pole of date
  r_tod = qtvrot(r_iers, temporary(qter))

  ;; Velocity referred to true celestial pole of date
  if arg_present(v_eci) then begin
      v_tod = r_tod*0
      v_tod(0,*) = -omega*r_tod(1,*)
      v_tod(1,*) = +omega*r_tod(0,*)
  endif

  if keyword_set(noprec) then begin
      zeta = 0d & theta = 0d & z = 0d 
  endif
  if keyword_set(nonut) then begin
      eps0 = 0d & dpsi = 0d & eps = 0d
  endif

  ;; Application of earth precession and nutation from true equator
  ;; and equinox of date, to mean equator and equinox of J2000
  ;; (celestial)
  ;;              Precession        Nutation 
  qcel = qteuler(['z','y','z',      'x','z','x'], $
                 -zeta, +theta, -z, +eps0, -dpsi, -eps)
  zeta = 0 & theta = 0 & z = 0 & eps0 = 0 & dpsi = 0 & eps = 0

  ;; Compute station position in Earth Centered Inertial coordinates
  DO_PRNUTROT:
  r_eci = qtvrot(temporary(r_tod), qcel)
  if arg_present(v_eci) then $
    v_eci = qtvrot(temporary(v_tod), temporary(qcel))

  return
end