1
2
3 """
4 Module for calculating density structures and molecular abundances as a
5 function of radius.
6
7 Author: R. Lombaert
8
9 """
10
11 import numpy as np
12 from scipy.interpolate import InterpolatedUnivariateSpline as spline1d
13 from astropy import constants as cst
14 from astropy import units as u
15
16 from cc.modeling.profilers import Velocity, Profiler, Mdot
17
18
20
21 '''
22 Calculate the mean molecular weight, given the fractional abundances
23 fH = n(H)/n(H2)
24 and
25 fHe = n(He)/n(Htot)
26 with
27 n(Htot) = n(H) + n(H2)
28
29 Derived from: rho_gas / n(gas) = mu * m_H
30
31 @param fH: The fractional abundance of atomic hydrogen with respect to
32 molecular hydrogen. fH = n(H)/n(H2). Can be a constant or
33 a radially dependent value.
34 @type fH: float/array
35 @param fHe: The fractional abundance of helium with respect to total
36 hydrogen number density. fHe = n(He)/(n(H)+2*n(H2)). Can be
37 a constant or a radially dependent value.
38 @type fHe: float/array
39
40 @return: The mean molecular weight in units of m_H
41 @rtype: float
42
43 '''
44
45 mu = (fH+2.+4.*fHe*(fH+2.))/(fH+1.+fHe*(fH+2.))
46
47 return mu
48
49
50
52
53 '''
54 Calculate the mass density profile based on mass-loss rate (Msun/yr) and
55 expansion velocity (cm/s).
56
57 Based on conservation of mass.
58
59 @param r: The radial points (cm)
60 @type r: array
61 @param mdot: The mass-loss-rate profile or as a single value (constant mass
62 loss) (Msun/yr)
63 @type mdot: Mdot()/float
64 @param v: The velocity profile or as a single value (constant velocity)
65 (cm/s)
66 @type v: Velocity()/float
67
68 @return: The density profile (g/cm3)
69 @rtype: array/float
70
71 '''
72
73
74
75 vi = v.eval(r) if isinstance(v,Velocity.Velocity) else v
76 mdoti = mdot.eval(r) if isinstance(mdot,Mdot.Mdot) else mdot
77
78
79 denominator = 4.*np.pi*r**2.*vi
80 nominator = (mdoti*u.Msun/u.yr).to(u.g/u.s).value
81 return nominator/denominator
82
83
84
86
87 '''
88 An interface for the density profile as a function of radius.
89
90 '''
91
94
95 '''
96 Create an instance of the Density() class. Requires a radial grid, a
97 function for calculating the density profile, and its arguments.
98
99 The default is the mass-loss rate density profile. Requires mdot and v
100 to be given either as constants or Profile()s.
101
102 @param r: The radial points (cm)
103 @type r: array
104
105 @keyword func: The function for calculating the density profile.
106
107 (default: dens_mdot)
108 @type func: func
109 @keyword dfunc: The function that describes the derivative of the
110 profile with respect to r. Can be given as an interp1d
111 object. If None, a generic central difference is taken
112 and interpolated.
113
114 (default: None)
115 @type dfunc: function/interp1d object
116 @keyword order: Order of the spline interpolation. Only relevant if
117 spline is requested instead of linear. Default is cubic
118 interpolation.
119
120 (default: 3)
121 @type order: int
122 @keyword dust: This is a dust density profile rather than a gas density
123
124 (default: 0)
125 @type dust: bool
126
127 @keyword args: Additional parameters passed to the functions when eval
128 or diff are called.
129
130 (default: [])
131 @type args: tuple
132 @keyword kwargs: Additional keywords passed to the functions when eval
133 or diff are called.
134
135 (default: {})
136 @type kwargs: dict
137
138 '''
139
140
141 super(Density, self).__init__(r,func,dfunc,order,*args,**kwargs)
142
143 self.r = self.x
144 self.rho = self.y
145 self.drhodr = self.dydx
146
147 self.dust = dust
148
149
150 self.mh = cst.m_p.cgs.value
151
152
153 self.n = dict()
154 self.fac = dict()
155
156
157
159
160 '''
161 Return the mean molecular weight for this Density profile.
162
163 Requires calcNumberDensity to have been called before.
164
165 @return: The mean molecular weight in units of m_H
166 @rtype: float
167
168 '''
169
170 return self.mu
171
172
173
175
176 '''
177 In case of gas, calculate the H_2, H, He and total number density
178 profiles. In case of dust, calculate the dust number density.
179
180 Takes into account fractional abundances of atomic hydrogen and helium,
181 and the specific density/average grain size, respectively.
182
183 Note the specific definitions of fH and fHe. These are implemented this
184 way to work well with the energy balance calculation and the definition
185 of dT/dR as well as H_dg.
186
187 @keyword fH: The fractional abundance of atomic hydrogen with respect to
188 molecular hydrogen. fH = n(H)/n(H2). Can be a constant or
189 a radially dependent value.
190
191 (default: 0)
192 @type fH: float/array
193 @keyword fHe: The fractional abundance of helium with respect to total
194 hydrogen number density. fHe = n(He)/(n(H)+2*n(H2)). Can be
195 a constant or a radially dependent value.
196
197 (default: 0)
198 @type fHe: float/array
199 @keyword sd: The average specific density of the dust grains in g/cm3.
200 Only relevant if this is a dust density profile.
201
202 (default: None)
203 @type sd: float
204 @keyword a: The average grain size (cm)
205
206 (default: None)
207 @type a: float
208
209 '''
210
211 if self.dust:
212 self.sd = sd
213 self.a = a
214 denominator = sd*4./3.*np.pi*a**3
215 self.n['Dust'] = self.rho/denominator
216 self.fac['Dust'] = 1./denominator
217 return
218
219
220 self.fH = fH
221 self.fHe = fHe
222
223
224 self.mu = mean_molecular_weight(fH,fHe)
225
226
227
228 self.n['Gas'] = self.rho/self.mh/self.mu
229 self.fac['Gas'] = 1./self.mh/self.mu
230
231
232 denominator = self.mh*(fH+2.)*(1.+4.*fHe)
233 self.n['H2'] = self.rho/denominator
234 self.fac['H2'] = 1./denominator
235
236
237 self.n['H'] = fH*self.n['H2']
238 self.fac['H'] = fH/denominator
239
240
241 self.n['Htot'] = self.n['H'] + 2.*self.n['H2']
242 self.fac['Htot'] = fH/denominator + 2./denominator
243
244
245 self.n['He'] = fHe*self.n['Htot']
246 self.fac['He'] = fHe*(fH/denominator + 2./denominator)
247
248
249
250 - def eval(self,r=None,dtype='dens',warn=1):
251
252 '''
253 Evaluate the density function at a coordinate point.
254
255 r can be any value or array. If func is an interp1d object, it is
256 in principle limited by the r-range of the interpolator. However,
257 linear extrapolation is enabled, but it is advised not to extend much
258 beyond the given r-range.
259
260 @keyword r: The coordinate point(s). If None, the default
261 coordinate grid is used.
262
263 (default: None)
264 @type r: array/float
265 @keyword dtype: The type of density requested. Default is the mass
266 density. Alternatives: nh2, nh, nhtot, nhe, ngas.
267
268 (default: dens)
269 @type dtype: str
270 @keyword warn: Warn when extrapolation occurs.
271
272 (default: 1)
273 @type warn: bool
274
275 @return: The profile evaluated at x
276 @rtype: array/float
277
278 '''
279
280
281 dtype = dtype.lower()
282 rho = super(Density,self).eval(x=r,warn=warn)
283 if not dtype in ['nh2', 'nh', 'nhtot', 'nhe', 'ngas', 'ndust']:
284 return rho
285
286
287 if not self.n:
288 raise ValueError('Number densities not yet calculated.')
289
290
291 dtype = dtype[1:].capitalize()
292 if r is None:
293 return self.n[dtype]
294
295
296 return self.fac[dtype]*rho
297
298
299
301
302 '''
303 An interface for the dust-to-gas ratio profile as a function of radius.
304
305 '''
306
309
310 '''
311 Create an instance of the DustToGas() class. Requires a radial grid, a
312 function, and its arguments.
313
314 Default is a constant dust-to-gas ratio, in which case d2g must be given
315
316 @param r: The radial points (cm)
317 @type r: array
318 @keyword func: The function for calculating the density profile.
319
320 (default: dens_mdot)
321 @type func: func
322 @keyword dfunc: The function that describes the derivative of the
323 profile with respect to r. Can be given as an interp1d
324 object. If None, a generic central difference is taken
325 and interpolated.
326
327 (default: None)
328 @type dfunc: function/interp1d object
329 @keyword order: Order of the spline interpolation. Default is cubic.
330
331 (default: 3)
332 @type order: int
333
334 @keyword args: Additional parameters passed to the functions when eval
335 or diff are called.
336
337 (default: [])
338 @type args: tuple
339 @keyword kwargs: Additional keywords passed to the functions when eval
340 or diff are called.
341
342 (default: {})
343 @type kwargs: dict
344
345 '''
346
347
348 super(DustToGas, self).__init__(r,func,dfunc,order,*args,**kwargs)
349
350 self.r = self.x
351 self.d2g = self.y
352 self.dd2gdr = self.dydx
353