xyz.pro

basic median filter simulation

pro xyz,date,x,y,z,xvel,yvel,zvel,equinox=equinox

;+
; NAME:
;       XYZ
; PURPOSE:
;       Calculate geocentric X,Y, and Z  and velocity coordinates of the Sun
; EXPLANATION:
;       Calculates geocentric X,Y, and Z vectors and velocity coordinates 
;       (dx, dy and dz) of the Sun.   (The positive X axis is directed towards 
;       the equinox, the y-axis, towards the point on the equator at right 
;       ascension 6h, and the z axis toward the north pole of the equator).
;       Typical position accuracy is <1e-4 AU (15000 km).
;
; CALLING SEQUENCE:
;       XYZ, date, x, y, z, [ xvel, yvel, zvel, EQUINOX = ]
;
; INPUT:
;       date: reduced julian date (=JD - 2400000), scalar or vector
;
; OUTPUT:
;       x,y,z: scalars or vectors giving heliocentric rectangular coordinates
;                 (in A.U) for each date supplied.    Note that sqrt(x^2 + y^2
;                 + z^2) gives the Earth-Sun distance for the given date.
;       xvel, yvel, zvel: velocity vectors corresponding to X, Y and Z.
;
; OPTIONAL KEYWORD INPUT:
;       EQUINOX: equinox of output. Default is 1950.
;
; EXAMPLE:
;       What were the rectangular coordinates and velocities of the Sun on 
;       Jan 22, 1999 0h UT (= JD 2451200.5) in J2000 coords? NOTE:
;       Astronomical Almanac (AA) is in TDT, so add 64 seconds to 
;       UT to convert.
;
;       IDL> xyz,51200.5+64.d/86400.d,x,y,z,xv,yv,zv,equinox = 2000
;
;       Compare to Astronomical Almanac (1999 page C20)
;                   X  (AU)        Y  (AU)     Z (AU)
;       XYZ:      0.51456871   -0.76963263  -0.33376880
;       AA:       0.51453130   -0.7697110   -0.3337152
;       abs(err): 0.00003739    0.00007839   0.00005360
;       abs(err)
;           (km):   5609          11759         8040 
;
;       NOTE: Velocities in AA are for Earth/Moon barycenter
;             (a very minor offset) see AA 1999 page E3
;                  X VEL (AU/DAY) YVEL (AU/DAY)   Z VEL (AU/DAY)
;       XYZ:      -0.014947268   -0.0083148382    -0.0036068577
;       AA:       -0.01494574    -0.00831185      -0.00360365
;       abs(err):  0.000001583    0.0000029886     0.0000032077
;       abs(err)
;        (km/sec): 0.00265        0.00519          0.00557
;
; PROCEDURE CALLS:
;       PRECESS_XYZ
; REVISION HISTORY
;       Original algorithm from Almanac for Computers, Doggett et al. USNO 1978
;       Adapted from the book Astronomical Photometry by A. Henden
;       Written  W. Landsman   STX       June 1989
;       Correct error in X coefficient   W. Landsman HSTX  January 1995
;       Added velocities, more terms to positions and EQUINOX keyword,
;          some minor adjustments to calculations 
;          P. Plait/ACC March 24, 1999
;-

   On_error,2
  
   if (n_params() eq 0) then begin
      print,'Syntax - XYZ, date, x, y, z, [ xvel, yvel, zvel, EQUINOX= ]'
      print,'     (date is REDUCED Julian date (JD - 2400000.0) )'
      return
   endif

   picon = !DPI/180.0d
   t = (date - 15020.0d0)/36525.0d0         ;Relative Julian century from 1900

; NOTE: longitude arguments below are given in *equinox* of date.
;   Precess these to equinox 1950 to give everything an even footing.
;   Compute argument of precession from equinox of date back to 1950
   pp = (1.396041d + 0.000308d*(t + 0.5d))*(t-0.499998d)

; Compute mean solar longitude, precessed back to 1950
   el = 279.696678D + 36000.76892D*t + 0.000303d*t*t - pp

; Compute Mean longitude of the Moon
   c = 270.434164d + 480960.d*t + 307.883142d*t - 0.001133d*t*t - pp

; Compute longitude of Moon's ascending node
   n = 259.183275d - 1800.d*t - 134.142008d*t + 0.002078d*t*t - pp

