-
Notifications
You must be signed in to change notification settings - Fork 15
/
Copy pathexofast_keplereq.pro
executable file
·166 lines (146 loc) · 4.86 KB
/
exofast_keplereq.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
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)<thresh
; OUTPUTS:
; the function returns the eccentric anomaly
; KEYWORD OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
; 2002/05/29 - Marc W. Buie, Lowell Observatory. Ported from fortran routines
; supplied by Larry Wasserman and Ted Bowell.
; http://www.lowell.edu/users/buie/
;
; 2002-09-09 -- Joern Wilms, IAA Tuebingen, Astronomie.
; use analytical values obtained for the low eccentricity case as
; starting values for a Newton-Raphson method to allow high
; eccentricity values as well
;
; $Log: keplereq.pro,v $
; Revision 1.4 2009/11/04 Jason Eastman (OSU)
; Changed name from keplereq to exofast_keplereq.pro
; Fix infinite loop when diff is near 2pi +/- epsilon.
; Fix infinite loop when oldval and eccanom[i] differ by integer
; multiples of 2pi.
; Consistent implementation of THRESH
; Made 1/6 and 1/24 take full advantage of IDL double precision
; Vectorized for loop -- factors of several faster for (rare) cases
; with many elements to be refined
;
; Revision 1.3 2005/05/25 16:11:35 wilms
; speed up: Newton Raphson is only done if necessary
; (i.e., almost never)
;
; Revision 1.2 2004/08/05 10:02:05 wilms
; now also works for more than 32000 time values
;
; Revision 1.1 2002/09/09 14:54:11 wilms
; initial release into aitlib
;-
;; set default values
IF n_elements(thresh) EQ 0 THEN thresh=1d-10
IF (ecc LT 0. OR ecc GE 1.) THEN BEGIN
message,'Eccentricity must be 0<= ecc. < 1'
ENDIF
;; so we don't modify the inputs
mx=m
;; Range reduction of m to -pi < m <= pi
;; ... m > 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
; eccanom = exofast_mod(eccanom,2d0*!dpi)
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