pip install pdrtpy
git clone https://github.com/mpound/pdrtpy-nb.git
This notebook is in notebooks/PDRToolboxDemo.ipynb
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
# 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
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.
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]
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.
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 | [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 |
The results are stored in member variables as Measurements and in an lmfit.ModelResult
object.
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
p.fit_result[0]
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 |
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 |
density | radiation_field | -0.3572 |
Create a plotter and feed it the LineRatioFit object. Then you can call various methods to plot the fit. Here are three examples:
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)
<matplotlib.collections.LineCollection at 0x7f714ab5a070>
Here we use the N22 [C II], [O I], and FIR maps 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()
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
plot = LineRatioPlot(p)
plot.density(contours=True,norm="log")
plot.radiation_field(units="Habing",contours=True,norm="simple")
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.
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.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
hplot = ExcitationPlot(h,"H_2")
hplot.ex_diagram(show_fit=True,ymax=21)
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