FUNCTION KeplerOrbit, T         , $
        elements    = OrbitElements , $
        plane       = OrbitPlane    , $
        m1          = m1            , $
        m2          = m2            , $
        degrees     = Degrees       , $
        velocity    = velocity      , $
        angularvelocity=angularvelocity
    @compile_opt.pro        ; On error, return to caller

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