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")