University of Maryland
A python toolkit for analyzing observations of PDRs from SOFIA , JWST, ALMA, APEX, Herschel, and many more ¶
- 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
SOFIA School 2023
What does the PDR Toolbox do? ¶
The PDR Toolbox lets you analyze your observations of photodissociation regions by
-
Fitting spectral line and continuum data to state-of-the art PDR models to deterimine radiation field and $H$ density
- Both Wolfire-Kaufman and KOSMA-tau models are supported
- Plotting data on models in various ways
- Fitting $H_2$ intensities on excitation diagrams to determine temperature, column density, and ortho-to-para ratio
- Determining electron temperature and density in ionized gas regions from observations
How to get the PDR Toolbox
Install
pdrtpy
(or upgrade if you already have it)
¶
pip install pdrtpy -U
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 about PDRT: Pound & Wolfire 2023, AJ, 165, 25 : https://ui.adsabs.harvard.edu/abs/2023AJ....165...25P/abstract
Our AAS Youtube Author video! https://www.youtube.com/watch?v=tYomntouCT8&list=PLFhVT3VzlwKrArGdcNhtSCF4NbDBd6K5_&index=14
# 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
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") ¶
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 example, 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.
All models are stored as intensity or 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.
ms = ModelSet("wk2020",z=1)
ms.supported_ratios.show_in_notebook()
| 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 | H$_2$2-1S(1) 2.24 $\mu$m / H$_2$1-0S(1) 2.12 $\mu$m | H221S1/H210S1 |
| 66 | [Fe II] 1.60 $\mu$m/[Fe II] 1.64 $\mu$m | FEII_1.60/FEII_1.64 |
| 67 | [Fe II] 1.64 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_1.64/FEII_5.34 |
| 68 | [Fe II] 5.67 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_5.67/FEII_5.34 |
| 69 | [Fe II] 17.94 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_17.94/FEII_5.34 |
| 70 | [Fe II] 22.90 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_22.90/FEII_5.34 |
| 71 | [Fe II] 25.99 $\mu$m/[Fe II] 5.34 $\mu$m | FEII_26/FEII_5.34 |
| 72 | [Ar III] 21.83 $\mu$m/[Ar III] 8.99 $\mu$m | ARIII_21.83/ARIII_8.99 |
| 73 | [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. ¶
This example uses the default least-squares fitting; MCMC is also available. The results are stored in member variables as
Measurements
and in an
lmfit.ModelResult
object.
p = LineRatioFit(ms,measurements=a)
p.run()
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
p.fit_result[0]
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.21840 | (6.06%) | 42169.65034285822 | 10.0000000 | 10000000.0 | True |
| radiation_field | 0.35817960 | 0.03950398 | (11.03%) | 0.37941973904039467 | 5.0596e-04 | 5059.64354 | True |
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 reduced $\chi^2$ plot of the fit as a function of $(n,G_0)$
- An overlay of all your observed ratios with errors in the model $(n,G_0)$ space
Note we choose to plot radiation field in Habing units, but common alternatives Draine, Mathis, and cgs units are supported.
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")
plot.reduced_chisq(cmap='gray_r',norm='log',label=True,
colors='white', legend=True,vmax=8E4,
yaxis_unit='Habing',
figsize=(5,7))
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)
You can also fit maps! ¶
Here we use [C II], [O I], and FIR maps of the SMC molecular cloud N22 from
Jameson et al. 2018
and a custom
ModelSet
for the SMC.
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()
fitted 4768 of 11259 pixels got 0 exceptions
plot = LineRatioPlot(p)
plot.density(contours=True,norm="log",figsize=(5,8))
plot.radiation_field(units="Habing",contours=True,
norm="simple",figsize=(5,8))
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.
¶
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)
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.62 K T(hot) = 633.39 +/- 59.82 K N(cold) = 5.64E+21 +/- 1.67E+22 1 / cm2 N(hot) = 2.25E+20 +/- 1.10E+20 1 / cm2 N(total) = +5.9e+21 +/- +1.7e+22 1 / cm2
hplot = ExcitationPlot(h,"H_2")
hplot.ex_diagram(show_fit=True,ymax=21)
The first fit looked ok, but we can do better by allowing the OPR to vary as well. ¶
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.
# 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"))
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)
You can also run the fit on every pixel in a map ¶
Here is an example with maps of 6 $H_2$ lines. The maps are synthetic data made up of two Gaussian "blobs", with differing hot and cold temperatures and OPRs.
intensity = ['H200S0', 'H200S1','H200S2','H200S3','H200S4','H200S5']
meas = []
for j in intensity:
file = f"../data/{j}two_blobs_measurement.fits"
meas.append(Measurement.read(file,identifier=j))
hmap = H2ExcitationFit(meas)
hmap.run(fit_opr=True)
100%|████████████████████████████████████████████████████████████████████████| 10000/10000 [00:54<00:00, 183.95it/s]
fitted 10000 of 10000 pixels got 0 exceptions and 0 bad fits
You can examine the results with the
explore
method
¶
explore() plots data on the left panel and the fit at the selected pixel on the right panel. Click on anywhere in the left panel to see the fit for that pixel. This example uses the total column density for the left panel.
hplot = ExcitationPlot(hmap,"H_2")
hplot.explore(hmap.colden("total"),figsize=(10,6),fmt='c+')
Explore using position: (50, 50) size=1
Example 3: Phase space and isocontour plots ¶
As phase space plot shows any model ratio against any other model ratio (or model intensity) on which you can overlay observations. In these examples we use the [H II] region diagnostic lines of Fe and Ar. You create a
ModelPlot
with your chosen
ModelSet
and away you go!
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'],
nax1_clip=[2E3,8E3]*u.Unit("K"),
nax2_clip=[10,1E6]*u.Unit("cm-3"),
measurements=[m1,m2],errorbar=True)
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"],
nax1_clip=nax1clip,
nax2_clip=nax2clip)
You can make isocontour plots of individual model parameters ¶
Related to the phase space plot is the
isoplot
, which will plot one of the model axes versus model value (intensity or intensity ratio) and draw on it iso-lines of the other model axis.
In the example below, we plot lines of constant electron gas temperature $T_e$ (isotherms) as a function of electron gas density $n_e$ and [Ar III] 21.83 $\mu$m / [Ar III] 8.99 $\mu$m intensity ratio. In our ionized gas models, $T_e$ is on FITS
NAXIS1
and $n_e$ is on
NAXIS2
.
nax1clip = [2000,8000]*u.Unit("K")
# step=N means plot every N-th line in the given range.
# This helps avoid crowding on the plot.
mp.isoplot("ARIII_21.83/ARIII_8.99",
plotnaxis=1,logx=True,
nax_clip=nax1clip,step=2)
# Plotting the other model axis is easy,
# just change plotnaxis and the clip values and units
nax1clip = [100,1E7]*u.Unit("cm-3")
mp.isoplot("ARIII_21.83/ARIII_8.99",
plotnaxis=2,nax_clip=nax1clip,logx=False)
What's next for the PDR Toolbox? ¶
- additional metallicities in PDR models: z=0.03 to 5
-
expand excitation fitting
- other molecules, e.g., CO
- vibrational level fitting to compare $T_{vib}$ to $T_{rot}$, which tells you if UV-pumping of $H_2$ is important
- additional viewing angles
- deuterium chemistry (HD line)
- more robust fitting for image data
- deal with [CII] and [OI] self-absorption
- and more!
We encourage you to download it and see for yourself!