In [1]:
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

In [2]:
from astrodust import distlib

Examples of creating a power law distribution

The distlib.MRN_dist() functions return a distlib.DustSpectrum object.

The object contains an array of grain radii (a), the number density (nd), and total dust mass (md).

In [3]:
mrn_test1 = distlib.MRN_dist(0.005, 0.3, 3.5)
mrn_test2 = distlib.MRN_dist(0.005, 0.25, 3.5)
mrn_test3 = distlib.MRN_dist(0.005, 0.3, 4.0)
mrn_test4 = distlib.MRN_dist(0.005, 0.3, 3.5, na=10, log=True)

Here's a quick way to see all the keys in the DustSpectrum object.

In [4]:

<class 'astrodust.distlib.sizedist.DustSpectrum'>
['a', 'md', 'nd', 'rho']

In [5]:
plt.plot(mrn_test1.a, mrn_test1.nd, label='amax=%.3f, p=%.2f' % (mrn_test1.a[-1], 3.5))
plt.plot(mrn_test4.a, mrn_test4.nd, 'ko', label='')
plt.plot(mrn_test2.a, mrn_test2.nd, label='amax=%.3f, p=%.2f' % (mrn_test2.a[-1], 3.5))
plt.plot(mrn_test3.a, mrn_test3.nd, label='amax=%.3f, p=%.2f' % (mrn_test3.a[-1], 4.0))
plt.legend(loc='upper right', frameon=False)
plt.xlabel('grain radius (micron)')
plt.ylabel(r'$dn/da$ (cm$^{-2}$ micron$^{-1}$)')

(0.005, 0.3)

Play with WD01 dust distributions

Dust grain size distributions from Weingartner & Draine (2001)

Use the function make_WD01_DustSpectrum to return a DustSpectrum object.

In [6]:
wd_MW_gra = distlib.make_WD01_DustSpectrum(type='Graphite', verbose=False)
wd_MW_sil = distlib.make_WD01_DustSpectrum(type='Silicate', verbose=False)

wd_MW_gra_bc6 = distlib.make_WD01_DustSpectrum(type='Graphite', bc=6.0, verbose=False)
wd_MW_sil_bc6 = distlib.make_WD01_DustSpectrum(type='Silicate', bc=6.0, verbose=False)

In [7]:
plt.plot(wd_MW_gra.a, wd_MW_gra.nd*wd_MW_gra.a**4, label='Graphite MW dust')
plt.plot(wd_MW_sil.a, wd_MW_sil.nd*wd_MW_sil.a**4, label='Silicate MW dust')

plt.plot(wd_MW_gra_bc6.a, wd_MW_gra_bc6.nd*wd_MW_gra_bc6.a**4, 'b--', label='bc=6')
plt.plot(wd_MW_sil_bc6.a, wd_MW_sil_bc6.nd*wd_MW_sil_bc6.a**4, 'g--', label='bc=6')

plt.xlabel('grain radius (micron)')
plt.ylabel(r'$a^4 dn/da$ (cm$^{-2}$ um$^{3}$)')

plt.legend(loc='lower left', frameon=False)
plt.xlim(0.005, 1)

(1e-20, 1e-15)

The DustSpectrum objects also contain an integrate_dust_mass() function for calculating the total mass column of dust (g cm^-2).

In [8]:
print("Graphite dust mass = %.3e" %(wd_MW_gra.integrate_dust_mass()))
print("Silicate dust mass = %.3e" %(wd_MW_sil.integrate_dust_mass()))

Graphite dust mass = 4.931e-27
Silicate dust mass = 1.379e-26