FUNCTION exofast_keplereq,m,ecc,thresh=thresh ;+ ; NAME: ; exofast_keplereq ; PURPOSE: ; Solve Kepler's Equation ; DESCRIPTION: ; Solve Kepler's Equation. Method by S. Mikkola (1987) Celestial ; Mechanics, 40 , 329-334. ; result from Mikkola then used as starting value for ; Newton-Raphson iteration to extend the applicability of this ; function to higher eccentricities ; ; CATEGORY: ; Celestial Mechanics ; CALLING SEQUENCE: ; eccanom=exofast_keplereq(m,ecc) ; INPUTS: ; m - Mean anomaly (radians; can be an array) ; ecc - Eccentricity ; OPTIONAL INPUT PARAMETERS: ; ; KEYWORD INPUT PARAMETERS: ; thresh: stopping criterion for the Newton Raphson iteration; the ; iteration stops once abs(E-Eold) pi zz=where(mx GT !dpi) IF zz[0] NE -1 THEN BEGIN mx[zz]=mx[zz] MOD (2*!dpi) zz=where(mx GT !dpi) IF zz[0] NE -1 THEN mx[zz]=mx[zz]-2.0D0*!dpi ENDIF ;; ... m < -pi zz=where(mx LE -!dpi) IF zz[0] NE -1 THEN BEGIN mx[zz]=mx[zz] MOD (2*!dpi) zz=where(mx LE -!dpi) IF zz[0] ne -1 THEN mx[zz]=mx[zz]+2.0D0*!dpi ENDIF ;; Bail out for circular orbits... IF (ecc EQ 0.) THEN return,mx ;; equation 9a aux = 4.d0*ecc+0.5d0 alpha = (1.d0-ecc)/aux beta=mx/(2.d0*aux) ;; equation 9b (except not really; is there an errata?) ;; the actual equation 9b is much much slower, but gives the same ;; answer (probably because more refinement necessary) aux=sqrt(beta*beta+alpha*alpha*alpha) z=beta+aux zz=where(z LE 0.0d0) if zz[0] ne -1 THEN z[zz]=beta[zz]-aux[zz] z = z^0.3333333333333333d0 s0=z-alpha/z s1=s0-(0.078d0*s0*s0*s0*s0*s0)/(1.d0+ecc) e0=mx+ecc*(3.d0*s1-4.d0*s1*s1*s1) se0=sin(e0) ce0=cos(e0) f = e0-ecc*se0-mx f1 = 1.d0-ecc*ce0 f2 = ecc*se0 f3 = ecc*ce0 u1 = -f/f1 u2 = -f/(f1+0.5d0*f2*u1) u3 = -f/(f1+0.5d0*f2*u2+0.16666666666666667d0*f3*u2*u2) u4 = -f/(f1+0.5d0*f2*u3+0.16666666666666667d0*f3*u3*u3-0.041666666666666667d0*f2*u3*u3*u3) eccanom=e0+u4 zz = where(eccanom GE 2.0d0*!dpi) IF zz[0] NE -1 THEN eccanom[zz]=eccanom[zz]-2.0d0*!dpi zz = where(eccanom LT 0.0d0) IF zz[0] NE -1 THEN eccanom[zz]=eccanom[zz]+2.0d0*!dpi ;; Now get more precise solution using Newton Raphson method ;; for those times when the Kepler equation is not yet solved ;; to better than 1e-10 ;; (modification J. Wilms) ndx=where(mx LT 0.) IF (ndx[0] NE -1 ) THEN mx[ndx]=mx[ndx]+2.*!dpi ;; calculate the differences diff = abs(eccanom-ecc*sin(eccanom) - mx) diff = diff < abs(diff - 2*!dpi) ;; 0 - epsilon = 2*pi - epsilon ndx = where(diff GT thresh) WHILE NDX[0] NE -1 DO BEGIN ;; E-e sinE-M fe=(eccanom[ndx]-ecc*sin(eccanom[ndx])-mx[ndx]) MOD (2*!dpi) ;; f' = 1-e*cosE fs=(1.d0-ecc*cos(eccanom[ndx])) MOD (2*!dpi) oldval=eccanom[ndx] eccanom[ndx]=(oldval-fe/fs) tmp = where(abs(oldval-eccanom[ndx]) GT thresh) if tmp[0] ne -1 then ndx = ndx[tmp] $ else break ENDWHILE ;; range reduction toohigh = where(eccanom GE 2*!dpi) if toohigh[0] ne -1 then eccanom[toohigh] = eccanom[toohigh] MOD (2*!dpi) toolow = where(eccanom LT 0) if toolow[0] ne -1 then eccanom[toolow] = (eccanom[toolow] MOD (2*!dpi)) + 2*!dpi return,eccanom END