Skip to content

Commit

Permalink
Support for a log-normal size distirbution with sigma_g=1 (i.e. a del…
Browse files Browse the repository at this point in the history 
…ta function)
  • Loading branch information
Tomas Stolker committed Dec 16, 2020
1 parent 1974ebe commit 5da6aa0
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 77 deletions.
17 changes: 11 additions & 6 deletions species/plot/plot_mcmc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -595,16 +595,21 @@ def plot_size_distributions(tag: str,

for i in range(samples.shape[0]):
if 'lognorm_radius' in box.parameters:
dn_dr, _, radii = dust_util.log_normal_distribution(10.**log_r_g[i], sigma_g[i], 1000)
dn_grains, r_width, radii = \
dust_util.log_normal_distribution(10.**log_r_g[i], sigma_g[i], 1000)

# Set the number density to zero for grain radii smaller than 1 nm
dn_dr[radii < 1e-3] = 0.
# Exclude radii smaller than 1 nm
indices = np.argwhere(radii >= 1e-3)

dn_grains = dn_grains[indices]
r_width = r_width[indices]
radii = radii[indices]

elif 'powerlaw_max' in box.parameters:
dn_dr, _, radii = dust_util.power_law_distribution(
exponent[i], 1e-3, 10.**r_max[i], 1000)
dn_grains, r_width, radii = \
dust_util.power_law_distribution(exponent[i], 1e-3, 10.**r_max[i], 1000)

ax.plot(radii, dn_dr, ls='-', lw=0.5, color='black', alpha=0.5)
ax.plot(radii, dn_grains/r_width, ls='-', lw=0.5, color='black', alpha=0.5)

plt.savefig(os.getcwd()+'/'+output, bbox_inches='tight')
plt.clf()
Expand Down
112 changes: 41 additions & 71 deletions species/util/dust_util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,64 +69,42 @@ def log_normal_distribution(radius_g: float,
Returns
-------
np.ndarray
Number of grains per radius bin, normalized to an integrated value of 1 grain.
Number of grains in each radius bin, normalized to a total of 1 grain.
np.ndarray
Widths of the radius bins (um).
np.ndarray
Grain radii (um).
"""

# # Create bins across a broad radius range to make sure that the full distribution is captured
# r_test = np.logspace(-20., 20., 1000) # (um)
#
# # Create a size distribution for extracting the approximate minimum and maximum radius
# dn_dr = np.exp(-np.log(r_test/radius_g)**2./(2.*np.log(sigma_g)**2.)) / \
# (r_test*np.sqrt(2.*np.pi)*np.log(sigma_g))
#
# # Select the radii for which dn/dr is larger than 0.1% of the maximum dn/dr
# indices = np.where(dn_dr/np.amax(dn_dr) > 0.001)[0]
#
# # Create bin boundaries (um), so +1 because there are n_sizes+1 bin boundaries
# r_bins = np.logspace(np.log10(r_test[indices[0]]),
# np.log10(r_test[indices[-1]]),
# n_bins+1) # (um)
#
# # Width of the radius bins (um)
# r_width = np.diff(r_bins)
#
# # Grains radii (um) at which the size distribution is sampled
# radii = (r_bins[1:]+r_bins[:-1])/2.
#
# # Number of grains per radius bin, normalized to an integrated value of 1 grain
# dn_dr = np.exp(-np.log(radii/radius_g)**2./(2.*np.log(sigma_g)**2.)) / \
# (radii*np.sqrt(2.*np.pi)*np.log(sigma_g))
#
# # Total of grains to one
# # n_grains = np.sum(r_width*dn_dr)
#
# # Normalize the size distribution to 1 grain
# # dn_dr /= n_grains
#
# return dn_dr, r_width, radii

# Get the radius interval which contains 99.999% of the distribution
interval = lognorm.interval(1.-1e-5, np.log(sigma_g), loc=0., scale=radius_g)

# Create bin boundaries (um), so +1 because there are n_bins+1 bin boundaries
r_bins = np.logspace(np.log10(interval[0]), np.log10(interval[1]), n_bins+1)
if sigma_g == 1.:
# The log-normal distribution is equal to a delta function with sigma_g = 1
radii = np.array([radius_g])
r_width = np.array([np.nan])
dn_grains = np.array([1.])

# Width of the radius bins (um)
r_width = np.diff(r_bins)
else:
# Get the radius interval which contains 99.999% of the distribution
interval = lognorm.interval(1.-1e-5, np.log(sigma_g), loc=0., scale=radius_g)

# Grain radii (um) at which the size distribution is sampled
radii = (r_bins[1:]+r_bins[:-1])/2.
# Create bin boundaries (um), so +1 because there are n_bins+1 bin boundaries
r_bins = np.logspace(np.log10(interval[0]), np.log10(interval[1]), n_bins+1)

# Width of the radius bins (um)
r_width = np.diff(r_bins)

# Grain radii (um) at which the size distribution is sampled
radii = (r_bins[1:]+r_bins[:-1])/2.

# Number of grains per radius bin, normalized to an integrated value of 1 grain
# The log-normal distribution from Ackerman & Marley 2001 gives the same
# result as scipy.stats.lognorm.pdf with s = log(sigma_g) and scale=radius_g
dn_dr = lognorm.pdf(radii, s=np.log(sigma_g), loc=0., scale=radius_g)
# Number of grains per radius bin width, normalized to an integrated value of
# 1 grain, that is, np.sum(dn_dr*r_width) = 1
# The log-normal distribution from Ackerman & Marley 2001 gives the same
# result as scipy.stats.lognorm.pdf with s = log(sigma_g) and scale=radius_g
dn_dr = lognorm.pdf(radii, s=np.log(sigma_g), loc=0., scale=radius_g)

return dn_dr, r_width, radii
# Number of grains for each radius bin
dn_grains = dn_dr*r_width

return dn_grains, r_width, radii


@typechecked
Expand All @@ -151,7 +129,7 @@ def power_law_distribution(exponent: float,
Returns
-------
np.ndarray
Number of grains per radius bin, normalized to an integrated value of 1 grain.
Number of grains in each radius bin, normalized to a total of 1 grain.
np.ndarray
Widths of the radius bins (um).
np.ndarray
Expand All @@ -167,21 +145,20 @@ def power_law_distribution(exponent: float,
# Grains radii (um) at which the size distribution is sampled
radii = (r_bins[1:]+r_bins[:-1])/2.

# Number of grains per radius bin
# Number of grains per radius bins size
dn_dr = radii**exponent

# Total of grains to one
n_grains = np.sum(r_width*dn_dr)

# Normalize the size distribution to 1 grain
dn_dr /= n_grains
dn_dr /= np.sum(r_width*dn_dr)

# Number of grains for each radius bin
dn_grains = dn_dr*r_width

return dn_dr, r_width, radii
return dn_grains, r_width, radii


@typechecked
def dust_cross_section(dn_dr: np.ndarray,
r_width: np.ndarray,
def dust_cross_section(dn_grains: np.ndarray,
radii: np.ndarray,
wavelength: float,
n_index: float,
Expand All @@ -191,10 +168,8 @@ def dust_cross_section(dn_dr: np.ndarray,
Parameters
----------
dn_dr : np.ndarray
Number of grains per radius bin, normalized to an integrated value of 1 grain.
r_width : np.ndarray
Widths of the radius bins (um).
dn_grains : np.ndarray
Number of grains in each radius bin, normalized to a total of 1 grain.
radii : np.ndarray
Grain radii (um).
wavelength : float
Expand All @@ -213,8 +188,6 @@ def dust_cross_section(dn_dr: np.ndarray,
c_ext = 0.

for i, item in enumerate(radii):
# mean_radius = (r_lognorm[i+1]+r_lognorm[i]) / 2. # (um)

# From the PyMieScatt documentation: When using PyMieScatt, pay close attention to
# the units of the your inputs and outputs. Wavelength and particle diameters are
# always in nanometers, efficiencies are unitless, cross-sections are in nm2,
Expand All @@ -226,7 +199,7 @@ def dust_cross_section(dn_dr: np.ndarray,
asCrossSection=False)

if 'Qext' in mie:
c_ext += np.pi*item**2*mie['Qext']*dn_dr[i]*r_width[i] # (um2)
c_ext += np.pi*item**2*mie['Qext']*dn_grains[i] # (um2)

else:
raise ValueError('Qext not found in the PyMieScatt dictionary.')
Expand Down Expand Up @@ -272,7 +245,7 @@ def calc_reddening(filters_color: Tuple[str, str],

filters = [extinction[0], filters_color[0], filters_color[1]]

dn_dr, r_width, radii = log_normal_distribution(radius_g, 2., 100)
dn_grains, _, radii = log_normal_distribution(radius_g, 2., 100)

c_ext = {}

Expand All @@ -289,16 +262,14 @@ def calc_reddening(filters_color: Tuple[str, str],
# Average cross section of the three axes

if i == 0:
c_ext[item] = dust_cross_section(dn_dr,
r_width,
c_ext[item] = dust_cross_section(dn_grains,
radii,
data[wavel_index, 0],
data[wavel_index, 1],
data[wavel_index, 2]) / 3.

else:
c_ext[item] += dust_cross_section(dn_dr,
r_width,
c_ext[item] += dust_cross_section(dn_grains,
radii,
data[wavel_index, 0],
data[wavel_index, 1],
Expand All @@ -316,8 +287,7 @@ def calc_reddening(filters_color: Tuple[str, str],

wavel_index = (np.abs(data[:, 0] - filter_wavel)).argmin()

c_ext[item] += dust_cross_section(dn_dr,
r_width,
c_ext[item] += dust_cross_section(dn_grains,
radii,
data[wavel_index, 0],
data[wavel_index, 1],
Expand Down

0 comments on commit 5da6aa0

Please sign in to comment.