; Compute mean solar anomaly
   g = 358.475833d + 35999.04975d*t - 0.00015d*t*t

; Compute the mean jupiter anomaly
   j = 225.444651d + 2880.0d*t + 154.906654d*t*t

; Compute mean anomaly of Venus
   v = 212.603219d + 58320.d*t + 197.803875d*t + 0.001286d*t*t

; Compute mean anomaly of Mars
   m = 319.529425d + 19080.d*t + 59.8585d*t + 0.000181d*t*t

; Convert degrees to radians for trig functions
   el = el*picon
   g  = g*picon
   j =  j*picon
   c  = c*picon
   v  = v*picon
   n  = n*picon
   m  = m*picon

; Calculate X,Y,Z using trigonometric series
   X =   0.999860d*cos(el)                          $
       - 0.025127d*cos(g - el)                      $
       + 0.008374d*cos(g + el)                      $
       + 0.000105d*cos(g + g + el)                  $
       + 0.000063d*t*cos(g - el)                    $
       + 0.000035d*cos(g + g - el)                  $
       - 0.000026d*sin(g - el - j)                  $
       - 0.000021d*t*cos(g + el)                    $
       + 0.000018d*sin(2.d*g + el - 2.d*v)          $
       + 0.000017d*cos(c)                           $
       - 0.000014d*cos(c - 2.d*el)                  $
       + 0.000012d*cos(4.d*g + el - 8.d*m + 3.d*j)  $
       - 0.000012d*cos(4.d*g - el - 8.d*m + 3.d*j)  $
       - 0.000012d*cos(g + el - v)                  $
       + 0.000011d*cos(2.d*g + el - 2.d*v)          $
       + 0.000011d*cos(2.d*g - el - 2.d*j)         
  

   Y =   0.917308d*sin(el)                             $
       + 0.023053d*sin(g - el)                         $
       + 0.007683d*sin(g + el)                         $
       + 0.000097d*sin(g + g + el)                     $
       - 0.000057d*t*sin(g - el)                       $
       - 0.000032d*sin(g + g - el)                     $
       - 0.000024d*cos(g - el - j)                     $
       - 0.000019d*t*sin(g + el)                       $
       - 0.000017d*cos(2.d0*g + el - 2.d0*v)           $
       + 0.000016d*sin(c)                              $
       + 0.000013d*sin(c - 2.d0*el )                   $
       + 0.000011d*sin(4.d0*g + el - 8.d0*m + 3.d0*j)  $
       + 0.000011d*sin(4.d0*g - el - 8.d0*m + 3.d0*j)  $
       - 0.000011d*sin(g + el - v)                     $
       + 0.000010d*sin(2.d0*g + el - 2.d0*v )          $
       - 0.000010d*sin(2.d0*g - el - 2.d0*j )         


   Z =   0.397825d*sin(el)        $
       + 0.009998d*sin(g-el)      $
       + 0.003332d*sin(g+el)      $
       + 0.000042d*sin(g+g+el)    $
       - 0.000025d*t*sin(g-el)    $
       - 0.000014d*sin(g+g-el)    $
       - 0.000010d*cos(g-el-j)    

;Precess_to new equator?
   if keyword_set(equinox) then precess_xyz, x, y, z, 1950, equinox

   if N_params() LE 3 then return
   
   XVEL = -0.017200d * sin(el)           $
          -0.000288d * sin(g + el)       $
          -0.000005d * sin(2.d0*g + el)  $
          -0.000004d * sin(c)            $
          +0.000003d * sin(c - 2.d0*el)  $
          +0.000001d *t * sin(g+el)      $
          -0.000001d * sin(2.d0*g-el)           
 
   YVEL =  0.015780 * cos(el)            $
          +0.000264 * cos(g + el)        $
          +0.000005 * cos(2.d0*g + el)   $
          +0.000004 * cos(c)             $
          +0.000003 * cos(c - 2.d0*el)   $
          -0.000001 * t * cos(g + el)    

   ZVEL = 0.006843 * cos(el)             $
         +0.000115 * cos(g  + el)        $
         +0.000002 * cos(2.d0*g + el)    $
         +0.000002 * cos(c)              $
         +0.000001 * cos(c - 2.d0*el)    

;Precess to new equator?

   if keyword_set(equinox) then precess_xyz, xvel, yvel, zvel, 1950, equinox

   return
   end