sdist.py 9.12 KB
 findesgh committed Aug 21, 2018 1 ``````############################################################################### `````` Martin Glatzle committed Jan 31, 2019 2 ``````import astropy.units as _u `````` findesgh committed Aug 21, 2018 3 ``````import numpy as _np `````` Martin Glatzle committed Feb 21, 2019 4 5 ``````import os as _os from scipy.special import erf as _erf `````` findesgh committed Aug 21, 2018 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 ``````############################################################################### class SizeDist(object): """ Grain size distribution. Computes 1/n_H * dn_gr/da [1/L], where n_H is the hydrogen nuclei number density, n_gr(a) is the number density of grains <= a and a is a grain size. Note that this is not a distribution in the mathematical sense. Since it is commonly referred to as such in the community, though, and for lack of a better term, this convention is maintained here. It might, however, be changed in the future. See [1]_ and specifically Sec. 1 for an example. Parameters ---------- sizeMin : scalar Quantity [L] `````` findesgh committed Sep 24, 2018 23 `````` Low end size cutoff of the distribution. `````` findesgh committed Aug 21, 2018 24 `````` sizeMax : scalar Quantity [L] `````` findesgh committed Sep 24, 2018 25 `````` High end size cutoff of the distribution. `````` findesgh committed Aug 21, 2018 26 `````` func : callable `````` findesgh committed Sep 24, 2018 27 28 29 `````` Must take a single Quantity with array value of sizes and return 1/n_H * dn_gr/da [1/L]. It is called in the limits set by `sizeMin` and `sizeMax` when an instance of this class is called. `````` findesgh committed Aug 21, 2018 30 31 32 33 `````` Attributes ---------- sizeMin : scalar Quantity [L] `````` findesgh committed Sep 24, 2018 34 `````` Directly taken from parameters. `````` findesgh committed Aug 21, 2018 35 `````` sizeMax : scalar Quantity [L] `````` findesgh committed Sep 24, 2018 36 `````` Directly taken from parameters. `````` findesgh committed Aug 21, 2018 37 `````` func : callable `````` findesgh committed Sep 24, 2018 38 `````` Directly taken from parameters. `````` findesgh committed Aug 21, 2018 39 40 41 42 43 44 45 46 47 48 49 50 `````` References ---------- .. [1] Weingartner & Draine 2001ApJ...548..296W Examples -------- >>> import astropy.units as u >>> def f(s): \ return 1.0/s.unit >>> a = SizeDist(3.5*u.angstrom, 1*u.micron, f) >>> a(_np.logspace(-11, -5, 10)*u.m) `````` findesgh committed Sep 24, 2018 51 `````` `````` findesgh committed Aug 21, 2018 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 `````` """ def __init__(self, sizeMin, sizeMax, func): """ See class docstring. """ self.sizeMin = sizeMin self.sizeMax = sizeMax self.func = func return def __call__(self, sizes): """ Compute 1/n_H * dn_gr/da [1/L]. Returns func(sizes) in limits set by sizeMin and sizeMax and zero elsewhere. Parameters ---------- sizes : array-valued Quantity [L] `````` findesgh committed Sep 24, 2018 72 `````` Grain sizes at which to evaluate func. `````` findesgh committed Aug 21, 2018 73 74 75 76 77 78 79 80 81 82 `````` Returns ------- r : array-valued Quantity [1/L] """ r = _np.zeros(len(sizes))/sizes.unit ind = _np.where((sizes >= self.sizeMin) & (sizes <= self.sizeMax)) r[ind] = self.func(sizes[ind]) return r `````` findesgh committed Dec 06, 2018 83 84 85 86 87 88 `````` def __add__(self, other): """ Commutative addition of size distributions. Overloading of ``+`` operator. Can add an instance of ``SizeDist`` (or subclass thereof), any callable or a scalar to ``self.func``. No check `````` Martin Glatzle committed Feb 21, 2019 89 `````` on units is performed. Returns an instance of ``SizeDist``, the `````` findesgh committed Dec 06, 2018 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 `````` ``func`` attribute of which is defined to return the corresponding sum. If ``other`` is an instance of ``SizeDist`` (or subclass), take maximum of the two ``sizeMin`` attributes and minimum of the two ``sizeMax`` attributes. Parameters ---------- other : SizeDist, callable or scalar-valued Quantity [1/L] Is added to ``self.func``. Returns ------- : ``self.__class__`` Instance of own class with corresponding ``func`` attribute. Examples -------- >>> import astropy.units as u >>> def f(s): \ return 1.0/s.unit >>> a = SizeDist(1*u.angstrom, 1*u.micron, f) >>> b = SizeDist(10*u.angstrom, 10*u.micron, f) >>> c = a + b >>> c(_np.logspace(-11, -4, 10)*u.m) """ if issubclass(other.__class__, SizeDist): # find new size limits sizeMin = max(self.sizeMin, other.sizeMin) sizeMax = min(self.sizeMax, other.sizeMax) # new differential number density is sum def func(sizes): return self.func(sizes) + other.func(sizes) `````` Martin Glatzle committed Feb 21, 2019 127 `````` return SizeDist(sizeMin, sizeMax, func) `````` findesgh committed Dec 06, 2018 128 129 130 `````` elif callable(other): def func(sizes): return self.func(sizes) + other(sizes) `````` Martin Glatzle committed Feb 22, 2019 131 `````` return SizeDist(self.sizeMin, self.sizeMax, func) `````` findesgh committed Dec 06, 2018 132 133 134 `````` else: def func(sizes): return other + self.func(sizes) `````` Martin Glatzle committed Feb 22, 2019 135 `````` return SizeDist(self.sizeMin, self.sizeMax, func) `````` findesgh committed Dec 06, 2018 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 `````` # make addition commutative __radd__ = __add__ def __mul__(self, other): if callable(other): def func(sizes): return self.function(sizes) * other(sizes) return self.__class__(self.sizeMin, self.sizeMax, func) else: def func(sizes): return other * self.function(sizes) return self.__class__(self.sizeMin, self.sizeMax, func) # make multiplication commutative __rmul__ = __mul__ `````` findesgh committed Jan 08, 2019 154 155 156 157 158 159 160 161 162 163 ``````class PowerLaw(SizeDist): """ Parameters ---------- sizeMin : scalar Quantity [L] Low end size cutoff of the distribution. sizeMax : scalar Quantity [L] High end size cutoff of the distribution. power : float Log-slope of the size distribution. `````` findesgh committed Jan 24, 2019 164 `````` C : scalar Quantity [L**(-1-power)] `````` findesgh committed Jan 08, 2019 165 166 167 168 169 `````` Normalization of the size distribution. """ def __init__(self, sizeMin, sizeMax, power, C): def f(a): return C*a**power `````` Martin Glatzle committed Feb 20, 2019 170 `````` super().__init__(sizeMin, sizeMax, f) `````` findesgh committed Jan 08, 2019 171 172 `````` return `````` findesgh committed Aug 21, 2018 173 `````` `````` Martin Glatzle committed Feb 21, 2019 174 175 176 ``````class DoubleLogNormal(SizeDist): """ Sum of two log normal distributions as in Eq. (2) of [1]_. `````` Martin Glatzle committed Jan 31, 2019 177 `````` `````` Martin Glatzle committed Feb 21, 2019 178 179 180 181 182 `````` References ---------- .. [1] Weingartner & Draine 2001ApJ...548..296W """ def __init__(self, sizeMin, sizeMax, rho, sgma, bc, a0): `````` Martin Glatzle committed Jan 31, 2019 183 184 185 186 187 188 189 190 `````` # mass of carbon atom m_C = 12.0107*_u.u def B_val(i): nominator = (3 * _np.exp(-4.5 * sgma**2) * bc[i] * m_C) denominator = ( (2*_np.pi**2)**1.5 * rho * a0[i]**3 * sgma * `````` Martin Glatzle committed Feb 20, 2019 191 `````` (1 + `````` Martin Glatzle committed Feb 21, 2019 192 193 194 `````` _erf((3 * sgma/_np.sqrt(2)) + (_np.log((a0[i]/3.5/_u.angstrom).decompose().value) / (sgma * _np.sqrt(2))) `````` Martin Glatzle committed Feb 20, 2019 195 196 197 `````` ) ) ) `````` Martin Glatzle committed Jan 31, 2019 198 199 200 201 202 203 204 205 `````` # FIXME: what do we do if the denominator is zero if denominator != 0: B = (nominator/denominator).decompose().value return B _B = [B_val(i) for i in range(2)] def f(a): `````` Martin Glatzle committed Feb 21, 2019 206 `````` return sum([_B[i]/a * _np.exp(-0.5*( `````` Martin Glatzle committed Jan 31, 2019 207 208 209 `````` _np.log((a/a0[i]).decompose().value)/sgma)**2) for i in range(2)]) `````` Martin Glatzle committed Feb 20, 2019 210 `````` super().__init__(sizeMin, sizeMax, f) `````` Martin Glatzle committed Jan 31, 2019 211 212 213 `````` return `````` Martin Glatzle committed Feb 21, 2019 214 215 216 217 ``````class WD01ExpCutoff(SizeDist): """ Power law distribution with exponential cutoff as in Eqs. (4) and (5) of [1]_. `````` Martin Glatzle committed Jan 31, 2019 218 `````` `````` Martin Glatzle committed Feb 21, 2019 219 220 221 222 223 `````` References ---------- .. [1] Weingartner & Draine 2001ApJ...548..296W """ def __init__(self, sizeMin, sizeMax, alpha, beta, a_t, a_c, C): `````` Martin Glatzle committed Jan 31, 2019 224 `````` `````` Martin Glatzle committed Feb 21, 2019 225 226 227 228 229 230 `````` if beta >= 0: def F(a): return 1 + (beta * a) / a_t else: def F(a): return 1 / (1 - (beta * a) / a_t) `````` Martin Glatzle committed Jan 31, 2019 231 `````` `````` Martin Glatzle committed Feb 21, 2019 232 233 234 235 236 `````` def exp_func(a): r = _np.ones_like(a.value) ind = _np.where(a > a_t) r[ind] = _np.exp(-(((a[ind] - a_t)/a_c).decompose().value)**3) return r `````` Martin Glatzle committed Jan 31, 2019 237 `````` `````` Martin Glatzle committed Feb 21, 2019 238 239 240 241 242 243 `````` def f(a): return C*(a / a_t)**alpha/a \ * F(a) * exp_func(a) super().__init__(sizeMin, sizeMax, f) return `````` Martin Glatzle committed Jan 31, 2019 244 245 `````` `````` Martin Glatzle committed Feb 21, 2019 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 ``````def WD01(RV, bc, case): # pandas dataframe would be far better suited for this but don't know if # pulling in another dependency makes sense just for this data = _np.genfromtxt( _os.path.join( _os.path.dirname(__file__), 'data/sdist/WD01-2001ApJ...548..296W-Table1.txt' ), dtype=None, encoding=None ) # doing == float comparison here, might come back to bite us at some point params = data[ [ j for j, x in enumerate(data) if (x[0] == RV and x[1] == bc and x[3] == case) ] ] if len(params) > 1: # this should never happen raise ValueError("Could not uniquely identify parameter set.") elif len(params) == 0: raise ValueError("Could not identify parameter set.") params = params[0] sizeMin = 3.5*_u.angstrom sizeMax = 10*_u.micron rho = 2.24*_u.g/_u.cm**3 `````` Martin Glatzle committed Feb 21, 2019 276 `````` s_car = DoubleLogNormal( `````` Martin Glatzle committed Feb 21, 2019 277 278 279 280 281 282 283 `````` sizeMin, sizeMax, rho, 0.4, [0.75*params[2], 0.25*params[2]], [3.5*_u.angstrom, 30*_u.angstrom] ) `````` Martin Glatzle committed Feb 21, 2019 284 `````` l_car = WD01ExpCutoff( `````` Martin Glatzle committed Feb 21, 2019 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 `````` sizeMin, sizeMax, params[4], params[5], params[6]*_u.m, params[7]*_u.m, params[8] ) sil = WD01ExpCutoff( sizeMin, sizeMax, params[9], params[10], params[11]*_u.m, params[12]*_u.m, params[13] ) `````` Martin Glatzle committed Feb 21, 2019 302 `````` return s_car + l_car, sil `````` findesgh committed Aug 21, 2018 303 `````` `````` Martin Glatzle committed Feb 23, 2019 304 `````` `````` findesgh committed Aug 21, 2018 305 306 307 308 309 ``````############################################################################### if __name__ == "__main__": import doctest doctest.testmod() ###############################################################################``````