sdist.py 9.12 KB
Newer Older
1
###############################################################################
2
import astropy.units as _u
3
import numpy as _np
Martin Glatzle's avatar
Martin Glatzle committed
4
5
import os as _os
from scipy.special import erf as _erf
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's avatar
findesgh committed
23
        Low end size cutoff of the distribution.
24
    sizeMax : scalar Quantity [L]
findesgh's avatar
findesgh committed
25
        High end size cutoff of the distribution.
26
    func : callable
findesgh's avatar
findesgh committed
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.
30
31
32
33

    Attributes
    ----------
    sizeMin : scalar Quantity [L]
findesgh's avatar
findesgh committed
34
        Directly taken from parameters.
35
    sizeMax : scalar Quantity [L]
findesgh's avatar
findesgh committed
36
        Directly taken from parameters.
37
    func : callable
findesgh's avatar
findesgh committed
38
        Directly taken from parameters.
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's avatar
findesgh committed
51
    <Quantity [0., 0., 0., 1., 1., 1., 1., 1., 0., 0.] 1 / m>
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's avatar
findesgh committed
72
            Grain sizes at which to evaluate func.
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's avatar
findesgh committed
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's avatar
Martin Glatzle committed
89
        on units is performed. Returns an instance of ``SizeDist``, the
findesgh's avatar
findesgh committed
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)
        <Quantity [0., 0., 0., 2., 2., 2., 2., 0., 0., 0.] 1 / 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's avatar
Martin Glatzle committed
127
            return SizeDist(sizeMin, sizeMax, func)
findesgh's avatar
findesgh committed
128
129
130
        elif callable(other):
            def func(sizes):
                return self.func(sizes) + other(sizes)
131
            return SizeDist(self.sizeMin, self.sizeMax, func)
findesgh's avatar
findesgh committed
132
133
134
        else:
            def func(sizes):
                return other + self.func(sizes)
135
            return SizeDist(self.sizeMin, self.sizeMax, func)
findesgh's avatar
findesgh committed
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's avatar
findesgh committed
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's avatar
findesgh committed
164
    C : scalar Quantity [L**(-1-power)]
findesgh's avatar
findesgh committed
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's avatar
Martin Glatzle committed
170
        super().__init__(sizeMin, sizeMax, f)
findesgh's avatar
findesgh committed
171
172
        return

173

Martin Glatzle's avatar
Martin Glatzle committed
174
175
176
class DoubleLogNormal(SizeDist):
    """
    Sum of two log normal distributions as in Eq. (2) of [1]_.
177

Martin Glatzle's avatar
Martin Glatzle committed
178
179
180
181
182
    References
    ----------
    .. [1] Weingartner & Draine 2001ApJ...548..296W
    """
    def __init__(self, sizeMin, sizeMax, rho, sgma, bc, a0):
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's avatar
Martin Glatzle committed
191
                (1 +
Martin Glatzle's avatar
Martin Glatzle committed
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's avatar
Martin Glatzle committed
195
196
197
                 )
                )
            )
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's avatar
Martin Glatzle committed
206
            return sum([_B[i]/a * _np.exp(-0.5*(
207
208
209
                _np.log((a/a0[i]).decompose().value)/sgma)**2)
                        for i in range(2)])

Martin Glatzle's avatar
Martin Glatzle committed
210
        super().__init__(sizeMin, sizeMax, f)
211
212
213
        return


Martin Glatzle's avatar
Martin Glatzle committed
214
215
216
217
class WD01ExpCutoff(SizeDist):
    """
    Power law distribution with exponential cutoff as in Eqs. (4) and (5) of
    [1]_.
218

Martin Glatzle's avatar
Martin Glatzle committed
219
220
221
222
223
    References
    ----------
    .. [1] Weingartner & Draine 2001ApJ...548..296W
    """
    def __init__(self, sizeMin, sizeMax, alpha, beta, a_t, a_c, C):
224

Martin Glatzle's avatar
Martin Glatzle committed
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)
231

Martin Glatzle's avatar
Martin Glatzle committed
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
237

Martin Glatzle's avatar
Martin Glatzle committed
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
244
245


Martin Glatzle's avatar
Martin Glatzle committed
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's avatar
Martin Glatzle committed
276
    s_car = DoubleLogNormal(
Martin Glatzle's avatar
Martin Glatzle committed
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's avatar
Martin Glatzle committed
284
    l_car = WD01ExpCutoff(
Martin Glatzle's avatar
Martin Glatzle committed
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's avatar
Martin Glatzle committed
302
    return s_car + l_car, sil
303

Martin Glatzle's avatar
PEP8    
Martin Glatzle committed
304

305
306
307
308
309
###############################################################################
if __name__ == "__main__":
    import doctest
    doctest.testmod()
###############################################################################