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