PRO MOONPOS, dte, ra, dec, x, y, z
;+
; NAME:                                     
;	MOONPOS
; PURPOSE:
;	To compute the RA and Dec of the Moon at a given date.
;
; CALLING SEQUENCE:
;	MOONPOS, Dte, Ra, Dec, [X, Y, Z]
;
; INPUTS:
;	Dte - The Julian date of the day (and time) in question.
;
; OUTPUTS:
;	Ra  - The right ascension of the moon at that date in radians.
;	Dec - The declination of the moon at that date in radians.
;	X, Y, Z - the geocentric equatorial inertial coordinates of the
;	          center of the Moon w.r.t. the center of the Earth (Re).
;
; SIDE EFFECTS:
;	The user should be aware that the equations used are not highly 
;	accurate.    The error in right ascension will rarely exceed 0.3 
;	degrees; the error in declination will rarely exceed 0.2 degrees.
;
; PROCEDURE:
;     The lunar position is computed from the low precision formulae provided
;     in the Astronomical Almanac.
;
; EXAMPLE:
;	IDL> juldate,[1982,4,6],jd    ;Get reduced julian date
;	IDL> jd = jd + 2400000.       ;Convert to full julian date
;	IDL> moonpos, jd, ra ,dec     ;Get RA and Dec of moon
;	IDL> print,adstring(ra*!RADEG,dec*!RADEG)
;		==> 11 17  6.1  +09 17 56.0
;
;	This is about 5' from the position given in the Astronomical Almanac
; MODIFICATION HISTORY:
;	Written by Michael R. Greason, STX, 31 October 1988.
;       Modified by David P. Steele, CNSR, 17 February 1994:
;       Added code to calculate and return the GEI cartesian
;       coordinates of the Moon, using formulae from the Astronomical
;	Almanac.
;-
 On_error,2
 if N_params() EQ 0 then begin
     print,'Syntax - MOONPOS, dte, ra, dec'
     print,'dte = Julian Date, Output Ra and Dec are in radians'
     return
 endif
 t = (double(dte) - 2451545.0D0) / 36525.D0
 d2r = !DPI / 180.D0
 sz = size(dte)
 scalar = sz(0) eq 0               ;Scalar input?
 if scalar then t = dblarr(1) + t
;
 theta1 = (134.9D0 + (477198.85D0 * t)) * d2r
 lamda = 218.32D0 + (481267.883D0 * t) + (6.29D0 * sin(theta1))
 theta1 = (259.2D0 - (413335.38D0 * t)) * d2r
 theta2 = (235.7D0 + (890534.23D0 * t)) * d2r
 lamda = lamda - (1.27D0 * sin(theta1)) + (0.66D0 * sin(theta2))
 theta1 = (269.9D0 + (954397.70D0 * t)) * d2r
 theta2 = (357.5D0 + (35999.05D0 * t)) * d2r
 lamda = lamda + (0.21D0 * sin(theta1)) - (0.19D0 * sin(theta2))
 theta1 = (186.6D0 + (966404.05D0 * t)) * d2r
 lamda = lamda - (0.11D0 * sin(theta1))
 neg = where( lamda lt 0.D0, Nneg)
 if ( Nneg GT 0 ) then lamda(neg) = 360.D0 + ( lamda(neg) mod 360)
 lamda = (lamda mod 360.D0) * d2r
;
 theta1 = (93.3D0 + (483202.03D0 * t)) * d2r
 theta2 = (228.2D0 + (960400.87D0 * t)) * d2r
 beta = (5.13D0 * sin(theta1)) + (0.28D0 * sin(theta2))
 theta1 = (318.3D0 + (6003.18D0 * t)) * d2r
 theta2 = (217.6D0 - (407332.20D0 * t)) * d2r
 beta = beta - (0.28D0 * sin(theta1)) - (0.17D0 * sin(theta2))
 neg = where( beta LT 0.D0, Nneg )
 if ( Nneg GT 0 ) then beta(neg) = 360.D0 + ( beta(neg) mod 360)
 beta = (beta mod 360.D0) * d2r

 theta1 = (134.9 + 477198.85D0 * t)*d2r
 theta2 = (259.2 - 413335.38D0 * t)*d2r
 pi = (0.9508 + 0.0518*COS(theta1) + 0.0095*COS(theta2))
 theta1 = (235.7 + 890534.23D0 * t)*d2r
 theta2 = (269.9 + 954397.70D0 * t)*d2r
 pi = pi + 0.0078*COS(theta1) + 0.0028*COS(theta2)

 r = 1.0D0/SIN(pi*d2r)

 l = cos(beta) * cos(lamda)
 m = (0.9175D0 * cos(beta) * sin(lamda)) - (0.3978D0 * sin(beta))
 n = (0.3978D0 * cos(beta) * sin(lamda)) + (0.9175D0 * sin(beta))

 x = r*l
 y = r*m
 z = r*n
 
 ra = atan(m, l)
 neg = where( ra LT 0.0D, Nneg )
 if ( Nneg GT 0 ) then ra(neg) = 2.D0*!DPI + ra(neg)
 dec = asin(n)
 if scalar then begin
    ra = ra(0) & dec = dec(0) & x = x(0) & y = y(0) & z = z(0)
 endif

 return
 end