forked from shergreen/SatGen
-
Notifications
You must be signed in to change notification settings - Fork 0
/
orbit.py
208 lines (165 loc) · 7.07 KB
/
orbit.py
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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
############################# orbit class ###############################
# Arthur Fangzhou Jiang 2016, HUJI --- original version
# Arthur Fangzhou Jiang 2019, HUJI, UCSC --- revisions:
# - no longer convert speed unit from kpc/Gyr to km/s
# - improved dynamical-friction (DF) (see profiles.py for more details)
#########################################################################
from profiles import ftot
import numpy as np
from scipy.integrate import ode
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
#########################################################################
#---
class orbit(object):
"""
Class for orbit and orbit integration in axisymetric potential
Syntax:
o = orbit(xv,potential=None)
which initializes an orbit object "o", where
xv: the phase-space coordinate in a cylindrical frame
[R,phi,z,VR,Vphi,Vz] [kpc,radian,kpc,kpc/Gyr,kpc/Gyr,kpc/Gyr]
(numpy array)
potential: host potential (a density profile object, as defined
in profiles.py, or a list such objects that constitute a
composite potential)
Attributes:
o.t: integration time [Gyr] (float, list, or array)
o.xv: phase-space coordinates in a cylindrical frame
[R,phi,z,VR,Vphi,Vz] [kpc,radian,kpc,kpc/Gyr,kpc/Gyr,kpc/Gyr]
which are initialized to be the value put in by hand when
the orbit object is created, and are updated once the method
o.integrate is called
(numpy array)
o.tList: time sequence
o.xvList: coordinates along the time sequence
and conditionally also: (only available if potential is not None and
spherically symmetric)
o.rperi: peri-center distance [kpc] (float)
o.rapo: apo-center distance [kpc] (float)
Methods:
o.integrate(t,potential,m=None,CoulombLogChoice=None):
updates the coordinates o.xv by integrates over time "t" in
the "potential", using scipy.integrate.ode, and considering
dynamical friction if "m" and "CoulombLogChoice" are provided
Arthur Fangzhou Jiang, 2016-10-27, HUJI
Arthur Fangzhou Jiang, 2019-08-21, UCSC
"""
def __init__(self,xv,potential=None):
r"""
Initialize an orbit by specifying the phase-space coordinates
in a cylindrical frame.
Syntax:
o = orbit(xv, potential=None)
where
xv: phase-space coordinates in a cylindrical frame
[R,phi,z,VR,Vphi,Vz]
[kpc,radian,kpc,kpc/Gyr,kpc/Gyr,kpc/Gyr]
(numpy array)
potential: host potential (a density profile object, as
defined in profiles.py, or a list such objects that
constitute a composite potential)
(dafault=None, i.e., when initializing an orbit, do not
specify the potential, but if provided, a few more
attributes attributes are triggered, including o.rperi
and o.rapo, and maybe more, to be determined)
"""
self.xv = xv # instantaneous coordinates, initialized by input
self.t = 0. # instantaneous time, initialized here to be 0.
self.tList = [] # initialize time sequencies
self.xvList = []
if potential is not None:
pass # <<< to be added: self.rperi and self.rapo etc
def integrate(self,t,potential,m=None):
r"""
Integrate orbit over time "t" [Gyr], using methods that comes
with scipy.integrate.ode; and update coorinates to be the values
at the end of t.
Syntax:
o.integrate(t,potential,m=None)
where
t: time [Gyr] (float, list, or numpy array)
potential: host potential (a profile object, as defined
in profile.py, or a list of such objects which
altogether constitute the host potential)
m: satellite mass [Msun]
(default is None; if provided, dynamical friction is
triggered)
Note that in case when t is list or array, attributes such as
.tList
.xvList
which are lists, will start from empty and get appended new
value for each timestep; while attributes
.t
.xv
store the instantaneous time and coordinates atthe end of t.
"""
solver = ode(f,jac=None).set_integrator(
#'vode',
'dopri5',
nsteps=500, # default=500
max_step = 0.1, # default=0.0
rtol=1e-5, # default = 1e-6
atol=1e-5, # default = 1e-12
)
solver.set_initial_value(self.xv, self.t)
solver.set_f_params(potential,m,)
if isinstance(t, list) or isinstance(t,np.ndarray):
for tt in t:
solver.integrate(tt)
self.t = tt
self.xv = solver.y
self.tList.append(self.t)
self.xvList.append(self.xv)
else: # i.e., if t is a scalar
solver.integrate(t)
self.xv = solver.y
self.t = solver.t
self.tList.append(self.t)
self.xvList.append(self.xv)
self.tArray = np.array(self.tList)
self.xvArray = np.array(self.xvList)
def f(t,y,p,m):
r"""
Returns right-hand-side functions of the EOMs for orbit integration:
d R / d t = VR
d phi / d t = Vphi / R
d z / d t = Vz
d VR / dt = Vphi^2 / R + fR
d Vphi / dt = - VR * Vphi / R + fphi
d Vz / d t = fz
for use in the method ".integrate" of the "orbit" class.
Syntax:
f(t,y,p,m)
where
t: integration time [Gyr] (float)
y: phase-space coordinates in a cylindrical frame
[R,phi,z,VR,Vphi,Vz] [kpc,radian,kpc,kpc/Gyr,kpc/Gyr,kpc/Gyr]
(numpy array)
p: host potential (a density profile object, as defined
in profiles.py, or a list of such objects that constitute a
composite potential)
m: satellite mass [Msun] (float or None)
(If m is None, DF is ignored)
Note that fR, fphi, fz are the R-, phi- and z-components of the
acceleration at phase-space location y, computed by the function
"ftot" defined in profiles.py.
Return: the list of
[VR,
Vphi / R
Vz,
Vphi^2 / R + fR ,
- VR * Vphi / R + fphi,
fz]
i.e., the right-hand side of the EOMs describing the evolution of the
phase-space coordinates in a cylindrical frame
"""
R, phi, z, VR, Vphi, Vz = y
fR, fphi, fz = ftot(p,y,m)
R = max(R,1e-6) # safety
return [VR,
Vphi/R,
Vz,
Vphi**2./R + fR,
- VR*Vphi/R + fphi,
fz]