The PhotoDissociation Region Toolbox
Marc Pound & Mark Wolfire
University of Maryland
University of Maryland
A python toolkit for analyzing observations of PDRs from JWST, SOFIA, ALMA, Herschel, and many more¶
Thanks to our sponsors:
- NASA ADAP #80NSSC19K0573 https://dustem.astro.umd.edu
- SOFIA FEEDBACK Legacy Program SOF070077 https://feedback.astro.umd.edu
- JWST-ERS PDRs4All program ID 1288 https://pdrs4all.org
Reliable Astrophysics at Everyday Low, Low Prices! ®
JWST Webbinar 12/06/2022
Getting the PDR Toolbox¶
pip install pdrtpy
Download the example notebooks¶
git clone https://github.com/mpound/pdrtpy-nb.git
This notebook is in notebooks/PDRToolboxDemo.ipynb
To Learn More¶
Visit our website https://dustem.astro.umd.edu
Full documentation at https://pdrtpy.readthedocs.io
Our AJ paper in press: https://arxiv.org/abs/2210.08062
In [1]:
# all the imports needed for this demo
from pdrtpy.measurement import Measurement
from pdrtpy.modelset import ModelSet
from pdrtpy.tool.lineratiofit import LineRatioFit
from pdrtpy.tool.h2excitation import H2ExcitationFit
from pdrtpy.plot.modelplot import ModelPlot
from pdrtpy.plot.excitationplot import ExcitationPlot
from pdrtpy.plot.lineratioplot import LineRatioPlot
import pdrtpy.pdrutils as utils
from astropy.nddata import StdDevUncertainty
import astropy.units as u
from astropy.table import Table
import numpy as np
/home/mpound/src/pdrtpy/pdrtpy/pbar.py:10: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console) import tqdm.autonotebook as tqdm
Example 1: Finding density and radiation field¶
This example shows use the PDRT Toolbox to determine the PDR radiation field $G_0$ and hydrogen nucleus volume density $n$ from your spectral line and far-infrared (FIR) continuum data into the PDR Toolbox. The case is for single pointing observations.
Step 1. Input your intensity observations ("Measurements")¶
In [2]:
myunit = "erg s-1 cm-2 sr-1" # default unit for value and error
m1 = Measurement(data=3.6E-4,uncertainty = StdDevUncertainty(1.2E-4),
identifier="OI_63",unit=myunit)
m2 = Measurement(data=1E-6,uncertainty = StdDevUncertainty([3E-7]),
identifier="CI_609",unit=myunit)
m3 = Measurement(data=26,uncertainty = StdDevUncertainty([5]),
identifier="CO_43",restfreq="461.04077 GHz", unit="K km/s")
m4 = Measurement(data=8E-5,uncertainty = StdDevUncertainty([8E-6]),
identifier="CII_158",unit=myunit)
a = [m1,m2,m3,m4]
Step 2. Choose a set of models¶
In this case we choose Wolfire-Kaufman 2020 set of models with metallicity z=1. The models in this ModelSet are listed. We also have KOSMA-Tau models available.
The models are stored as intensity ratio maps and cover [C I], [C II], [O I], $H_2$, CO and $^{13}$CO up to J=20, and more. The PDR Toolbox calculates appropriate ratios from your intensity observations, propogating errors.
In [3]:
ms = ModelSet("wk2020",z=1)
ms.supported_ratios.show_in_notebook()
Out[3]:
Table length=73
| idx | title | ratio label |
|---|---|---|
| 0 | [O I] 63 $\mu$m / [C II] 158 $\mu$m | OI_63/CII_158 |
| 1 | ([O I] 63 $\mu$m + [C II] 158 $\mu$m) / I$_{FIR}$ | OI_63+CII_158/FIR |
| 2 | ([O I] 145 $\mu$m + [C II] 158 $\mu$m) / I$_{FIR}$ | OI_145+CII_158/FIR |
| 3 | [O I] 145 $\mu$m / [O I] 63 $\mu$m | OI_145/OI_63 |
| 4 | [C I] 370 $\mu$m / [C I] 609 $\mu$m | CI_370/CI_609 |
| 5 | [C II] 158 $\mu$m / [C I] 609 $\mu$m | CII_158/CI_609 |
| 6 | [C II] 158 $\mu$m / [O I] 145 $\mu$m | CII_158/OI_145 |
| 7 | [C II] 158 $\mu$m / I$_{FIR}$ | CII_158/FIR |
| 8 | [C II] 158 $\mu$m / CO(J=1-0) | CII_158/CO_10 |
| 9 | [C II] 158 $\mu$m / CO(J=2-1) | CII_158/CO_21 |
| 10 | [C II] 158 $\mu$m / CO(J=3-2) | CII_158/CO_32 |
| 11 | [C II] 158 $\mu$m / CO(J=6-5) | CII_158/CO_65 |
| 12 | [C II] 158 $\mu$m / $^{13}$CO(J=1-0) | CII_158/13CO_10 |
| 13 | CO(J=2-1) / CO(J=1-0) | CO_21/CO_10 |
| 14 | CO(J=3-2) / CO(J=1-0) | CO_32/CO_10 |
| 15 | CO(J=3-2) / CO(J=2-1) | CO_32/CO_21 |
| 16 | [C I] 609 $\mu$m / CO(J=4-3) | CI_609/CO_43 |
| 17 | [C I] 609 $\mu$m / CO(J=2-1) | CI_609/CO_21 |
| 18 | [C I] 609 $\mu$m / CO(J=1-0) | CI_609/CO_10 |
| 19 | [C I] 609 $\mu$m / $^{13}$CO(J=1-0) | CI_609/13CO_10 |
| 20 | [C I] 609 $\mu$m / $^{13}$CO(J=2-1) | CI_609/13CO_21 |
| 21 | [C I] 609 $\mu$m / $^{13}$CO(J=3-2) | CI_609/13CO_32 |
| 22 | CO(J=4-3) / CO(J=2-1) | CO_43/CO_21 |
| 23 | CO(J=6-5) / CO(J=1-0) | CO_65/CO_10 |
| 24 | CO(J=7-6) / CO(J=1-0) | CO_76/CO_10 |
| 25 | CO(J=7-6) / CO(J=2-1) | CO_76/CO_21 |
| 26 | CO(J=7-6) / CO(J=4-3) | CO_76/CO_43 |
| 27 | CO(J=7-6) / CO(J=5-4) | CO_76/CO_54 |
| 28 | CO(J=7-6) / CO(J=6-5) | CO_76/CO_65 |
| 29 | CO(J=8-7) / CO(J=5-4) | CO_87/CO_54 |
| 30 | CO(J=8-7) / CO(J=6-5) | CO_87/CO_65 |
| 31 | CO(J=9-8) / CO(J=5-4) | CO_98/CO_54 |
| 32 | CO(J=9-8) / CO(J=6-5) | CO_98/CO_65 |
| 33 | CO(J=10-9) / CO(J=5-4) | CO_109/CO_54 |
| 34 | CO(J=10-9) / CO(J=6-5) | CO_109/CO_65 |
| 35 | CO(J=10-9) / CO(J=7-6) | CO_109/CO_76 |
| 36 | CO(J=11-10) / CO(J=5-4) | CO_1110/CO_54 |
| 37 | CO(J=11-10) / CO(J=6-5) | CO_1110/CO_65 |
| 38 | CO(J=12-11) / CO(J=5-4) | CO_1211/CO_54 |
| 39 | CO(J=12-11) / CO(J=6-5) | CO_1211/CO_65 |
| 40 | CO(J=13-12) / CO(J=5-4) | CO_1312/CO_54 |
| 41 | CO(J=13-12) / CO(J=6-5) | CO_1312/CO_65 |
| 42 | CO(J=14-13) / CO(J=5-4) | CO_1413/CO_54 |
| 43 | CO(J=14-13) / CO(J=6-5) | CO_1413/CO_65 |
| 44 | CO(J=15-14) / CO(J=5-4) | CO_1514/CO_54 |
| 45 | CO(J=15-14) / CO(J=6-5) | CO_1514/CO_65 |
| 46 | CO(J=16-15) / CO(J=5-4) | CO_1615/CO_54 |
| 47 | CO(J=16-15) / CO(J=6-5) | CO_1615/CO_65 |
| 48 | CO(J=17-16) / CO(J=5-4) | CO_1716/CO_54 |
| 49 | CO(J=17-16) / CO(J=6-5) | CO_1716/CO_65 |
| 50 | CO(J=18-17) / CO(J=5-4) | CO_1817/CO_54 |
| 51 | CO(J=18-17) / CO(J=6-5) | CO_1817/CO_65 |
| 52 | $^{13}$CO(J=1-0) / CO(J=1-0) | 13CO_10/CO_10 |
| 53 | $^{13}$CO(J=2-1) / CO(J=2-1) | 13CO_21/CO_21 |
| 54 | $^{13}$CO(J=3-2) / CO(J=2-1) | 13CO_32/CO_21 |
| 55 | $^{13}$CO(J=4-3) / CO(J=2-1) | 13CO_43/CO_21 |
| 56 | $^{13}$CO(J=5-4) / CO(J=2-1) | 13CO_54/CO_21 |
| 57 | $^{13}$CO(J=6-5) / CO(J=2-1) | 13CO_65/CO_21 |
| 58 | $^{13}$CO(J=7-6) / CO(J=2-1) | 13CO_76/CO_21 |
| 59 | $^{13}$CO(J=3-2) / $^{13}$CO(J=1-0) | 13CO_32/13CO_10 |
| 60 | H$_2$0-0S(2) 12.3 $\mu$m / H$_2$0-0S(0) 28.2 $\mu$m | H200S2/H200S0 |
| 61 | H$_2$0-0S(2) 12.3 $\mu$m / H$_2$0-0S(1) 17 $\mu$m | H200S2/H200S1 |
| 62 | H$_2$0-0S(3) 9.7 $\mu$m / H$_2$0-0S(0) 28.2 $\mu$m | H200S3/H200S0 |
| 63 | H$_2$0-0S(3) 9.7 $\mu$m / H$_2$0-0S(1) 17 $\mu$m | H200S3/H200S1 |
| 64 | H$_2$0-0S(3) 9.7 $\mu$m / H$_2$0-0S(2) 12.3 $\mu$m | H200S3/H200S2 |
| 65 | [Fe II] 1.60 $\mu$m/[Fe II] 1.64 $\mu$m | FEII_1.60/FEII_1.64 |
| 66 | [Fe II] 1.64 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_1.64/FEII_5.34 |
| 67 | [Fe II] 5.67 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_5.67/FEII_5.34 |
| 68 | [Fe II] 17.94 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_17.94/FEII_5.34 |
| 69 | [Fe II] 22.90 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_22.90/FEII_5.34 |
| 70 | [Fe II] 25.99 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_26/FEII_5.34 |
| 71 | [Ar III] 21.83 $\mu$m/[Ar III] 8.99 $\mu$m | ARIII_21.83/ARIII_8.99 |
| 72 | [Ar V] 13.07 $\mu$m/[Ar V] 7.91 $\mu$m | ARV_13.07/ARV_7.91 |
Step 3. Create the fitter and run it.¶
The results are stored in member variables as Measurements and in an lmfit.ModelResult object.
In [4]:
p = LineRatioFit(ms,measurements=a)
p.run()
/home/mpound/src/pdrtpy/pdrtpy/tool/lineratiofit.py:384: UserWarning: LineRatioFit: No beam parameters in Measurement headers, assuming they are all equal! self._check_compatibility()
Converting K km/s to erg / (cm2 s sr) using Factor = +1.004E-07 g / (cm K s2) fitted 1 of 1 pixels got 0 exceptions
In [5]:
p.fit_result[0]
Out[5]:
Fit Statistics
| fitting method | leastsq | |
| # function evals | 13 | |
| # data points | 3 | |
| # variables | 2 | |
| chi-square | 0.04123761 | |
| reduced chi-square | 0.04123761 | |
| Akaike info crit. | -8.86105089 | |
| Bayesian info crit. | -10.6638263 |
Variables
| name | value | standard error | relative error | initial value | min | max | vary |
|---|---|---|---|---|---|---|---|
| density | 40730.7043 | 2467.35685 | (6.06%) | 42169.65034285822 | 10.0000000 | 10000000.0 | True |
| radiation_field | 0.35817960 | 0.03952220 | (11.03%) | 0.37941973904039467 | 5.0596e-04 | 5059.64354 | True |
Correlations (unreported correlations are < 0.100)
| density | radiation_field | -0.3572 |
Step 5: Visualize the results¶
Create a plotter and feed it the LineRatioFit object. Then you can call various methods to plot the fit. Here are three examples:
- Your observed ratios with errors on their matching models in $(n,G_0)$ space
- The $\chi^2$ plot of the fit
- An overlay of all your observed ratios with errors in in $(n,G_0)$ space
In [6]:
plot = LineRatioPlot(p)
plot.ratios_on_models(yaxis_unit="Habing",image=True,norm='simple',
ncols=3, cmap='plasma',
meas_color=['red'],label=True,colorbar=True)
# Save the figure to a PNG
plot.savefig("modelfits.png")
In [7]:
plot.reduced_chisq(cmap='gray_r',norm='log',label=True,
colors='white', legend=True,vmax=8E4,
yaxis_unit='Habing',
figsize=(5,7))
In [8]:
plot.overlay_all_ratios(yaxis_unit="Habing",figsize=(5,5))
# add a marker for G0 as determined from FIR observations
plot.text(1000,240,r"$G_{0,FIR}$",fontsize='large',color='k')
plot._plt.hlines(200,10,1E7,color='k',linewidth=1)
Out[8]:
<matplotlib.collections.LineCollection at 0x7f714ab5a070>
You can also fit maps!¶
Here we use the N22 [C II], [O I], and FIR maps from Jameson et al. 2018 and a custom ModelSet for the SMC.
In [9]:
cii_meas = Measurement.read("../data/n22_cii_flux_error.fits",
identifier="CII_158")
FIR_meas = Measurement.read("../data/n22_FIR_flux_error.fits",
identifier="FIR")
oi_meas = Measurement.read("../data/n22_oi_flux_error.fits",
identifier="OI_63")
smcmod = ModelSet("smc",z=0.1)
p = LineRatioFit(modelset=smcmod,
measurements=[cii_meas,FIR_meas,oi_meas])
p.run()
WARNING: VerifyWarning: Invalid 'BLANK' keyword in header. The 'BLANK' keyword is only applicable to integer data, and will be ignored in this HDU. [astropy.io.fits.hdu.image] WARNING: UnitsWarning: 'W/m2/sr' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic] WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 51544.500000 from DATE-OBS'. [astropy.wcs.wcs]
0%| | 0/11259 [00:00<?, ?it/s]
fitted 4768 of 11259 pixels got 0 exceptions
In [10]:
plot = LineRatioPlot(p)
plot.density(contours=True,norm="log")
In [11]:
plot.radiation_field(units="Habing",contours=True,norm="simple")
Example 2: $H_2$ Excitation Diagrams¶
This example will go through the steps of computing the temperature and column density using $H_2$ line intensities plotted in an excitation diagram. We will show fitting of an excitation diagram comprised of single pointing Measurements, both with fixed ortho-to-para ratio (OPR) and allowing OPR to vary. The data are from Sheffer et al 2011.
As before, you input intensities as Measurements, instantiate a fitter, run it, and look at results.¶
In [12]:
intensity = dict()
intensity['H200S0'] = 3.003e-05
intensity['H200S1'] = 3.143e-04
intensity['H200S2'] = 3.706e-04
intensity['H200S3'] = 1.060e-03
intensity['H200S4'] = 5.282e-04
intensity['H200S5'] = 5.795e-04
a = []
for i in intensity:
m = Measurement(data=intensity[i],
uncertainty=StdDevUncertainty(0.75*intensity[i]),
identifier=i,unit="erg cm-2 s-1 sr-1")
a.append(m)
In [13]:
h = H2ExcitationFit(a)
h.run()
print('T(cold) = {:>8.2f}'.format(h.tcold))
print('T(hot) = {:>8.2f}'.format(h.thot))
print(f'N(cold) = {h.cold_colden:3.2E}')
print(f'N(hot) = {h.hot_colden:3.2E}')
print(f'N(total) = {h.total_colden:+.1e}')
fitted 1 of 1 pixels got 0 exceptions and 0 bad fits T(cold) = 123.77 +/- 88.69 K T(hot) = 633.39 +/- 59.82 K N(cold) = 5.64E+21 +/- 1.68E+22 1 / cm2 N(hot) = 2.25E+20 +/- 1.10E+20 1 / cm2 N(total) = +5.9e+21 +/- +1.7e+22 1 / cm2
In [14]:
hplot = ExcitationPlot(h,"H_2")
hplot.ex_diagram(show_fit=True,ymax=21)
The first fit looked pretty good, but we can do better by allowing the OPR to vary as well.¶
In [15]:
h.run(fit_opr=True)
hplot.ex_diagram(show_fit=True,ymax=21)
fitted 1 of 1 pixels got 0 exceptions and 0 bad fits
You can also plot multiple vibrational levels.¶
This example uses Orion Bar data from Kaplan et al 2021.
In [ ]:
# read in the Kaplan data and make Measurements from it
ktab = Table.read(utils.get_testdata("Kaplan+2021.tab"),format='ipac')
morion=[]
for d,un,i in zip(ktab["Forionbar"],ktab["e_Forionbar"],ktab["Line"]):
if not np.ma.is_masked(d):
morion.append(Measurement(d,uncertainty=StdDevUncertainty(un),
identifier=i,unit="erg s-1 cm-2 sr-1"))
In [27]:
h = H2ExcitationFit(measurements=morion)
hplot = ExcitationPlot(h,"H_2")
# Give the parameter `label='v'` to show the vibrational levels
hplot.ex_diagram(ymin=15.5,ymax=20.5,label="v",xmax=56000)
Example 3: Phase space and isotemperature plots¶
You can plot an model ratio against any other model ratio (or model intensity) and overlay observations. In these examples we use the [H II] region diagnostic lines of Fe and Ar. You create a ModelPlot with your chose ModelSet and away you go!
In [17]:
mp = ModelPlot(ms)
m1 = Measurement(data=[0.1],
uncertainty = StdDevUncertainty(0.025),
identifier="FEII_1.60/FEII_1.64",unit="")
m2 = Measurement(data=[0.5],
uncertainty = StdDevUncertainty(0.1),
identifier="FEII_1.64/FEII_5.34",unit="")
mp.phasespace(['FEII_1.60/FEII_1.64','FEII_1.64/FEII_5.34'],
nax2_clip=[10,1E6]*u.Unit("cm-3"),
nax1_clip=[2E3,8E3]*u.Unit("K"),
measurements=[m1,m2],errorbar=True)
In [24]:
nax1clip = [1E3,5E3]*u.Unit("K")
nax2clip = [10,1E6]*u.Unit("cm-3")
mp.phasespace(["ARIII_21.83/ARIII_8.99","ARV_13.07/ARV_7.91"],
nax2_clip=nax2clip,
nax1_clip=nax1clip)
In [19]:
nax1clip = [2000,8000]*u.Unit("K")
# step=2 means plot every other line in the given range
mp.isoplot("ARIII_21.83/ARIII_8.99",
plotnaxis=1,logx=True,
nax_clip=nax1clip,step=2)
In [25]:
nax1clip = [100,1E7]*u.Unit("cm-3")
mp.isoplot("ARIII_21.83/ARIII_8.99",
plotnaxis=2,nax_clip=nax1clip,logx=False)
There is much more the PDR Toolbox can do.
We encourage you to download it and see for yourself!
We encourage you to download it and see for yourself!
We also welcome development. Fork us on github! https://github.com/mpound/pdrtpy¶
Questions?¶
In [ ]: