Skip to content

Commit

Permalink
Merge pull request #34 from eblur/astropy-units
Browse files Browse the repository at this point in the history 
Incorporate astropy units into the flow
  • Loading branch information
@eblur
eblur authored Nov 22, 2022
2 parents a4ba004 + 6c2f431 commit cd5a2ea
Show file tree
Hide file tree
Showing 38 changed files with 3,165 additions and 1,559 deletions.
301 changes: 240 additions & 61 deletions examples/plot_cmindex.ipynb
View file

Large diffs are not rendered by default.

205 changes: 161 additions & 44 deletions examples/plot_graindist.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,15 +9,22 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import astropy.units as u\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook describes how to set up a grain size distribution: `newdust.graindist.GrainDist`. Grain size distributions combine the size distribution (number density distribution of dust grains with different radii, $a$) with the material properties of the grain (such as `rho`, provided by `newdust.graindist.composition`. It does not matter what the composition of the dust grains are, except for the material density of those grains."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"from newdust import graindist"
Expand All @@ -29,8 +36,7 @@
"metadata": {},
"outputs": [],
"source": [
"# Some basic magic numbers\n",
"\n",
"# Some generic parameters\n",
"AMAX = 0.3 # maximum grain size, in um\n",
"RHO = 3.0 # grain material density, in g cm^-3\n",
"MD = 1.e-5 # dust mass column, g cm^-2"
Expand All @@ -40,19 +46,40 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Making a GrainDist object from scratch"
"## Making a GrainDist object from scratch\n",
"\n",
"There are a few shortcuts for setting up grain distributions (which connect a size distribution function in the `newdust.graindist.sizedist` module with the material properties of the grain, such as density, in th `newdust.graindist.composition` module). Here, we will show how to set them up together by hand.\n",
"\n",
"First, one must choose a function for the grain size distribution. There are three built-in `newdust.sizedist` classes:\n",
"\n",
"* **Grain** (`newdust.graindist.sizedist.Grain`): A single grain size distribution consisting solely of grains with radius `a`\n",
"* **Powerlaw** (`newdust.graindist.sizedist.Powerlaw`): A power law distribution of grains with radii between `amin` and `amax` and a power law slope of negative `p`. Use the keyword `na` to set the number of grid points to use for the grain radii, `a`.\n",
"* **ExpCutoff** (`newdust.graindist.sizedist.ExpCutoff`): A power law with an exponential cutoff, which describes a continuous distribution of grains with radii starting with `amin` and following a power law slope of negative `p`. The radius `acut` is used to taper the grain size distribution on the large end of the distribution. This grain size distibution follows the formula:\n",
"$$ \\frac{dn}{da} \\propto a^{-p} \\exp(-a/a_{\\rm cut})$$\n",
"The maximum grain size evaluated for the distribution is `nfold` $\\times$ `acut`\n",
"\n",
"Continuous distributions (`Powerlaw` and `ExpCutoff`) have units of `mdens` in g cm$^{-2}$ $\\mu$m$^{-1}$ and `ndens` in cm$^{-2}$ $\\mu$m$^{-1}$. The single grain size distribution, `Grain` has units of `mdens` in g cm$^{-2}$ and `ndens` in cm$^{-2}$.\n",
"\n",
"The input parameter `md` provides the total integrated dust mass column in g cm$^{-2}$. The grain size distributions are automatically normalized so that `mdens` will integrate over the grain radius values to return `md`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# To completely customize your grain size distribution, \n",
"# you must first set up a newdust.sizedist object\n",
"\n",
"# This only provides the functional form of the grain size distribution\n",
"# It does not connect to the composition or density of the dust grains\n",
"MRN = graindist.sizedist.Powerlaw(amax=AMAX)\n",
"ECUT = graindist.sizedist.ExpCutoff(acut=AMAX)\n",
"\n",
"# Next, you much specify a composition for the grains\n",
"# Composition provides the material density, `rho`, which is needed to translate\n",
"# dust mass column into number density of dust grains\n",
"SIL = graindist.composition.CmSilicate(rho=RHO)"
]
},
Expand All @@ -62,8 +89,9 @@
"metadata": {},
"outputs": [],
"source": [
"gd_mrn = graindist.GrainDist(MRN, SIL, md=MD, custom=True)\n",
"gd_ecut = graindist.GrainDist(ECUT, SIL, md=MD, custom=True)"
"# Here we combine the sizedist and composition objects to define a GrainDist object\n",
"gd_mrn = graindist.GrainDist(MRN, SIL, md=MD, custom=True) # typical powerlaw (MRN)\n",
"gd_ecut = graindist.GrainDist(ECUT, SIL, md=MD, custom=True) # MRN with an exponential cut-off"
]
},
{
Expand All @@ -72,29 +100,74 @@
"metadata": {},
"outputs": [],
"source": [
"# Use the built-in plotting function to visualize the results\n",
"# It is conventional to plot the grain size distribution as (dn/da) a^4\n",
"# Particularly useful for X-ray studies, because the scattering cross-section is propto a^4\n",
"ax = plt.subplot(111)\n",
"gd_mrn.plot(ax, color='k', lw=2, label='Powerlaw')\n",
"gd_ecut.plot(ax, color='b', lw=2, alpha=0.8, label='Exponential cut-off')\n",
"plt.legend(loc='upper left', frameon=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Or, you could plot the properties by hand to get rid of the a^4 component\n",
"plt.plot(gd_mrn.a, gd_mrn.ndens, 'k-', lw=2, label='Powerlaw')\n",
"plt.plot(gd_ecut.a, gd_ecut.ndens, 'b-', lw=2, alpha=0.8, label='Exponential cut-off')\n",
"plt.xlabel(gd_mrn.a.unit)\n",
"plt.ylabel(r'$dn/da$ (cm$^{-2}$ um$^{-1}$)')\n",
"plt.loglog()\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"## Note that the grain radii have units on them... \n",
"print(gd_mrn.a)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ... but the number densities do not\n",
"print(gd_mrn.ndens)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Making a GrainDist object with helper function\n",
"\n",
"This is a shortcut for making some common grain distributions\n",
"Instead of supplying the `newdust.graindist.sizedist` and `newdust.graindist.composition` objects directly, you can provide strings. \n",
"\n",
"You can change AMAX and RHO from here, too\n",
"Options for the size distribution function:\n",
"+ \"Grain\"\n",
"+ \"Powerlaw\"\n",
"+ \"ExpCutoff\"\n",
"\n",
"For the different grain size distributions, amax acts in a different way:\n",
"Options for the composition:\n",
"+ \"Silicate\" - uses Draine (2003) silicate properties\n",
"+ \"Graphite\" - uses Draine (2003) graphite properties, perpendicular orientation\n",
"+ \"Drude\" - uses the Drude approximation\n",
"\n",
"+ *graindist.sizedist.Grain:* amax sets the singular grain size\n",
"You can also use the keywords `amax` and `rho` to customize. The `rho` keyword value is used to initiate the composition objct. The `amax` keyword value is used to initiate the size distribution function, and acts differently depending on the size distribution you picked:\n",
"+ Grain: amax sets the singular grain size\n",
"+ Powerlaw: amax sets the maximum grain size in the distribution\n",
"+ ExpCutoff: amax sets the `acut` value\n",
"\n",
"+ *graindist.sizedist.Powerlaw:* amax sets the maximum grain size in the distribution\n",
"\n",
"+ *graindist.sizedist.ExpCutoff:* amax sets the \"acut\" value"
"Any other keywords provided to the `newdust.graindist.GrainDist` are used to initiate the size distribution."
]
},
{
Expand All @@ -103,6 +176,8 @@
"metadata": {},
"outputs": [],
"source": [
"# To use the shortcuts, provide the name of the sizedist class you want to use, as a string\n",
"# You can also provide the grain composition as a string: 'Silicate', 'Graphite', or 'Drude'\n",
"gd_mrn2 = graindist.GrainDist('Powerlaw', 'Silicate', amax=AMAX, rho=RHO, md=MD)\n",
"gd_ecut2 = graindist.GrainDist('ExpCutoff', 'Silicate', amax=AMAX, rho=RHO, md=MD)"
]
Expand All @@ -114,31 +189,34 @@
"outputs": [],
"source": [
"ax = plt.subplot(111)\n",
"\n",
"# We still get the same power law as before\n",
"gd_mrn.plot(ax, color='k', lw=2, label='')\n",
"gd_mrn2.plot(ax, color='r', lw=2, ls='--', label='Powerlaw')\n",
"\n",
"gd_ecut.plot(ax, color='k', lw=2, label='')\n",
"# And we get the same ExpCutoff distribution as before\n",
"gd_ecut.plot(ax, color='k', alpha=0.5, lw=2, label='')\n",
"gd_ecut2.plot(ax, color='b', lw=2, ls='--', label='Exponential Cutoff')\n",
"\n",
"plt.legend(loc='upper left', frameon=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Note that the values I've used in this example are different from the defaults\n",
"### Note that the values used in this example are different from the defaults\n",
"\n",
"Silicate has a default grain material density of 3.8 g cm^-3"
"Silicate has a default grain material density of 3.8 g cm$^{-3}$. If you don't define `rho`, there is a default value for each."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Set up the same size distributions, but fall back on the default silicate density\n",
"gd_mrn3 = graindist.GrainDist('Powerlaw', 'Silicate', md=MD)\n",
"gd_ecut3 = graindist.GrainDist('ExpCutoff', 'Silicate', md=MD)"
]
Expand All @@ -149,22 +227,31 @@
"metadata": {},
"outputs": [],
"source": [
"# Plot to compare the results of using a slightly larger `rho` (the default)\n",
"ax = plt.subplot(111)\n",
"gd_mrn.plot(ax, color='k', lw=2, label='')\n",
"gd_mrn3.plot(ax, color='r', lw=2, ls='--', label='Powerlaw')\n",
"\n",
"gd_ecut.plot(ax, color='k', lw=2, label='')\n",
"gd_ecut3.plot(ax, color='b', lw=2, ls='--', label='Exponential Cutoff')\n",
"plt.legend(loc='upper left', frameon=False)"
"plt.legend(loc='upper left', frameon=False)\n",
"\n",
"# The default grain density is larger than 3.0 g cm^-3, \n",
"# so the number density with the default `rho` is smaller\n",
"# (fewer grains needed to obtain the same dust mass)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## You can also get the mass density of dust grains\n",
"### Hybride method for initiating a Grain Distribution\n",
"\n",
"As described above, most of the keyword arguments provided to `newdust.graindist.GrainDist` are passed to the size distribution function. In the example below, I set up a power law distribution that extends to larger dust grain sizes than the default MRN function. I set the amax value to a number that is two times larger than the maximum value in `gd_mrn`.\n",
"\n",
"But the GrainDist.plot function only does number density"
"I also change the number of grain radii used to sample the distribution using the keyword `na` (i.e., the number of grain radius $a$ values used to evaluate the size distribution function). I also make it so that the size distribution is uniformly sampled over log-space (instead of linear) by setting the keyword `log=True`.\n",
"\n",
"The plot below shows how it compares to the original `gd_mrn` grain distribution. I plot each sample point in the size distribution with an 'o', to demonstrate the difference in sampling method."
]
},
{
Expand All @@ -173,8 +260,41 @@
"metadata": {},
"outputs": [],
"source": [
"plt.plot(gd_mrn.a, gd_mrn.mdens, label='Powerlaw')\n",
"plt.plot(gd_ecut.a, gd_ecut.mdens, label='ExpCutoff')\n",
"gd_powerlaw2 = graindist.GrainDist('Powerlaw', 'Silicate', \n",
" md=MD, rho=RHO, amax=2.0*AMAX, na=10, log=True)\n",
"\n",
"ax = plt.subplot(111)\n",
"gd_mrn.plot(ax, marker='o', color='k', label='MRN')\n",
"gd_powerlaw2.plot(ax, marker='o', color='b', label='Powerlaw2')\n",
"plt.loglog()\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the \"Powerlaw2\" distribution follows the same slope, but extends to largr radii, and overall there are fewer dust grains of each radius. That's because both distributions integrate to the same total mass."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Look at the mass density distribution of dust grains\n",
"\n",
"The `md` attribute contains the total dust mass, but the `mdens` attribute (for continuous distributions like `Powerlaw` and `ExpCutoff`) has units of g cm$^{-2}$ um$^{-1}$ and will integrate to the total mass column."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot mdens to show that it is a continuous distribution\n",
"plt.plot(gd_mrn.a, gd_mrn.mdens, 'k-', label='Powerlaw')\n",
"plt.plot(gd_ecut.a, gd_ecut.mdens, 'b-', label='ExpCutoff')\n",
"plt.loglog()\n",
"plt.xlabel('Radius (um)')\n",
"plt.ylabel('Mass density (g cm$^{-2}$ um$^{-1}$)')\n",
Expand All @@ -185,46 +305,43 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### You can't plot single grain sizes\n",
"### Single grain size plots are un-informative\n",
"\n",
"It prints the number density instead"
"You can plot them, but it's just a single point! Also nthe units displayed in the plot are technically wrong. The `ndens` units for a `newdust.graindist.sizedist.Grain` object is number density ($cm^{-2}$) while continuous size distributions have units of $cm^{-2}$ $\\mu$m$^{-1}$, because they are meant to be integrated."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"gd_single = graindist.GrainDist('Grain', 'Silicate', amax=AMAX)"
"gd_single = graindist.GrainDist('Grain', 'Silicate', amax=AMAX)\n",
"\n",
"ax = plt.subplot(111)\n",
"gd_single.plot(ax, marker='o')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gd_single.plot()"
]
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "newdust-dev",
"language": "python",
"name": "python3"
"name": "newdust-dev"
},
"language_info": {
"codemirror_mode": {
Expand All @@ -236,7 +353,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.2"
"version": "3.10.1"
}
},
"nbformat": 4,
Expand Down
Loading

0 comments on commit cd5a2ea

Please sign in to comment.