forked from ibackus/ICgen
-
Notifications
You must be signed in to change notification settings - Fork 1
/
sigma_profile.py
288 lines (205 loc) · 7.6 KB
/
sigma_profile.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
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 23 10:17:53 2014
@author: ibackus
"""
import numpy as np
import pynbody
SimArray = pynbody.array.SimArray
import isaac
def make_profile(ICobj):
"""
A wrapper for generating surface density profiles according to the IC object.
Settings for the profile are defined in ICobj.settings. Which profile gets
used is defined by ICobj.settings.sigma.kind
Currently available kinds are:
viscous
powerlaw
MQWS
**RETURNS**
r : SimArray
Radii at which sigma is calculated
sigma : SimArray
Surface density profile as a function of R
"""
kind = ICobj.settings.sigma.kind
if kind == 'powerlaw':
r, sigma = powerlaw(ICobj.settings, ICobj.T)
elif (kind == 'mqws') | (kind == 'MQWS'):
r, sigma = MQWS(ICobj.settings, ICobj.T)
elif (kind == 'viscous'):
r, sigma = viscous(ICobj.settings)
else:
raise TypeError, 'Could not make profile for kind {0}'.format(kind)
if hasattr(ICobj.settings.sigma, 'innercut'):
sigma = _applycut(r, sigma, ICobj.settings.sigma.innercut, False)
if hasattr(ICobj.settings.sigma, 'outercut'):
sigma = _applycut(r, sigma, ICobj.settings.sigma.outercut, True)
return r, sigma
def _applycut(r, sigma, rcut, outer=True):
"""
Applies a hard cut to a surface density profile (sigma). If outer=True,
sigma = 0 at r > rcut. Otherwise, sigma = 0 at r < rcut. If rcut is
None, inf, or nan no cut is performed.
"""
if rcut is None:
return sigma
elif np.isnan(rcut) or np.isinf(rcut):
return sigma
if outer:
mask = r > rcut
else:
mask = r < rcut
if np.any(mask):
sigma[mask] = 0
return sigma
def viscous(settings):
"""
Generates a surface density profile derived from a self-similarity solution
for a viscous disk, according to:
sigma ~ r^-gamma exp(-r^(2-gamma))
Where r is a dimensionless radius and gamma is a constant less than 2.
Rd (disk radius) is defined as the radius containing 95% of the disk mass
**ARGUMENTS**
settings : IC settings
settings like those contained in an IC object (see ICgen_settings.py)
**RETURNS**
R : SimArray
Radii at which sigma is calculated
sigma : SimArray
Surface density profile as a function of R
"""
Rd = settings.sigma.Rd
rin = settings.sigma.rin
rmax = settings.sigma.rmax
n_points = settings.sigma.n_points
gamma = settings.sigma.gamma
m_disk = settings.sigma.m_disk
# Define the fraction of mass contained within Rd
A = 0.95
# Normalization for r
R1 = Rd / (np.log(1/(1-A))**(1/(2-gamma)))
Rmax = rmax * Rd
Rin = rin * Rd
R = np.linspace(0, Rmax, n_points)
r = (R/R1).in_units('1')
sigma = (r**-gamma) * np.exp(-r**(2-gamma)) * (m_disk/(2*np.pi*R1*R1)) * (2-gamma)
# Deal with infinities at the origin with a hard cut off
sigma[0] = sigma[1]
# Apply interior cutoff
cut_mask = R < Rin
if np.any(cut_mask):
sigma[cut_mask] *= isaac.smoothstep(r[cut_mask],degree=21,rescale=True)
return R, sigma
def powerlaw(settings, T = None):
"""
Generates a surface density profile according to a powerlaw sigma ~ 1/r
with a smooth interior cutoff and smooth exterior exponential cutoff.
**ARGUMENTS**
settings : IC settings
settings like those contained in an IC object (see ICgen_settings.py)
T : callable function
Function that returns temperature of the disk as a function of radius
IF none, a powerlaw temperature is assumed
**RETURNS**
R : SimArray
Radii at which sigma is calculated
sigma : SimArray
Surface density profile as a function of R
"""
# Parse settings
Rd = settings.sigma.Rd
rin = settings.sigma.rin
rmax = settings.sigma.rmax
cutlength = settings.sigma.cutlength
Mstar = settings.physical.M
Qmin = settings.sigma.Qmin
n_points = settings.sigma.n_points
m = settings.physical.m
if T is None:
# If no callable object to calculate Temperature(R) is provided,
# default to a powerlaw T ~ R^-q
T0 = SimArray([129.0],'K') # Temperature at 1 AU
R0 = SimArray([1.0],'au')
q = 0.59
def T(x):
return T0 * np.power((x/R0).in_units('1'),-q)
Rd = isaac.match_units(pynbody.units.au, Rd)[1]
Mstar = isaac.match_units(pynbody.units.Msol, Mstar)[1]
# Molecular weight
m = isaac.match_units(m, pynbody.units.m_p)[0]
# Maximum R to calculate sigma at (needed for the exponential cutoff region)
Rmax = rmax*Rd
# Q calculation parameters:
G = SimArray([1.0],'G')
kB = SimArray([1.0],'k')
# Initialize stuff
A = SimArray(1.0,'Msol')/(2*np.pi*np.power(Rd,2))
R = np.linspace(0,Rmax,n_points)
r = np.array((R/Rd).in_units('1'))
# Calculate sigma
# Powerlaw
sigma = A/r
sigma[0] = 0.0
# Interior cutoff
sigma[r>1] *= np.exp(-(r[r>1] - 1)**2 / (2*cutlength**2))
# Exterior cutoff
sigma[r<rin] *= isaac.smoothstep(r[r<rin],degree=21,rescale=True)
# Calculate Q
Q = np.sqrt(Mstar*kB*T(R)/(G*m*R**3))/(np.pi*sigma)
Q.convert_units('1')
# Rescale sigma to meet the minimum Q requirement
sigma *= Q.min()/Qmin
# Calculate Q
Q = np.sqrt(Mstar*kB*T(R)/(G*m*R**3))/(np.pi*sigma)
Q.convert_units('1')
return R, sigma
#def MQWS(n_points=1000, rin=4.0, rout=20.0, rmax = None, m_disk=0.1):
def MQWS(settings, T):
"""
Generates a surface density profile as the per method used in Mayer, Quinn,
Wadsley, and Stadel 2004
** ARGUMENTS **
NOTE: if units are not supplied, assumed units are AU, Msol
settings : IC settings
settings like those contained in an IC object (see ICgen_settings.py)
T : callable
A function to calculate temperature as a function of radius
** RETURNS **
r : SimArray
Radii at which sigma is calculated
sigma : SimArray
Surface density profile as a function of R
"""
# Q calculation parameters:
G = SimArray([1.0],'G')
kB = SimArray([1.0],'k')
# Load in settings
n_points = settings.sigma.n_points
rin = settings.sigma.rin
rout = settings.sigma.rout
rmax = settings.sigma.rmax
Qmin = settings.sigma.Qmin
m = settings.physical.m
Mstar = settings.physical.M
#m_disk = settings.sigma.m_disk
rin = isaac.match_units(pynbody.units.au, rin)[1]
rout = isaac.match_units(pynbody.units.au, rout)[1]
#m_disk = isaac.match_units(pynbody.units.Msol, m_disk)[1]
if rmax is None:
rmax = 2.5 * rout
else:
rmax = isaac.match_units(pynbody.units.au, rmax)[1]
r = np.linspace(0, rmax, n_points)
#A = m_disk * np.exp(2 * (rin/rout).in_units('1'))/(rout * np.pi**1.5)
a = (rin/r).in_units('1')
b = (r/rout).in_units('1')
#sigma = A * np.exp(-a**2 - b**2)/r
sigma = (np.exp(-a**2 - b**2)/r) * Mstar.units/r.units
# Calculate Q
Q = np.sqrt(Mstar*kB*T(r)/(G*m*r**3))/(np.pi*sigma)
Q.convert_units('1')
sigma *= np.nanmin(Q)/Qmin
# Remove all nans
sigma[np.isnan(sigma)] = 0.0
return r, sigma