function ephem, elem, t, uvw

;+
; NAME:
;  ephem
; PURPOSE: (one line)
;  Return heliocentric position and velocity
; DESCRIPTION:
;  Given t in Julian date and orbital elements, return heliocentric XYZ in AU, and, optionally, 
;  velocity in AU/day
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;  xyz = ephem(elem, t, uvw)
; INPUTS:
;  elem - vector of 6 or 7 elements
;      e      Eccentricity (degrees)
;      i      Inclination w.r.t xy-plane (degrees)
;     LAN     Longitude of Ascending Node (degrees)
;     APF     Argument of Perifocus (degrees)
;     tau     Julian Date of last periapsis
;      n      Mean motion (au/day)
;      A      semimajor axis.  If not passed, assume M = solar
;
;  t - Julian date
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  xyz - vector of position in AU
;  uvw - vector of velocity in AU/day
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2000, by Leslie Young, SwRI
;-

; some constants
day = 86400.                      ; day (in s)
au = 149597870.691d3              ; astronomical unit (in m)
deg = !pi/180.                    ; degrees (in radians)
GM = 1.3271243994d20              ; m^3/s^2.  JPL Horizons
GM = au^3 * (0.01720209895D/day)^2 ; GM for sun, 1.327124520514e+20 m^3/s^2
                                   ;  Fundemental definition, AA suplement p. 696

e  = elem[0] 
i  = elem[1]
LAN= elem[2]
APF= elem[3] 
tau= elem[4]
n  = elem[5]
if n_elements(elem) lt 7 then a = ( ( gm/(n*deg/day)^2)^(1/3.) ) / au else a = elem(6)
                                  ; GM for sun, 1.327124520514e+20 m^3/s^2
                                  ;  Fundemental definition, AA suplement p. 696

m = (n*deg) * (t-tau)             ; mean anomaly

; Eccentric anomlaly.  Solve Kepler's eqn for bound orbits, M = E - e sin E
ecc = m + e * sin(m)
decc = 1.D ; force at least one iteration
n_iter = 0
while (abs(decc) gt 1.e-15 and n_iter le 20) do begin &$
   decc = (m - ecc + e*sin(ecc))/(1 - e*cos(ecc)) &$
   ecc = ecc + decc &$
   n_iter = n_iter + 1 &$
end

; calculate true anomaly.  Use 2-argument form of tan to perserve quadrant
v = 2 * atan( sqrt(1+e) * sin(ecc/2.), sqrt(1-e) * cos(ecc/2.) )
;print, 'true anom ', v/deg

; calculate the distance
r = a * (1 - e * cos(ecc))

; calculate the position (Green p. 165)
cosLAN = cos(LAN*deg) & sinLAN = sin(LAN*deg) ; long of ascending node, Omega
cosi = cos(i*deg) & sini = sin(i*deg)         ; inclination
cosarg = cos(v+APF*deg) & sinarg=sin(v+APF*deg) ; true anom + argument of perifocus, omega

x = r * (cosarg*cosLAN - sinarg*sinLAN*cosi)
y = r * (cosarg*sinLAN + sinarg*cosLAN*cosi)
z = r * (sinarg*sini)

IF N_PARAMS(0) LT 3 THEN return, [x, y, z]  ; if velocity not needed, just return

; and now the velocity.  First calculate V0
;v0 = sqrt(GM/(a*au*(1-e^2)))
v0 = n*(deg/day)*a*au/sqrt(1-e^2)
vx = - v0 * sin(v)
vy = v0 * (cos(v) + e)
cosAPF = cos(APF*deg) & sinAPF=sin(APF*deg) ; argument of perifocus, omega

; make the three rotation matrices
r1 = [ [cosAPF, -sinAPF, 0], [sinAPF, cosAPF, 0], [ 0, 0, 1] ]
r2 = [ [  1   ,    0,    0], [ 0, cosi,   -sini], [ 0, sini, cosi] ]
r3 = [ [cosLAN, -sinLAN, 0], [sinLAN, cosLAN, 0], [ 0, 0, 1] ]
mat = r3 ## (r2 ## r1)
uvw = reform( mat ## reform([vx, vy, 0], 1,3) )/(au/day)
return, [x, y, z]

end