...
 
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):
......