...

Commits (4)
 ... ... @@ -171,9 +171,9 @@ class PowerLaw(SizeDist): return class DoubleLogNormal(SizeDist): class LogNormal(SizeDist): """ Sum of two log normal distributions as in Eq. (2) of [1]_. Log normal distribution as in Eq. (2) of [1]_. References ---------- ... ... @@ -184,28 +184,23 @@ class DoubleLogNormal(SizeDist): # 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 * (1 + _erf((3 * sgma/_np.sqrt(2)) + (_np.log((a0[i]/3.5/_u.angstrom).decompose().value) / (sgma * _np.sqrt(2))) ) ) nominator = (3 * _np.exp(-4.5 * sgma**2) * bc * m_C) denominator = ( (2*_np.pi)**1.5 * rho * a0**3 * sgma * (1 + _erf((3 * sgma/_np.sqrt(2)) + (_np.log(a0/3.5/_u.angstrom) / (sgma * _np.sqrt(2))) ) ) # 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)] ) # FIXME: what do we do if the denominator is zero if denominator != 0: B = (nominator/denominator).decompose() def f(a): return sum([_B[i]/a * _np.exp(-0.5*( _np.log((a/a0[i]).decompose().value)/sgma)**2) for i in range(2)]) return B/a * _np.exp(-0.5*( _np.log(a/a0)/sgma)**2) super().__init__(sizeMin, sizeMax, f) return ... ... @@ -224,19 +219,19 @@ class WD01ExpCutoff(SizeDist): if beta >= 0: def F(a): return 1 + (beta * a) / a_t return 1.0 + (beta * a) / a_t else: def F(a): return 1 / (1 - (beta * a) / a_t) return 1.0 / (1 - (beta * a) / a_t) 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) r[ind] = _np.exp(-(((a[ind] - a_t)/a_c))**3) return r def f(a): return C*(a / a_t)**alpha/a \ return C/a*(a / a_t)**alpha \ * F(a) * exp_func(a) super().__init__(sizeMin, sizeMax, f) ... ... @@ -273,13 +268,21 @@ def WD01(RV, bc, case): sizeMin = 3.5*_u.angstrom sizeMax = 10*_u.micron rho = 2.24*_u.g/_u.cm**3 s_car = DoubleLogNormal( s_car = LogNormal( sizeMin, sizeMax, rho, 0.4, 0.75*params[2], 3.5*_u.angstrom ) + \ LogNormal( sizeMin, sizeMax, rho, 0.4, [0.75*params[2], 0.25*params[2]], [3.5*_u.angstrom, 30*_u.angstrom] 0.25*params[2], 30*_u.angstrom ) l_car = WD01ExpCutoff( sizeMin, ... ...
 ... ... @@ -81,6 +81,53 @@ class TestSdist(TestCase): return class TestWD01ExpCutoff(TestCase): def test_simple(self): sd = sdist.WD01ExpCutoff( 3.5*u.angstrom, 1*u.micron, 0, 0, 1*u.micron, 1*u.micron, 1) sizes = [ 1*u.micron, 1*u.nm, 1*u.angstrom, 2*u.micron, ] res = [ 1./(1*u.micron), 1./(1*u.nm), 0, 0, ] for s, r in zip(sizes, res): self.assertEqual( sd(np.array([s.value])*s.unit)[0], r ) return def test_simple_2(self): sd = sdist.WD01ExpCutoff( 3.5*u.angstrom, 1*u.micron, 0, 0, 1*u.nm, 1*u.micron, 1) sizes = [ 1*u.micron, 2*u.nm, 0.5*u.nm, 1*u.angstrom, 2*u.micron, ] res = [ 1./(1*u.micron) * np.exp(-(1-1e-3)**3), 1./(2*u.nm), 2/u.nm * np.exp(1e-9/8), u.Quantity(0), u.Quantity(0), ] for s, r in zip(sizes, res): self.assertAlmostEqual( sd(np.array([s.value])*s.unit)[0].value, r.value ) return class TestWD01(TestCase): def test_smoke_test(self): ... ...