;+
; NAME:
;	KeplerOrbit
; PURPOSE:
;	Calculate positions and velocities for a simple Kepler orbit
;	(of body m2 relative to body m1)
; CATEGORY:
;	smei/gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION KeplerOrbit, T 		, $
		elements	= OrbitElements , $
		plane		= OrbitPlane	, $
		m1			= m1			, $
		m2			= m2			, $
		degrees 	= Degrees		, $
		velocity	= velocity		, $
		angularvelocity=angularvelocity
; INPUTS:
;	T				array; type: standard time structure
;						times (UT)
;	elements=OrbitElements
;					array[4]; type: float
;						OrbitElement[0]	Semi-major axis
;						OrbitElement[1]	Eccentricity
;						OrbitElement[2]	Argument of perihelion
;										=angle(ascending node-perihelion)
;						OrbitElement[3]	Time of perihelion passage in days since 2000, Jan 1.5 (njd)
; OPTIONAL INPUT PARAMETERS:
;	plane=OrbitPlane
;					array[2]; type: float
;						OrbitElement[0]	Longitude ascending node
;						OrbitElement[1]	Inclination of orbital plane
;	m1=m1			scalar; type: float: default: 1 (solar mass)
;						mass of primary body in units of the solar mass
;	m2=m2			scalar; type: float: default: 0 (ignored)
;						mass of secondary body in units of the solar mass
; OUTPUTS:
;	R				array[3,*]; type: float
;						longitude (deg/rad), latitude (deg/rad) and distance (AU)
;						the angular units depend on the setting of /degrees
; OPTIONAL OUTPUT PARAMETERS:
;	velocity=velocity
;					array[3,*]; type: float
;						velocity, tangential velocity and radial velocity (AU/day)
;						in orbital plane
;	angularvelocity=angularvelocity
;						longitude (deg/rad), latitude (deg/rad) and magnitude (AU/day)
;						of velocity vector
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; EXTERNAL:
;	EqKepler
; CALLS:
;	InitVar, destroyvar, ToRadians, TimeGet, nrZBrent, EulerRotate
; PROCEDURE:
;	Standard Kepler orbit stuff
; MODIFICATION HISTORY:
;	SEP-1998, Paul Hick (UCSD/CASS)
;	SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		modified to work with time structures
;-

rpm = ToRadians(degrees=Degrees)

aa	= OrbitElements[0]					; Semi-major axis
ee	= OrbitElements[1]					; Eccentricity
ww	= OrbitElements[2]*rpm				; Argument of perihelion=angle(ascending node-perihelion)
tt	= OrbitElements[3]					; Time of perihelion passage (njd)

InitVar, m1, 1.0d0						; Default: solar mass
InitVar, m2, 0.0d0						; Default: ignore m2

mu		= 0.1d0*!sun.gm*(m1+m2)			; units 10^27 cm^3/s^2
AU		=       !sun.au					; units 10^13 cm
SpD		= 0.086400d0					; 10^6 Seconds per day
TwoPi	= 2*!dpi
HalfPi	= !dpi/2
pm1		= 2*(ee LT 1)-1					; +1 for ellipse; -1 for hyperbole

nn		= sqrt(mu/(AU*aa)^3)*SpD		; Mean (angular) motion (radians/day)
							; M should be array[1] to avoid problems with transpose
M = nn*(TimeGet(T, /njd)-tt)			; Mean anomaly (radians)

M = M mod TwoPi
M = M+(M lt 0)*TwoPi					; In [0,2*pi]

E = 0.0d0*M								; Eccentric anomaly (radians)
FOR i=0,n_elements(M)-1 DO	$			; Solve Kepler's equation
	E[i] = nrZBrent('EqKepler',0.0d0,TwoPi,0.00001d0,arg=[M[i],ee])

IF ee LT 1 THEN	BEGIN					; Elliptical orbit
	r = 1.0d0-ee*cos(E)					; Distance to m1 in units of aa
	vv = atan(sqrt(1.0d0-ee*ee)*sin(E),cos(E)-ee) ; True anomaly (radians)
ENDIF ELSE BEGIN						; Hyperbolic orbit
	r = ee*cosh(E)-1.0d0
	vv = atan(sqrt(ee*ee-1.0d0)*sinh(E),ee-cosh(E)) ; True anomaly (radians)
ENDELSE

lr = vv+ww
lr = lr mod TwoPi
lr += (lr LT 0)*TwoPi

dr = 0.0d0*lr

v  = 2.0d0/r-pm1						; velocity^2
vt = pm1*(1.0d0-ee*ee)/(r*r)			; (tangential velocity)^2
vr = (v-vt) > 0							; (radial velocity)^2

v  = sqrt(v)
vt = sqrt(vt)
vr = sqrt(vr)*(2*(vv GT 0)-1)			; Sign depends on sign of vv

lv = lr+atan(vt,vr)
lv = lv mod TwoPi
lv += (lv LT 0)*TwoPi

dv = 0.0d0*lv
										; Convert from orbital plane to ecliptic
IF n_elements(OrbitPlane) NE 0 THEN BEGIN
	O  = HalfPi-OrbitPlane[0]*rpm
	ii = -OrbitPlane[1]*rpm

	g = EulerRotate([-HalfPi, ii, O], lr, dr)
	g = EulerRotate([-HalfPi, ii, O], lv, dv)
ENDIF

g  = (SpD/AU)*sqrt(mu/(AU*aa))			; Converts velocities to AU/day

orbitvector		=   [transpose(lr/rpm)	,transpose(dr/rpm)	,aa*transpose(r) ]
velocity 		= g*[transpose(v )		,transpose(vt)		,   transpose(vr)]
angularvelocity =   [transpose(lv/rpm)	,transpose(dv/rpm)	,g *transpose(v) ]

IF n_elements(T) EQ 1 THEN BEGIN
	orbitvector		= reform(orbitvector)
	velocity		= reform(velocity)
	angularvelocity = reform(angularvelocity)
ENDIF

RETURN, orbitvector  &  END