- The aim is to find φ as a function of ionization parameter (and other parameters)
- Mainly for OB stars, but may also try RSG/AGB case
- This can then be used to study gas-grain decoupling in bow shocks
- [X] Write an input file by hand
- Start with something from test suite?
- First attempt is models/test-dust.in
- [X] Run it and see what output we need to save
export D=~/Work/CLOUDY/cloudy/source cd models $D/cloudy.exe -p test-dust
- [X] Decide which approach to use?
- Python script to write multiple input, run with
make -j 8 - Work out how to use Cloudy grid command (need to install MPI?)
- Python script to write multiple input, run with
- Do it
- Decided on approach 1
- Python script to write the Cloudy input files is in Expand template for input files below
- Using the data in the opacity files
- All sizes in microns
- 10 size bins
- 0.005 to 0.25 micron
- Factor of roughly 1.5 size increase between adjacent bins
- Bin 10, 0.21 micron, 2100 Å dominates in NIR
- graphite more important because silicate has a dip there
- Bin 04, 0.02 micron, 200 Å is what dominates in FUV
- But the contribution from many size bins is important in FUV
- see ../dust-wave-case-studies//xsec-dust-silicate_ism_10.pdf and the like
| a_min | a_max | a, μm | a, Å | |
|---|---|---|---|---|
| 01 | 0.005 | 0.00739 | 0.006 | 60. |
| 02 | 0.00739 | 0.01093 | 0.009 | 90. |
| 03 | 0.01093 | 0.01616 | 0.014 | 140. |
| 04 | 0.01616 | 0.02390 | 0.020 | 200. |
| 05 | 0.02390 | 0.03534 | 0.030 | 300. |
| 06 | 0.03534 | 0.05226 | 0.045 | 450. |
| 07 | 0.05226 | 0.07728 | 0.066 | 660. |
| 08 | 0.07728 | 0.11428 | 0.098 | 980. |
| 09 | 0.11428 | 0.16899 | 0.144 | 1440. |
| 10 | 0.16899 | 0.24990 | 0.213 | 2130. |
- The mean a is weighted by area
- These should be similar to the ones I have in my table
- Not the same!
- Problem is the
luminositycommand by default is in ionizing photons only- Need to add
totalkeyword
- Need to add
- Return to ionizing luminosities [2018-08-28 Tue]
- There are small discrepancies still in the ionizing luminosity
- Meyer (2016) Table 1 have very round numbers
- The ones in my table are similar to these, but with more sig figs
- Where did they come from?
- Possibly the Brott paper tables?
- The Cloudy models are different again
- For instance, for 10 Msun model we have S_49 = 0.00028, which is twice what I have in my table
from textwrap import dedent
import numpy as np
LSUN = 3.82e33
def star_input(id_, L4, T3, log_g):
s = f"# {id_} star" + "\n"
s += f"# T_eff = {1e3*T3:.0f} K, L = {1e4*L4:.2e} L_sun, log(g) = {log_g:.2f}" + "\n"
s += f"table star tlusty OBstar 3-dim {1e3*T3:.0f} {log_g} 0.0" + "\n"
s += f"luminosity total {np.log10(LSUN*1e4*L4):.2f}" + "\n"
return s
def hden_input(hden):
return f"hden {hden:.2f} # density of {10**hden} pcc" + "\n"
def file_stem(hden, id_):
return f"dustrad-n{int(hden):02d}-{id_}"
def intro_input(hden, id_):
s = f"title Dust radiative acceleration: star {id_}, density {10**hden} pcc" + "\n"
s += f"set save prefix \"{file_stem(hden, id_)}\"" + "\n"
return s
stars = [
# id_ L4 T3 log_g
["MS10", 0.63, 25.2, 4.2],
["MS20", 5.45, 33.9, 4.2],
["MS40", 22.2, 42.5, 4.2],
["BSG", 30.2, 23.5, 3.4],
]
def radius_input(L4, hden):
# Scale R_in to give same F_bol / N as the MS10 model with hden=1
R_in = 1e16 * np.sqrt((L4/0.63) * 10**(1.0 - hden))
s = f"# Start close in: {R_in/3.085677582e18:.5f} pc" + "\n"
s += f"radius {np.log10(R_in)}" + "\n"
return s
outro_input = dedent("""\
# Go into PDR a bit
stop temperature 4000 K linear
iterate
sphere
grains ism function sublimation
abundances HII region no grains
cosmic ray background
cmb
table ism
# save all the output we want
save overview last ".ovr"
save physical conditions last ".phys"
save continuum last ".cont"
# save continuum last every ".zcont"
save radius last ".rad"
save grain abundance last ".gabun"
save grain charge last ".gcharge"
save grain continuum last ".gcont"
save grain drift velocity last ".gdrift"
save grain potential last ".gpot"
save grain temperature last ".gtemp"
save lines, emissivity last ".emis"
H 1 6562.81A
Ca B 6562.81A
N 2 6583.45A
O 3 5006.84A
IRAC 3.60000m
IRAC 4.50000m
IRAC 5.80000m
IRAC 8.00000m
F12 12.0000m
MIPS 24.0000m
PAC1 70.0000m
PAC3 160.000m
end of lines
""")
for hden in 0.0, 1.0, 2.0, 3.0, 4.0:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden)
+ radius_input(L4, hden)
+ outro_input)
file_name = f"models/{file_stem(hden, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)- [2019-02-08 Fri] I want to compare the dust emission with the models that Kobulnicky are using, to see how important it is that they are using the wrong incident SED
- This could also explain why K17 (p 13) find higher T than predicted, because they are underestimating UV flux
- I have added a save lines emissivity command to the output section above, using the built-in bands
- [X] Unfortunately, I don’t have a working cloudy at the moment
- Now sorted, but only using gcc
export D=/Users/will/Work/CLOUDY/cloudy/source/sys_gcc
- Now sorted, but only using gcc
- [2019-04-13 Sat] I want to see how much difference this really makes at 24 micron
- Aigen Li says that 24 micron should be dominated by (stochastic heating) reported by Chip Kobulnicky
- But that contradicts what I say in sec 4.1.1 of Paper III
- Conclusions
- Even if we allow PAHs to exist in the bows, then stochastic heating makes a neglible contribution to the 24 micron emission for U > 100, which covers nearly all the observed bows. This is assuming q_PAH = 0.47%, which is the lowest in the DL07 models
- For the 8 micron emission, stochastic heating dominates unless U > 10^4. Again, assuming q_PAH = 0.47%
- Note that some bows do have U > 10^4. Most notably θ^1 D Ori has U = 5e6
- Put everything in new folder models-no-qheat/
import glob
infiles = glob.glob("models/dustrad*.in")
for infile in infiles:
with open(infile) as f:
s = f.read()
# Change the grains command to turn off quantum heating
s = s.replace(
"grains ism function sublimation",
"grains ism function sublimation no qheat",
)
with open(infile.replace("models/", "models-no-qheat/"), "w") as f:
f.write(s)ls -l models-no-qheat/*.in-rw-r--r-- 1 will staff 1105 Apr 13 18:55 models-no-qheat/dustrad-n00-BSG.in -rw-r--r-- 1 will staff 1094 Apr 13 18:55 models-no-qheat/dustrad-n00-MS10.in -rw-r--r-- 1 will staff 1107 Apr 13 18:55 models-no-qheat/dustrad-n00-MS20.in -rw-r--r-- 1 will staff 1107 Apr 13 18:55 models-no-qheat/dustrad-n00-MS40.in -rw-r--r-- 1 will staff 1107 Apr 13 18:55 models-no-qheat/dustrad-n01-BSG.in -rw-r--r-- 1 will staff 1096 Apr 13 18:55 models-no-qheat/dustrad-n01-MS10.in -rw-r--r-- 1 will staff 1109 Apr 13 18:55 models-no-qheat/dustrad-n01-MS20.in -rw-r--r-- 1 will staff 1109 Apr 13 18:55 models-no-qheat/dustrad-n01-MS40.in -rw-r--r-- 1 will staff 1109 Apr 13 18:55 models-no-qheat/dustrad-n02-BSG.in -rw-r--r-- 1 will staff 1098 Apr 13 18:55 models-no-qheat/dustrad-n02-MS10.in -rw-r--r-- 1 will staff 1112 Apr 13 18:55 models-no-qheat/dustrad-n02-MS20.in -rw-r--r-- 1 will staff 1111 Apr 13 18:55 models-no-qheat/dustrad-n02-MS40.in -rw-r--r-- 1 will staff 1111 Apr 13 18:55 models-no-qheat/dustrad-n03-BSG.in -rw-r--r-- 1 will staff 1100 Apr 13 18:55 models-no-qheat/dustrad-n03-MS10.in -rw-r--r-- 1 will staff 1114 Apr 13 18:55 models-no-qheat/dustrad-n03-MS20.in -rw-r--r-- 1 will staff 1114 Apr 13 18:55 models-no-qheat/dustrad-n03-MS40.in -rw-r--r-- 1 will staff 1113 Apr 13 18:55 models-no-qheat/dustrad-n04-BSG.in -rw-r--r-- 1 will staff 1102 Apr 13 18:55 models-no-qheat/dustrad-n04-MS10.in -rw-r--r-- 1 will staff 1116 Apr 13 18:55 models-no-qheat/dustrad-n04-MS20.in -rw-r--r-- 1 will staff 1116 Apr 13 18:55 models-no-qheat/dustrad-n04-MS40.in
- Folders models-pah/ and models-pah-no-qheat/
import glob
infiles = glob.glob("models/dustrad*.in")
for infile in infiles:
with open(infile) as f:
s = f.read()
# Add in PAH grains with constant abundance
s = s.replace(
"grains ism function sublimation",
"\n".join(["grains ism function sublimation",
"grains PAH",
"set PAH constant"]))
with open(infile.replace("models/", "models-pah/"), "w") as f:
f.write(s)
# Change both grains commands to turn off quantum heating
for grain_cmd in "grains ism function sublimation", "grains PAH":
s = s.replace(grain_cmd, grain_cmd + " no qheat")
with open(infile.replace("models/", "models-pah-no-qheat/"), "w") as f:
f.write(s)cd models-no-qheat
export D=~/Work/CLOUDY/cloudy/source/sys_gcc
make -j6- The models with
no qheatare noticeably faster, but less than twice as fast - Average speedup factor = 7694.88/4502.76 = 1.71
grep 'ExecTime(s)' models-no-qheat/dustrad*.out |cut -f 1,11 -d' '| models-no-qheat/dustrad-n04-MS40.out: | 235.93 |
|---|---|
| 4502.76 |
grep 'ExecTime(s)' models/dustrad*.out |cut -f 1,11 -d' '| models/dustrad-n04-MS40.out: | 1051.74 |
|---|---|
| 7694.88 |
- Compare with the models with qheat off and on
- I did this originally in Grain emissivity versus U below
- These are based on the programs there
- Summary combination figure
- First I do 70 micron
- No-qheat models are shown with dashed lines
- Very slight differences for U < 1, but otherwise identical
- Not a surprise
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'no-qheat-j70-vs-U.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4,4))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth70 = 1.2491e12
Jy = 1e-23
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
mqh = CloudyModel(f'models-pah/{prefix}')
mnoqh = CloudyModel(f'models-pah-no-qheat/{prefix}')
except:
continue
for m, ls in [mqh, '-'], [mnoqh, '--']:
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
U[eden < 0.5*hden] = np.nan
j70 = m.data['emis']['PAC1 70.0000m'] / (4*np.pi*hden*bandwidth70*Jy)
ax.plot(U, j70, alpha=0.8, ls=ls, lw=1, color=col)
ax.plot(DL07tab['U']/8.0, DL07tab['70'], '.', alpha=0.5, color='r')
ax.set(
yscale='log',
xscale='log',
xlabel=r'Radiation field: $U = L_*\, /\, (4 \pi \,R^{2} \, c \, u_{\mathrm{MMP83}})$',
ylabel=r'Emissivity at 70 $\mu$m: $j_\nu$ [Jy cm$^{2}$ sr$^{-1}$ H$^{-1}$]',
#xlim=[1.0, 1.e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'no-qheat-j24-vs-U.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4,4))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth24 = 2.927e12
Jy = 1e-23
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
mqh = CloudyModel(f'models-pah/{prefix}')
mnoqh = CloudyModel(f'models-pah-no-qheat/{prefix}')
except:
continue
for m, ls in [mqh, '-'], [mnoqh, '--']:
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
U[eden < 0.5*hden] = np.nan
j24 = m.data['emis']['MIPS 24.0000m'] / (4*np.pi*hden*bandwidth24*Jy)
ax.plot(U, j24, alpha=0.8, ls=ls, lw=1, color=col)
ax.plot(DL07tab['U']/8.0, DL07tab['24'], '.', alpha=0.5, color='r')
ax.set(
yscale='log',
xscale='log',
xlabel=r'Radiation field: $U = L_*\, /\, (4 \pi \,R^{2} \, c \, u_{\mathrm{MMP83}})$',
ylabel=r'Emissivity at 24 $\mu$m: $j_\nu$ [Jy cm$^{2}$ sr$^{-1}$ H$^{-1}$]',
#xlim=[1.0, 1.e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'no-qheat-j8-vs-U.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4,4))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth8 = 1.3886e13
Jy = 1e-23
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
mqh = CloudyModel(f'models-pah/{prefix}')
mnoqh = CloudyModel(f'models-pah-no-qheat/{prefix}')
except:
continue
for m, ls in [mqh, '-'], [mnoqh, '--']:
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
U[eden < 0.5*hden] = np.nan
j8 = m.data['emis']['IRAC 8.00000m'] / (4*np.pi*hden*bandwidth8*Jy)
ax.plot(U, j8, alpha=0.8, ls=ls, lw=1, color=col)
#ax.plot(DL07tab['U'], DL07tab['8'], '.', alpha=0.5, color='k')
ax.plot(DL07tab['U']/8.0, DL07tab['8'], '.', alpha=0.5, color='r')
ax.set(
yscale='log',
xscale='log',
xlabel=r'Radiation field: $U = L_*\, /\, (4 \pi \,R^{2} \, c \, u_{\mathrm{MMP83}})$',
ylabel=r'Emissivity at 8 $\mu$m: $j_\nu$ [Jy cm$^{2}$ sr$^{-1}$ H$^{-1}$]',
#xlim=[1.0, 1.e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)- At fixed wavelength, emissivity is proportional to B_ν(T)
- B_ν has peak at 5.879e10 (T/K) Hz => λ = (5100 / T) µm
- Note that this is 1.7 times longer wavelength than the peak in B_λ
- So Planck peak is at 70 µm when T = 5100/70 = 72.9 K
- From Fig 3 of Paper III, this is when U ~ 1e3
- Default Cloudy PAH grains is in 10 size bins
- Described in Abel:2008a sec 3.1.3
- Abundance goes as a^-3.5
- From 30 to 500 carbon atoms
- That is about 0.4 to 1.0 nm
- Total carbon atom abundance in PAHs is n_C/n_H = 3e-6
- Note that Cloudy ISM grain mixture has smallest bin that starts at 5 nm
- Therefore there are no grains from 1 to 5 nm
- This is going to be important at 24 micron because From Fig 14 of Draine:2001a, it is the biggest PAH molecules that emit most there (e.g., the C_4000 H_1000 grain in the figure)
- Also, Fig 6 of DL07 shows that for λ > 12 µm it is grains with a > 10 Å = 1 nm that dominate the emission
- PAH grains in DL07
- We are using the
MW3.1_00models, which have q_PAH = 0.47%- See DL07-data/README.txt
- q_PAH is “the fraction of the total dust mass that is contributed by PAH particles containing < 1000 C atoms”
- We are using the
- Try and find equivalent q_PAH for Cloudy models
- graphite grains
- graphite_ism_10.opc
- Material density 2.26 g/cm^3
- Total grain volume 2.456864e-27 cm^3/H
- Atomic weight of Carbon: 12.01078 amu
- => Carbon abundance in grains = 2.26 2.456864e-27 / 12.01078 amu = 2.784e-4
- Cross-check: this should be the same as “depletion efficiency” times “abundance of grain molecule at max depletion”
- 4.800000e-04 5.800000e-01 = 2.784e-4 Exactly the same
- I was going to compare that with the PAH carbon abundance n_C/n_H = 3e-6 from above, but I realised I need to add in silicates too
- It would actually be asier to just stick with masses
- PAH grains
- pah1_ab08_10.opc
- Density: 1.80 g/cm^3
- Total volume: 3.300793e-29 cm^3/H
- Silicate grains
- silicate_ism_10.opc
- Density: 3.30 g/cm^3
- Total volume: 2.842704e-27 cm^3/H
- PAH mass fraction
- q_PAH = m(PAH) / (m(PAH) + m(graphite) + m(silicate))
- 1.80 3.300793e-29 / (1.80 3.300793e-29 + 2.26 2.456864e-27 + 3.30 2.842704e-27)
- => q_PAH = 3.96284052387e-3 = 0.40%
- This is similar to DL07 value, but slightly smaller
- graphite grains
- [2019-04-15 Mon] Suggested by Chip Kobulnicky
- Adapt the stellar models above but using the ISRF
- One option would be to use
table ism extinguish - But I decided it would be better to use the exact SED from MMP83
- This is written to the file models-isrf/isrf_mmp1983.sed
- See Compare SED between OB stars and interstellar field
- One option would be to use
- Luminosity is set to 1e6 Lsun in range 0.0114 → 1.0 Rydberg, which is 912 Å to 8 µm
- This means that radiation field is 1470 / R^2 with R in pc
- So to cover U = 0.1 → 1e6, we use R = 0.0383 → 120 pc
- Also: don’t iterate, since we don’t really need the diffuse fields. In fact, it would make the comparison easier if there were no diffuse radiation
- The emissivity from this model is compared with DL07 below in Cloudy emissivity exposed to MMP83 ISRF SED
from textwrap import dedent
import numpy as np
LSUN = 3.82e33
def isrf_input(L4=100.0):
s = f"# MMP83 interstellar radiation field, L = {1e4*L4:.2e} Lsun" + "\n"
s += f"table sed \"isrf_mmp1983.sed\"" + "\n"
s += f"luminosity {np.log10(LSUN*1e4*L4):.2f} range 0.0114 1.0" + "\n"
return s
def hden_input(hden):
return f"hden {hden:.2f} # density of {10**hden} pcc" + "\n"
def file_stem(hden, id_):
return f"dustrad-n{int(hden):02d}-{id_}"
def intro_input(hden, id_):
s = f"title Dust radiative acceleration: star {id_}, density {10**hden} pcc" + "\n"
s += f"set save prefix \"{file_stem(hden, id_)}\"" + "\n"
return s
def radius_input(L4, hden):
# Scale R_in to give U = 1e9
R_in_pc = 0.0383
R_out_pc = 120.0
s = f"# Start close in: {R_in_pc:.3f} pc" + "\n"
s += f"radius {np.log10(3.085677582e18*R_in_pc)}" + "\n"
s += f"stop radius {np.log10(3.085677582e18*R_out_pc)}" + "\n"
return s
outro_input = dedent("""\
# We stop on radius, so make sure not on T
stop temperature 4 K linear
sphere
grains ism
abundances HII region no grains
cosmic ray background
cmb
# save all the output we want
save overview last ".ovr"
save physical conditions last ".phys"
save continuum last ".cont"
# save continuum last every ".zcont"
save radius last ".rad"
save grain abundance last ".gabun"
save grain charge last ".gcharge"
save grain continuum last ".gcont"
save grain drift velocity last ".gdrift"
save grain potential last ".gpot"
save grain temperature last ".gtemp"
save lines, emissivity last ".emis"
H 1 6562.81A
Ca B 6562.81A
N 2 6583.45A
O 3 5006.84A
IRAC 3.60000m
IRAC 4.50000m
IRAC 5.80000m
IRAC 8.00000m
F12 12.0000m
MIPS 24.0000m
PAC1 70.0000m
PAC3 160.000m
end of lines
""")
id_ = "isrf"
for hden in [0.0]:
for L4 in [100.0]:
cloudy_input = (intro_input(hden, id_)
+ isrf_input(L4)
+ hden_input(hden)
+ radius_input(L4, hden)
+ outro_input)
file_name = f"models-isrf/{file_stem(hden, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)- [2018-11-04 Sun] Try and reproduce the SEDs of the potential bow wave candidates
- See
- Originally I had:
- B1.5V
- T = 23000 K
- L = 5600 Lsun
- But Alecian:2013a do a detailed spectral fit and find
- T = 20,000 +/- 1000 K
- V = 8.46
- (B - V) = 0.09
- Intrinsic (B - V)_0 = 0.25?
- A_V = 1.55
- They got that from assuming R_V = 5
- But M42 foreground dust is generally held to have R_V = 5.5
- Which would give A_V = 1.705 +/- 0.155
- => V_0 = 6.91
- Or V_0 = 6.755 +/- 0.155 with R_V = 5.5 +/- 0.5
- They got that from assuming R_V = 5
- Assume D = 375 +/- 30 pc
- But 410 +/- 10 pc would be better
- That would make 20% difference to luminosity
- Distance modulus 5 log D - 5 = 8.064 +/- 0.053
- From Malagnini:1986a Fig 5 and Table 4
- @ 20,000 +/- 1000 K, BC = -1.96 +/- 0.1
- => m_bol = 6.755 +/- 0.155 - 1.96 +/- 0.1 = 4.795 +/- 0.18
- => M_bol = 4.795 +/- 0.18 - 8.064 +/- 0.053 = -3.269 +/- 0.188
- Sun’s bolometric magnitude is 4.74
- L/Lsun = 10**(0.4 (4.74 + 3.269 +/- 0.188)) = 1600 +/- 300
- This is pretty much identical to what Alecian:2013a have in their Table 2
- log L/Lsun = 3.22 +/- 0.07 => L/Lsun = 1660 +/- 270
- We think the shell is neutral, so we need to go deeper than previously
- We can try the
constant pressure set LOGVALUEoption to CloudyLOGVALUEis log_10 P/k at illuminated face in units of cm^-3 K- Density still needs to be set, but it is only used in initial thermal solution, then discarded
- This includes radiative acceleration in the pressure balance, so it should automatically produce a radiation-pressurized shell
- The question is, whether to try and include the internal ionized zone that is seen in Hα
- To start with, we won’t
- We know that column is optical depth of 0.2 to 0.3 in FUV and optical
- N = 0.3 / m κ = 1.38e21 / κ_100
- Maybe easiest option is to set density as 20 times lower than our estimate for the neutral density
- So this would be 1e4
- We use Orion grains
from textwrap import dedent
import numpy as np
LSUN = 3.82e33
def star_input(id_, L4, T3, log_g):
s = f"# {id_} star" + "\n"
s += f"# T_eff = {1e3*T3:.0f} K, L = {1e4*L4:.2e} L_sun, log(g) = {log_g:.2f}" + "\n"
s += f"table star tlusty OBstar 3-dim {1e3*T3:.0f} {log_g} 0.0" + "\n"
s += f"luminosity total {np.log10(LSUN*1e4*L4):.2f}" + "\n"
return s
def hden_input(hden):
s = f"hden {hden:.2f} # density of {10**hden} pcc" + "\n"
s += "constant pressure" + "\n"
s += "stop AV 0.3" + "\n"
s += "stop temperature 100 K linear" + "\n"
return s
def file_stem(hden, R_pc, id_):
return f"shell-R{int(1000*R_pc):03d}-n{int(10*hden):02d}-{id_}"
def intro_input(hden, R_pc, id_):
s = f"title Dusty shell: star {id_}, R_in = {R_pc} pc, density {10**hden} pcc" + "\n"
s += f"set save prefix \"{file_stem(hden, R_pc, id_)}\"" + "\n"
return s
stars = [
# id_ L4 T3 log_g
["LP_Ori", 0.16, 20.0, 4.0],
]
def radius_input(R_pc):
R_in = 3.085677582e18*R_pc
s = f"# Start at LP Ori inner radius: {R_pc:.5f} pc" + "\n"
s += f"radius {np.log10(R_in)}" + "\n"
return s
outro_input = dedent("""\
iterate
grains orion function sublimation
abundances HII region no grains
cosmic ray background
cmb
table ism
# save all the output we want
save overview last ".ovr"
save physical conditions last ".phys"
save continuum last ".cont"
# save continuum last every ".zcont"
save radius last ".rad"
save grain abundance last ".gabun"
save grain charge last ".gcharge"
save grain continuum last ".gcont"
save grain drift velocity last ".gdrift"
save grain potential last ".gpot"
save grain temperature last ".gtemp"
""")
for R_pc, hden in [
[0.01, 3.5], [0.01, 4.0], [0.01, 4.5],
[0.005, 3.0], [0.003, 2.5]
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)- So we have some models that start at 0.01 pc, and some that start at smaller radii
- But we want them all to get to the right density in the shell
- New models, with some tweaks
- [X] Add a hotter star
- [X] Add PAHs
grains PAH- By default, this is just in atomic H zone (we have no H_2 so that doesn’t matter)
- [X] Try and save some emissivity profiles
- All line labels are saved in models/cloudy-line-labels.dat
- This includes continuum bands
- [X] Refine the densities, so they all have shell at the same radius
- We have a few models that all end up with the shell at 0.01 pc
R003-n29-LP_Ori20,R005-n30-LP_Ori20,R001-n25-LP_Ori22- The last one is with a hotter star and has a higher ionized density, but identical neutral shell density
- We have a few models that all end up with the shell at 0.01 pc
from textwrap import dedent
import numpy as np
LSUN = 3.82e33
def star_input(id_, L4, T3, log_g):
s = f"# {id_} star" + "\n"
s += f"# T_eff = {1e3*T3:.0f} K, L = {1e4*L4:.2e} L_sun, log(g) = {log_g:.2f}" + "\n"
s += f"table star tlusty OBstar 3-dim {1e3*T3:.0f} {log_g} 0.0" + "\n"
s += f"luminosity total {np.log10(LSUN*1e4*L4):.2f}" + "\n"
return s
def hden_input(hden, AV=0.3):
s = f"hden {hden:.2f} # density of {10**hden} pcc" + "\n"
s += "constant pressure" + "\n"
s += f"stop AV {AV:.2f}" + "\n"
s += "stop temperature 10 K linear" + "\n"
s += "iterate" + "\n"
return s
def file_stem(hden, R_pc, id_, extra=""):
"""Construct file name from radius, density and star id"""
return f"shell-R{int(1000*R_pc):03d}-n{int(10*hden):02d}-{id_}" + extra
def intro_input(hden, R_pc, id_, extra=""):
s = f"title Dusty shell: star {id_}, R_in = {R_pc} pc, density {10**hden} pcc" + "\n"
s += f"set save prefix \"{file_stem(hden, R_pc, id_, extra)}\"" + "\n"
return s
def radius_input(R_pc):
R_in = 3.085677582e18*R_pc
s = f"# Start at LP Ori inner radius: {R_pc:.5f} pc" + "\n"
s += f"radius {np.log10(R_in)}" + "\n"
return s
def magnetic_input(hden, hden1=3.0, vA=2.0, gamma_m=4.0/3.0):
"""
Set B field at illuminated face (where n = 10**`hden`) so that
Alfven speed is `vA` km/s for typical densities 10**`hden1`,
assuming Pmag ~ rho**`gamma_m`
"""
m = 1.3*1.67262158e-24
# This should give 33 micro G for vA=2 and hden1=3
B1 = np.sqrt(4*np.pi*m*10**hden1)*vA*1.0e5
B0 = B1*10**(0.5*gamma_m*(hden - hden1))
s = dedent(f"""\
# Magnetic field of {1e6*B0:.3f} microG to give Alfven speed
# of {vA:.3f} km/s for density of {10**hden1} pcc
""")
s += f"magnetic field, log(B) = {np.log10(B0):.3f}, tangled {gamma_m:.5f}" + "\n"
# s += "# turbulence equipartition" + "\n"
return s
extra_input = dedent("""\
cosmic ray background
cmb
table ism
""")
save_input = dedent("""\
# save all the output we want
save overview last ".ovr"
save pressure last ".pre"
save physical conditions last ".phys"
save continuum last ".cont"
# save continuum last every ".zcont"
save radius last ".rad"
save grain abundance last ".gabun"
save grain charge last ".gcharge"
save grain continuum last ".gcont"
save grain drift velocity last ".gdrift"
save grain potential last ".gpot"
save grain temperature last ".gtemp"
save lines, emissivity last ".emis"
H 1 6562.81A
Ca B 6562.81A
N 2 6583.45A
O 3 5006.84A
IRAC 3.60000m
IRAC 4.50000m
IRAC 5.80000m
IRAC 8.00000m
F12 12.0000m
F25 25.0000m
PAC1 70.0000m
PAC3 160.000m
end of lines
""")
dust_input = dedent("""\
grains orion function sublimation
grains PAH
abundances HII region no grains
""")
outro_input = extra_input + dust_input + save_input<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20", 0.16, 20.0, 4.0],
["LP_Ori22", 0.16, 22.0, 4.0],
]
for R_pc, hden in [
[0.005, 3.0], [0.003, 2.9], [0.001, 2.5],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)Finally, try different shell thickness
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20thick", 0.16, 20.0, 4.0],
["LP_Ori22thick", 0.16, 22.0, 4.0],
]
for R_pc, hden in [
[0.005, 3.0], [0.003, 2.9], [0.001, 2.5],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, 0.5)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)- And lower dust opacity too
- We need to put down the stopping temperature, since it gets to 100 K very soon in the PDR
- Also ditch the hotter star
- And compensate for less of a radiation-pressure hole in the low-radius model
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20lowZ", 0.16, 20.0, 4.0],
]
# Multiply all grain component abundances by 0.1
grain_commands = [_ for _ in outro_input.split("\n") if _.startswith("grains")]
for cmd in grain_commands:
outro_input = outro_input.replace(cmd, cmd + " 0.1")
for R_pc, hden in [
[0.005, 3.0], [0.003, 2.9], [0.001, 2.8],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)- And add a magnetic field with γ = 4/3
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20B", 0.16, 20.0, 4.0],
]
for R_pc, hden in [
[0.005, 3.0], [0.003, 2.9], [0.001, 2.8],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden)
+ magnetic_input(hden)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)- That makes a shell that is too thick, so try again with magnetic γ = 1, which corresponds to constant Alfvén speed
- For some reason γ = 1 is not allowed by Cloudy, so use γ = 1.001
- Also remove the turbulence and fix the density of the low radius model
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20BB", 0.16, 20.0, 4.0],
]
for R_pc, hden in [
[0.005, 3.0], [0.003, 2.9],
[0.001, 2.3],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden)
+ magnetic_input(hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)And one final model where we go to A_V = 1.7 just to get the extinguished incident spectrum
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20AV", 0.16, 20.0, 4.0],
]
for R_pc, hden in [
[0.001, 2.3],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, 1.7)
+ magnetic_input(hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)Still another: combo model with increased thickness (A_V = 0.5), reduced dust cross section, and magnetic field
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20BZ5", 0.16, 20.0, 4.0],
]
# Multiply all grain component abundances by 0.1
grain_commands = [_ for _ in outro_input.split("\n") if _.startswith("grains")]
for cmd in grain_commands:
outro_input = outro_input.replace(cmd, cmd + " 0.1")
for R_pc, hden in [
[0.001, 2.8],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, 0.5)
+ magnetic_input(hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)And one with increased thickness and B field, but normal dust
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20BB5", 0.16, 20.0, 4.0],
]
for R_pc, hden in [
[0.001, 2.5],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, 0.5)
+ magnetic_input(hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)And one with intermediate dust reduction: three times smaller instead of ten
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["LP_Ori20Bz5", 0.16, 20.0, 4.0],
]
# Multiply all grain component abundances by 0.333
grain_commands = [_ for _ in outro_input.split("\n") if _.startswith("grains")]
for cmd in grain_commands:
outro_input = outro_input.replace(cmd, cmd + " 0.333")
for R_pc, hden in [
[0.001, 2.7],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, 0.5)
+ magnetic_input(hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)- Stellar data from Simon-Diaz:2006b
- Use the same infrastructure as for LP Ori
First attempt
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D", 2.95, 32.0, 4.2],
]
for R_pc, hden, AV in [
[0.003, 3.3, 0.02], [0.003, 3.5, 0.02], [0.003, 3.7, 0.02],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)Put up the density to 1e5 and decrease the dust opacity by factor of 3, and put up the shell thickness a bit
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-z", 2.95, 32.0, 4.2],
]
# Multiply all grain component abundances by 0.333
grain_commands = [_ for _ in outro_input.split("\n") if _.startswith("grains")]
for cmd in grain_commands:
outro_input = outro_input.replace(cmd, cmd + " 0.333")
for R_pc, hden, AV in [
[0.003, 4.5, 0.03], [0.003, 5.0, 0.03], [0.003, 5.0, 0.05],
]:
for id_, L4, T3, log_g in stars:
extra = f"-AV{int(100*AV):1d}"
cloudy_input = (intro_input(hden, R_pc, id_, extra)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ outro_input)
file_name = f"models/{file_stem(hden, R_pc, id_, extra)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)Try single size dust - first 1 micron
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-1m000", 2.95, 32.0, 4.2],
]
dust_input = dedent("""\
grains "silicate_1m000.opc" function sublimation
abundances HII region no grains
""")
for R_pc, hden, AV in [
[0.003, 4.5, 0.03], [0.003, 5.0, 0.03], [0.003, 5.0, 0.05],
]:
for id_, L4, T3, log_g in stars:
extra = f"-AV{int(100*AV):1d}"
cloudy_input = (intro_input(hden, R_pc, id_, extra)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ extra_input
+ dust_input
+ save_input)
file_name = f"models/{file_stem(hden, R_pc, id_, extra)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)Now 0.1 micron
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-0m100", 2.95, 32.0, 4.2],
]
dust_input = dedent("""\
grains "silicate_0m100.opc" function sublimation 0.3
abundances HII region no grains
""")
for R_pc, hden, AV in [
[0.003, 4.5, 0.03], [0.003, 5.0, 0.03], [0.003, 5.0, 0.05],
]:
for id_, L4, T3, log_g in stars:
extra = f"-AV{int(100*AV):1d}"
cloudy_input = (intro_input(hden, R_pc, id_, extra)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ extra_input
+ dust_input
+ save_input)
file_name = f"models/{file_stem(hden, R_pc, id_, extra)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)Now both 0.1 micron (0.1 abundance) and 1 micron
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-twin", 2.95, 32.0, 4.2],
]
dust_input = dedent("""\
grains "silicate_1m000.opc" function sublimation 1.0
grains "silicate_0m100.opc" function sublimation 0.1
abundances HII region no grains
""")
for R_pc, hden, AV in [
[0.003, 4.5, 0.03], [0.003, 5.0, 0.03], [0.003, 5.0, 0.05],
]:
for id_, L4, T3, log_g in stars:
extra = f"-AV{int(100*AV):1d}"
cloudy_input = (intro_input(hden, R_pc, id_, extra)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ extra_input
+ dust_input
+ save_input)
file_name = f"models/{file_stem(hden, R_pc, id_, extra)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)And add in a small amount of 0.01 micron too for the IRAC bands, at the same reducing the 0.1 micron by 2
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-triple", 2.95, 32.0, 4.2],
]
dust_input = dedent("""\
grains "silicate_1m000.opc" function sublimation 1.0
grains "silicate_0m100.opc" function sublimation 0.03
grains "silicate_0m010.opc" function sublimation 0.002
abundances HII region no grains
""")
for R_pc, hden, AV in [
[0.001, 4.5, 0.03], [0.003, 5.0, 0.03], [0.003, 5.0, 0.05],
]:
for id_, L4, T3, log_g in stars:
extra = f"-AV{int(100*AV):1d}"
cloudy_input = (intro_input(hden, R_pc, id_, extra)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ extra_input
+ dust_input
+ save_input)
file_name = f"models/{file_stem(hden, R_pc, id_, extra)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)And one final model with the foreground extinction: AV = 2.0
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-AV200", 2.95, 32.0, 4.2],
]
dust_input = dedent("""\
grains orion function sublimation 0.333
grains PAH 0.333
abundances HII region no grains
""")
for R_pc, hden, AV in [
[0.003, 5.0, 2.0], [0.003, 4.0, 2.0],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ extra_input
+ dust_input
+ save_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)That was not the end! Now, I will do some lower luminosity models, since doing Audit of th1D fluxes in ../dust-wave-case-studies
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-L25-triple", 2.5, 32.0, 4.2],
]
dust_input = dedent("""\
grains "silicate_1m000.opc" function sublimation 1.0
grains "silicate_0m100.opc" function sublimation 0.015
grains "silicate_0m010.opc" function sublimation 0.002
grains "graphite_0m010.opc" function sublimation 0.0003
abundances HII region no grains
""")
for R_pc, hden, AV in [
[0.003, 4.7, 0.05], [0.003, 5.0, 0.05], [0.003, 5.0, 0.08],
]:
for id_, L4, T3, log_g in stars:
extra = f"-AV{int(100*AV):1d}"
cloudy_input = (intro_input(hden, R_pc, id_, extra)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ extra_input
+ dust_input
+ save_input)
file_name = f"models/{file_stem(hden, R_pc, id_, extra)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)And the AV = 2 model for the foreground
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-L25-AV200", 2.5, 32.0, 4.2],
]
dust_input = dedent("""\
grains orion function sublimation
grains PAH
abundances HII region no grains
""")
for R_pc, hden, AV in [
[0.003, 5.0, 2.0],
]:
for id_, L4, T3, log_g in stars:
cloudy_input = (intro_input(hden, R_pc, id_)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ extra_input
+ dust_input
+ save_input)
file_name = f"models/{file_stem(hden, R_pc, id_)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)And a few more, slightly higher density
<<lp-ori-cloudy-functions>>
stars = [
# id_ L4 T3 log_g
["th1D-L25-triple", 2.5, 32.0, 4.2],
]
dust_input = dedent("""\
grains "silicate_1m000.opc" function sublimation 1.0
grains "silicate_0m100.opc" function sublimation 0.015
grains "silicate_0m010.opc" function sublimation 0.002
grains "graphite_0m010.opc" function sublimation 0.0003
abundances HII region no grains
""")
for R_pc, hden, AV in [
[0.003, 5.1, 0.05], [0.003, 5.2, 0.05], [0.003, 5.3, 0.05],
]:
for id_, L4, T3, log_g in stars:
extra = f"-AV{int(100*AV):1d}"
cloudy_input = (intro_input(hden, R_pc, id_, extra)
+ star_input(id_, L4, T3, log_g)
+ hden_input(hden, AV)
+ magnetic_input(hden, hden1=hden, vA=2.0, gamma_m=1.001)
+ radius_input(R_pc)
+ extra_input
+ dust_input
+ save_input)
file_name = f"models/{file_stem(hden, R_pc, id_, extra)}.in"
with open(file_name, "w") as f:
f.write(cloudy_input)- This makes
.xsecfiles in the dust-opacity folder, which have just the cross-section data from the opacity files, for ease of reading
import glob
import sys
import os
opc_files = glob.glob("/Users/will/Work/CLOUDY/cloudy/data/*.opc")
ABS_HEADER = "# anu (Ryd) abs_cs_01 (cm^2/H) abs_cs_02....."
SCA_HEADER = "# anu (Ryd) sct_cs_01 (cm^2/H) sct_cs_02....."
GGG_HEADER = "# anu (Ryd) (1-g)_bin_01 (1-g)_bin_02....."
LAST_LINE = "# anu (Ryd) inverse attenuation length (cm^-1)"
for ofile in opc_files:
with open(ofile) as f:
text = f.read()
i1 = text.find(ABS_HEADER)
i2 = text.find(SCA_HEADER)
i3 = text.find(GGG_HEADER)
i4 = text.find(LAST_LINE)
nfile = os.path.join("dust-opacity",
os.path.basename(ofile))
with open(nfile.replace(".opc", ".abs"), "w") as f:
f.write(text[i1:i2])
with open(nfile.replace(".opc", ".sca"), "w") as f:
f.write(text[i2:i3])
with open(nfile.replace(".opc", ".ggg"), "w") as f:
f.write(text[i3:i4])from astropy.table import Table
from astropy.io.ascii import InconsistentTableError
import glob
# File extensions that might be present, but which are NOT Cloudy save files
IGNORE_EXTS = ["pdf", "png", "jpg"]
class CloudyModel(object):
"""Lightweight wrapper for output from Cloudy run
For example:
>>> from cloudytab import CloudyModel
>>> m = CloudyModel("myfolder/mymodel")
`m.files` contains a list of all the files that were found:
`['myfolder/mymodel.in', 'myfolder/mymodel.ovr', ETC]`
`m.data` contains dict of astropy.Table's, one for each save file:
`{'ovr': <Table length=289> ..., ETC}`
`m.io['in']` and `m.io['out']` contain the input and output streams
"""
def __init__(self, prefix):
self.files = glob.glob(prefix + ".*")
self.data = {}
self.io = {}
for file_ in self.files:
saveid = file_.split(".")[-1]
if saveid in IGNORE_EXTS:
# Figure files, etc need to be skipped
pass
elif saveid in ["in", "out"]:
# Special case of input and output files
with open(file_) as f:
# Just save the whole file as a string
self.io[saveid] = f.read()
else:
# Assume all else are save files
try:
self.data[saveid] = Table.read(
file_, delimiter="\t", guess=False, fast_reader=False,
format="ascii.commented_header")
except UnicodeDecodeError:
# Binary files can raise this error - ignore them
pass
except InconsistentTableError:
# The "save heating" files can raise this error - skip them
pass
- Estimate ionization parameter from H neutral fraction
- Do a single star at a time, since there may be a secondary dependence on the spectral shape
import glob
from matplotlib import pyplot as plt
import seaborn as sns
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
figfile = f"phi-ipar-{STAR}-{GRAIN}.pdf"
infiles = glob.glob(f"models/dustrad-*-{STAR}.in")
models = []
for infile in infiles:
prefix = infile.replace(".in", "")
modelid = prefix.replace("models/dustrad-", "") # e.g., n03-MS10
models.append([modelid, CloudyModel(prefix)])
fig, ax = plt.subplots()
for label, m in sorted(models):
# Ionization parameter, estimated as x^2 / (1 - x)
# (Initially, we neglect correction for alpha(T) and sigma(tau))
ipar = m.data["ovr"]["HII"]**2 / m.data["ovr"]["HI"]
# Grain potential divided by kT
gpot = m.data["gpot"][GRAIN]*u.eV / (m.data["ovr"]["Te"]*u.K*k_B).to(u.eV)
ax.plot(ipar, gpot, label=label)
ax.axvspan(0.0111, 8.1, color='k', alpha=0.1) # x = 0.1 -> 0.9
ax.axhspan(-1.0, 1.0, color='k', alpha=0.1) # |phi| < 1
ax.legend(title=GRAIN)
ax.set(
xscale='log',
yscale='symlog',
xlabel="Ionization parameter",
ylabel="Grain potential / k T",
xlim=[3e-5, 3e6],
ylim=[-5.0, 50.0],
)
sns.despine()
fig.savefig(figfile)import glob
from matplotlib import pyplot as plt
import seaborn as sns
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
figfile = f"phi-ipar-{STAR}-allgrain.pdf"
infiles = glob.glob(f"models/dustrad-*-{STAR}.in")
models = []
for infile in infiles:
prefix = infile.replace(".in", "")
modelid = prefix.replace("models/dustrad-", "") # e.g., n03-MS10
models.append([modelid, CloudyModel(prefix)])
sns.set_color_codes("deep")
fig, ax = plt.subplots()
colors = sns.color_palette(palette="magma_r", n_colors=len(models))
fastlabel = r"$w_\mathrm{drift} > 10$ km/s"
fastishlabel = r"$w_\mathrm{drift} > 1$ km/s"
for color, [label, m] in zip(colors, sorted(models)):
# Ionization parameter, estimated as x^2 / (1 - x)
# (Initially, we neglect correction for alpha(T) and sigma(tau))
ipar = m.data["ovr"]["HII"]**2 / m.data["ovr"]["HI"]
for grain in m.data["gpot"].colnames[1:]:
# Grain potential divided by kT
gpot = m.data["gpot"][grain]*u.eV / (m.data["ovr"]["Te"]*u.K*k_B).to(u.eV)
ls = '--' if grain.startswith("gra") else '-'
# Only plot where grains are not sublimated
mm = m.data["gabun"][grain] >= 0.3*m.data["gabun"][grain].max()
ax.plot(ipar[mm], gpot[mm], alpha=0.8, color=color, ls=ls, lw=0.4, label=label)
m2 = m.data["gdrift"][grain] >= 10.0
m1 = (m.data["gdrift"][grain] >= 1.0) & ~m2
ax.scatter(ipar[m2 & mm], gpot[m2 & mm], label=fastlabel,
marker='.', s=60, alpha=0.8, color='c', edgecolors='none')
ax.scatter(ipar[m1 & mm], gpot[m1 & mm], label=fastishlabel,
marker='.', s=40, alpha=0.8, color='b', edgecolors='none')
label = '_nolabel_' # Only label first grain component
fastlabel = '_nolabel_'
fastishlabel = '_nolabel_'
ax.legend()
ax.axvspan(0.0111, 8.1, color='k', alpha=0.1) # x = 0.1 -> 0.9
ax.axhspan(-1.0, 1.0, color='k', alpha=0.1) # |phi| < 1
#ax.legend(title=GRAIN)
ax.text(0.0015, -3.0, "PDR", ha="center")
ax.text(0.3, -3.0, "Ionization\nfront", ha="center")
ax.text(200, -3.0, "H II region", ha="center")
ax.set(
xscale='log',
yscale='symlog',
xlabel="Hydrogen ionization: $x^{2} / (1 - x)$",
ylabel="Grain potential / $k T$",
xlim=[3e-4, 3e8],
ylim=[-5.0, 50.0],
)
sns.despine()
fig.savefig(figfile)- First of all, look at τ_ν in the final zone by using the continuum
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'tau-{DENID}.pdf'
fig, ax = plt.subplots()
for star in "MS10", "MS20", "MS40", "BSG":
prefix = f"dustrad-{DENID}-{star}"
m = CloudyModel(f'models/{prefix}')
nu = m.data['cont']['Cont nu']
nuFnu_inc = m.data['cont']['incident']
nuFnu_trans = m.data['cont']['trans']
nuFnu_tot = m.data['cont']['total']
tau_nu = -np.log(m.data['cont']['trans'] / m.data['cont']['incident'] )
ax.plot(nu, tau_nu, label=star)
ax.legend(title=DENID)
ax.axvspan(912.0/2000.0, 1.0, color='0.9')
ax.axvspan(1.0, 4.0, color='0.95')
ax.set(
xlim=[0.05, 4.0],
ylim=[0.004, 200.0],
yscale='log',
xscale='log',
xlabel='Photon energy, Rydberg',
ylabel=r'$\tau_{\nu}$',
)
fig.savefig(pltfile)Then we use these to find the local fluxes and determine radiation pressure
import numpy as np
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
stars = [
["MS10", 0.63],
["MS20", 5.45],
["MS40", 22.2],
["BSG", 30.2]
]
denids = [f"n0{_}" for _ in range(5)]
for star, L4 in stars:
L = 1e4*3.82e33*L4
for denid in denids:
# Read model
prefix = f"dustrad-{denid}-{star}"
m = CloudyModel(f'models/{prefix}')
# Find tau
nu = m.data['cont']['Cont nu']
nuFnu_inc = m.data['cont']['incident']
nuFnu_trans = m.data['cont']['trans']
tau_nu = -np.log(m.data['cont']['trans'] / m.data['cont']['incident'] )
# F_nu spectrum that is normalized to unit integral
Fnu_0 = nuFnu_inc / nu
Fnu_0 /= np.trapz(Fnu_0, nu)
# Masks for non-ionizing and ionizing radiation
mfuv = nu < 1.0
meuv = ~mfuv
# Scale of tau with radius, normalized on [0, 1]
# For FUV, it is just column density - proprtional to depth at constant density
depth = m.data['ovr']['depth']
tau_rscale_fuv = depth / depth[-1]
# For EUV, it is neutral column density
nzones = len(depth)
nH0 = m.data['ovr']['hden']*m.data['ovr']['HI']
tau_rscale_euv = np.array([np.trapz(nH0[:i], depth[:i]) for i in range(nzones)])
tau_rscale_euv /= tau_rscale_euv[-1]
# Extinction factor e^(-tau) as function of depth for fuv and euv
extinct_fuv = np.array(
[np.trapz(Fnu_0[mfuv]*np.exp(-tau_nu[mfuv]*tau_rscale_fuv[i]), nu[mfuv])
for i in range(nzones)])
extinct_euv = np.array(
[np.trapz(Fnu_0[meuv]*np.exp(-tau_nu[meuv]*tau_rscale_euv[i]), nu[meuv])
for i in range(nzones)])
radius = m.data['rad']['radius']
# Flux in each band
F_fuv = L * extinct_fuv / (4*np.pi*radius**2)
F_euv = L * extinct_euv / (4*np.pi*radius**2)
F_bol = F_fuv + F_euv
tab = Table(
[radius, F_bol, F_fuv, F_euv, extinct_fuv, extinct_euv, tau_rscale_fuv, tau_rscale_euv],
names=('R', 'F', 'F_F', 'F_E', 'E_F', 'E_E', 'T_F', 'T_E')
)
tab.write(f'models/{prefix}.flux',
format='ascii.commented_header',
formats={_: "%.4g" for _ in tab.colnames},
delimiter='\t', overwrite=True)
This works OK, but it has the disadvantage that it ignores the dust extinction in the EUV
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'fluxes-{PREFIX}.pdf'
fig, ax = plt.subplots()
F_habing = 1.6e-3
prefix = f"dustrad-{PREFIX}"
m = CloudyModel(f'models/{prefix}')
R = m.data['flux']['R'] / 3.085677582e18
hden = m.data['ovr']['hden']
ax.plot(R, m.data['flux']['F_F'] / (F_habing*hden), label='G(FUV) / n')
ax.plot(R, m.data['flux']['F_E'] / (F_habing*hden), label='G(EUV) / n')
ax.legend(title=PREFIX)
ax.set(
yscale='log',
xscale='log',
xlabel='Radius, pc',
ylabel=r'Flux',
ylim=[2e-5, 2e5],
xlim=[0.8*R[0], 3000*R[0]],
)
fig.savefig(pltfile)Check ion fraction versus F(EUV)/n
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'ion-params.pdf'
fig, ax = plt.subplots()
F_habing = 1.6e-3
light_speed = 2.99792458e10
kB = 1.3806503e-16
eV = 1.602176462e-12
stars = [
["MS10", 0.63, 1.3e-4, 'Purples_d' ],
["MS20", 5.45, 0.16, 'Oranges_d'],
["MS40", 22.2, 1.41, 'Blues_d' ],
["BSG", 30.2, 0.016, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
for star, L4, S49, cmap in stars:
L = 1e4*3.82e33*L4
L_EUV = 1e49*S49*13.6*eV
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
label = f"{star} {denid}"
m = CloudyModel(f'models/{prefix}')
R = m.data['rad']['radius']
G_n = m.data['flux']['F_E']/(F_habing*m.data['ovr']['hden'])
F_E0 = L_EUV / (4*np.pi*R**2)
Rm = R.max()
F_E0 *= (1. - (R/Rm)**3)
G_n0 = F_E0 / (F_habing*m.data['ovr']['hden'])
x = m.data['ovr']['HII']
ax.plot(G_n, x**2/(1 - x), alpha=0.8, lw=0.6, color=col, label=label)
#ax.plot(G_n0, x**2/(1 - x), alpha=0.6, lw=0.5)
ax.legend(ncol=2, fontsize="x-small")
ax.set(
yscale='log',
xscale='log',
xlabel='G(EUV) / n',
ylabel=r'$x^{2} / (1 - x)$',
)
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'drift-pratio.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(6,4))
F_habing = 1.6e-3
light_speed = 2.99792458e10
kB = 1.3806503e-16
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
m = CloudyModel(f'models/{prefix}')
Prad = m.data['flux']['F']/light_speed
R = m.data['rad']['radius']
Prad0 = L / (4*np.pi*light_speed*R**2)
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
Pgas[eden < 0.5*hden] = np.nan
for grain in m.data['gdrift'].colnames[1::5]:
ax.plot(Prad0/Pgas, m.data['gdrift'][grain], alpha=0.7, lw=0.5, color=col)
ax.set(
yscale='log',
xscale='log',
xlabel='Prad / Pgas',
ylabel=r'V drift',
xlim=[1.0e-2, 4e4],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'phi-pratio.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4, 3))
F_habing = 1.6e-3
light_speed = 2.99792458e10
kB = 1.3806503e-16
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
m = CloudyModel(f'models/{prefix}')
kT_eV = (m.data["ovr"]["Te"]*u.K*k_B).to(u.eV)
Prad = m.data['flux']['F']/light_speed
R = m.data['rad']['radius']
Prad0 = L / (4*np.pi*light_speed*R**2)
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
Pgas[eden < 0.5*hden] = np.nan
for igrain, grain in enumerate(m.data['gpot'].colnames[1::2]):
ls = '--' if 'gra' in grain else '-'
phi = m.data['gpot'][grain] / kT_eV
ax.plot(Prad0/Pgas, phi, alpha=0.3, ls=ls, lw=0.2+0.1*igrain, color=col)
p1, p2 = 1.0e-2, 4e4
pgrid = np.logspace(-2.0, 4.6)
phifit = 1.5*np.log(pgrid/0.1)
ax.plot(pgrid, phifit, lw=2, ls='--', color="k")
ax.plot(pgrid, 1.5*phifit, lw=1, color="k")
ax.plot(pgrid, phifit/1.5, lw=1, color="k")
ax.set(
yscale='linear',
xscale='log',
xlabel=r'Radiation parameter: $\Xi = P_\mathrm{rad} \,/\, P_\mathrm{gas}$',
ylabel=r'Grain potential: $\phi = U \,/\, k T$',
xlim=[pgrid[0], pgrid[-1]],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'drift-gn.pdf'
sns.set_color_codes("bright")
F_habing = 1.6e-3
light_speed = 2.99792458e10
kB = 1.3806503e-16
stars = [
["MS40", 22.2, 'Blues_d'],
["MS20", 5.45, 'Oranges_d'],
["MS10", 0.63, 'Purples_d'],
["BSG", 30.2, 'Greens_d']
]
select_grains = "sil-ism04", "sil-ism10", "gra-ism04", "gra-ism10"
glabel = {
"sil-ism04": "Silicate\n0.02 micron",
"sil-ism10": "Silicate\n0.2 micron",
"gra-ism04": "Graphite\n0.02 micron",
"gra-ism10": "Graphite\n0.2 micron"
}
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(10, 5))
denids = [f"n0{_}" for _ in range(5)]
for grain, ax in zip(select_grains, axes.flat):
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
m = CloudyModel(f'models/{prefix}')
Prad = m.data['flux']['F']/light_speed
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
Fbol = L / (4*np.pi*m.data['rad']['radius']**2)
G_n = m.data['flux']['F_F']/(F_habing*m.data['ovr']['hden'])
mm = m.data["gabun"][grain] >= 0.9*m.data["gabun"][grain].max()
msub = m.data['gdrift'][grain] < 20.0
if denid.endswith("00"):
label = f"{denid} {star}"
else:
label = f"{denid}"
ax.plot(G_n[mm & msub], m.data['gdrift'][grain][mm & msub],
alpha=0.85, lw=0.7, color=col, label=label)
ax.plot(G_n[mm & ~msub], m.data['gdrift'][grain][mm & ~msub],
alpha=0.85, lw=0.7, color=col, label="_nolabel_")
ax.axvline(1e4, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axvline(1e5, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axhline(1.0, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axhline(10.0, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.set_title(glabel[grain], fontsize="small", pad=-12)
axes[0, 0].legend(ncol=2, fontsize="xx-small", loc="left")
axes[1, 0].set(
yscale='log',
xscale='log',
xlabel=r'FUV radiation parameter: $G\, /\, n$, Habing cm$^3$',
ylabel=r'$V_\mathrm{drift}$, km/s',
xlim=[3.0e-1, 3.0e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'drift-pratio-4panel.pdf'
sns.set_color_codes("bright")
F_habing = 1.6e-3
light_speed = 2.99792458e10
kB = 1.3806503e-16
stars = [
["MS40", 22.2, 'Blues_d'],
["MS20", 5.45, 'Oranges_d'],
["MS10", 0.63, 'Purples_d'],
["BSG", 30.2, 'Greens_d']
]
select_grains = "sil-ism04", "sil-ism10", "gra-ism04", "gra-ism10"
glabel = {
"sil-ism04": "Silicate\n0.02 micron",
"sil-ism10": "Silicate\n0.2 micron",
"gra-ism04": "Graphite\n0.02 micron",
"gra-ism10": "Graphite\n0.2 micron"
}
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(10, 5))
denids = [f"n0{_}" for _ in range(5)]
for grain, ax in zip(select_grains, axes.flat):
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
m = CloudyModel(f'models/{prefix}')
Prad = m.data['flux']['F']/light_speed
R = m.data['rad']['radius']
Prad0 = L / (4*np.pi*light_speed*R**2)
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
Pgas[eden < 0.5*hden] = np.nan
Upsilon = Prad0/Pgas
Fbol = L / (4*np.pi*m.data['rad']['radius']**2)
G_n = m.data['flux']['F_F']/(F_habing*m.data['ovr']['hden'])
mm = m.data["gabun"][grain] >= 0.9*m.data["gabun"][grain].max()
msub = m.data['gdrift'][grain] < 20.0
if denid.endswith("00"):
label = f"{denid} {star}"
else:
label = f"{denid}"
ax.plot(Upsilon[mm & msub], m.data['gdrift'][grain][mm & msub],
alpha=0.85, lw=0.7, color=col, label=label)
ax.plot(Upsilon[mm & ~msub], m.data['gdrift'][grain][mm & ~msub],
alpha=0.85, lw=0.7, color=col, label="_nolabel_")
ax.axvline(300, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axvline(3000, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axhline(1.0, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axhline(10.0, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.set_title(glabel[grain], fontsize="small", pad=-12)
axes[0, 0].legend(ncol=2, fontsize="xx-small", loc="left")
axes[1, 0].set(
yscale='log',
xscale='log',
xlabel=r'Radiation parameter: $\Xi = P_{\mathrm{rad}} / P_{\mathrm{gas}}$',
ylabel=r'$w_\mathrm{drift}$, km/s',
xlim=[1.0e-2, 4e4],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)Same but for the smallest grains of all
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'drift-pratio-small-grains.pdf'
sns.set_color_codes("bright")
F_habing = 1.6e-3
light_speed = 2.99792458e10
kB = 1.3806503e-16
stars = [
["MS40", 22.2, 'Blues_d'],
["MS20", 5.45, 'Oranges_d'],
["MS10", 0.63, 'Purples_d'],
["BSG", 30.2, 'Greens_d']
]
select_grains = "sil-ism01", "sil-ism02", "gra-ism01", "gra-ism02"
glabel = {
"sil-ism01": "Silicate\n0.006 micron",
"sil-ism02": "Silicate\n0.009 micron",
"gra-ism01": "Graphite\n0.006 micron",
"gra-ism02": "Graphite\n0.009 micron"
}
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(10, 5))
denids = [f"n0{_}" for _ in range(5)]
for grain, ax in zip(select_grains, axes.flat):
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
m = CloudyModel(f'models/{prefix}')
Prad = m.data['flux']['F']/light_speed
R = m.data['rad']['radius']
Prad0 = L / (4*np.pi*light_speed*R**2)
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
Pgas[eden < 0.5*hden] = np.nan
Upsilon = Prad0/Pgas
Fbol = L / (4*np.pi*m.data['rad']['radius']**2)
G_n = m.data['flux']['F_F']/(F_habing*m.data['ovr']['hden'])
mm = m.data["gabun"][grain] >= 0.9*m.data["gabun"][grain].max()
msub = m.data['gdrift'][grain] < 20.0
if denid.endswith("00"):
label = f"{denid} {star}"
else:
label = f"{denid}"
ax.plot(Upsilon[mm & msub], m.data['gdrift'][grain][mm & msub],
alpha=0.85, lw=0.7, color=col, label=label)
ax.plot(Upsilon[mm & ~msub], m.data['gdrift'][grain][mm & ~msub],
alpha=0.85, lw=0.7, color=col, label="_nolabel_")
ax.axvline(300, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axvline(3000, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axhline(1.0, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.axhline(10.0, lw=0.5, ls='--', color='k', alpha=0.5, zorder=0)
ax.set_title(glabel[grain], fontsize="small", pad=-12)
axes[0, 0].legend(ncol=2, fontsize="xx-small", loc="left")
axes[1, 0].set(
yscale='log',
xscale='log',
xlabel=r'Radiation parameter: $\Xi = P_{\mathrm{rad}} / P_{\mathrm{gas}}$',
ylabel=r'$w_\mathrm{drift}$, km/s',
xlim=[1.0e-2, 4e4],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)- This is for comparison with the Kobulnicky papers, where they used Draine models that are meant for the ISRF, not O stars
- Plan:
- [X] Look at 70 micron emissivity versus radiation field U
- [X] Look at grain T vs U (but this will be different for different sizes)
- [X] Look at flux ratios between bands: 70/24, 24/8, 160/70
- For comparison with K17 fig 4 and 5
- We will define U as in K18, as F/F_0
- F = L / 4 π R^2
- F_0 = 0.0217 erg/s/cm^2 from Mathis:03a
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'grain-T-vs-U.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4,4))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
m = CloudyModel(f'models/{prefix}')
except:
continue
Prad = m.data['flux']['F']/light_speed
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
# Drop the zones where gas is neutral
U[eden < 0.5*hden] = np.nan
for grain in m.data['gtemp'].colnames[1::5]:
ls = '--' if grain.startswith('gra') else '-'
igrain = int(grain[-2:])
ax.plot(U, m.data['gtemp'][grain],
ls=ls,
alpha=0.6 - 0.05*igrain, lw=0.3 + 0.3*igrain, color=col)
ax.set(
yscale='log',
xscale='log',
xlabel=r'Radiation field: $U = L_*\, /\, (4 \pi \,R^{2} \, c \, u_{\mathrm{MMP83}})$',
ylabel=r'Grain Temperature, K',
#xlim=[1.0, 1.e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)- We divide by number density to put the emissivity as per nucleon
- Cloudy reports integral across frequency band, so we need to divide by filter width Δ\nu to put it in per Hz. From the following table, Δ\nu = 1.3405e12 Hz
Band λ_1 λ_2 ν_1 ν_2 Δ\nu IRAC 8 6.5 9.3 4.6122e13 3.2236e13 1.3886e13 MIPS 24 20.8 26.1 1.4413e13 1.1486e13 2.927e12 F 25 19 30 1.5779e13 9.9931e12 5.7859e12 PAC 70 60 82 4.9965e12 3.6560e12 1.3405e12 PAC 160 130 198 2.3061e12 1.5141e12 7.92e11 - [2019-04-16 Tue] Where do the numbers come from in the previous table?
- Chip says that Hazy has slightly different numbers
- He is right: in
continuum_bands.iniit has 60 → 82 for PAC 70, whereas originally I had 60 → 80 - All the others are fine
- [2019-04-16 Tue] Where do the numbers come from in the previous table?
- First, just look at the 70 micron emission
- Cloudy gives emissivity in erg/cm^3/s
- So we need to divide by 4π Δ\nu n 10^-23 to put it in Jy.cm^2/sr/H
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'grain-j70-vs-U.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4,4))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth70 = 1.3405e12
Jy = 1e-23
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
m = CloudyModel(f'models/{prefix}')
except:
continue
Prad = m.data['flux']['F']/light_speed
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
U[eden < 0.5*hden] = np.nan
j70 = m.data['emis']['PAC1 70.0000m'] / (4*np.pi*hden*bandwidth70*Jy)
ax.plot(U, j70, alpha=0.8, lw=1, color=col)
ax.plot(DL07tab['U'], DL07tab['70'], 'o', alpha=0.5, color='k')
ax.plot(DL07tab['U']/8.0, DL07tab['70'], 'o', alpha=0.2, color='k')
ax.set(
yscale='log',
xscale='log',
xlabel=r'Radiation field: $U = L_*\, /\, (4 \pi \,R^{2} \, c \, u_{\mathrm{MMP83}})$',
ylabel=r'Emissivity at 70 $\mu$m: $j_\nu$ [Jy cm$^{2}$ sr$^{-1}$ H$^{-1}$]',
#xlim=[1.0, 1.e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'grain-j24-vs-U.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4,4))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth24 = 2.927e12
Jy = 1e-23
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
m = CloudyModel(f'models-pah/{prefix}')
except:
continue
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
U[eden < 0.5*hden] = np.nan
j24 = m.data['emis']['MIPS 24.0000m'] / (4*np.pi*hden*bandwidth24*Jy)
ax.plot(U, j24, alpha=0.8, lw=1, color=col)
ax.plot(DL07tab['U'], DL07tab['24'], 'o', alpha=0.5, color='k')
ax.plot(DL07tab['U']/8.0, DL07tab['24'], 'o', alpha=0.1, color='k')
ax.set(
yscale='log',
xscale='log',
xlabel=r'Radiation field: $U = L_*\, /\, (4 \pi \,R^{2} \, c \, u_{\mathrm{MMP83}})$',
ylabel=r'Emissivity at 24 $\mu$m: $j_\nu$ [Jy cm$^{2}$ sr$^{-1}$ H$^{-1}$]',
#xlim=[1.0, 1.e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'grain-j8-vs-U.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4,4))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth8 = 1.3886e13
Jy = 1e-23
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
m = CloudyModel(f'models-pah/{prefix}')
except:
continue
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
U[eden < 0.5*hden] = np.nan
j8 = m.data['emis']['IRAC 8.00000m'] / (4*np.pi*hden*bandwidth8*Jy)
ax.plot(U, j8, alpha=0.8, lw=1, color=col)
ax.plot(DL07tab['U'], DL07tab['8'], 'o', alpha=0.5, color='k')
ax.plot(DL07tab['U']/8.0, DL07tab['8'], 'o', alpha=0.1, color='k')
ax.set(
yscale='log',
xscale='log',
xlabel=r'Radiation field: $U = L_*\, /\, (4 \pi \,R^{2} \, c \, u_{\mathrm{MMP83}})$',
ylabel=r'Emissivity at 8 $\mu$m: $j_\nu$ [Jy cm$^{2}$ sr$^{-1}$ H$^{-1}$]',
#xlim=[1.0, 1.e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)Now look at flux ratios between bands
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'grain-jratios-vs-U.pdf'
sns.set_color_codes("bright")
fig, axes = plt.subplots(3, 1, sharex=True, figsize=(6,9))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth70 = 1.3405e12
Jy = 1e-23
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
dnu8, dnu24, dnu70, dnu160 = 1.3886e13, 2.927e12, 1.2491e12, 7.92e11
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
m = CloudyModel(f'models/{prefix}')
except:
continue
Prad = m.data['flux']['F']/light_speed
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
U[eden < 0.5*hden] = np.nan
r24_8 = m.data['emis']['MIPS 24.0000m'] / m.data['emis']['IRAC 8.00000m']
r70_24 = m.data['emis']['PAC1 70.0000m'] / m.data['emis']['MIPS 24.0000m']
r160_70 = m.data['emis']['PAC3 160.000m'] / m.data['emis']['PAC1 70.0000m']
r24_8 *= dnu8/dnu24
r70_24 *= dnu24/dnu70
r160_70 *= dnu160/dnu70
axes[0].plot(U, r24_8, alpha=0.8, lw=1, color=col)
axes[1].plot(U, r70_24, alpha=0.8, lw=1, color=col)
axes[2].plot(U, r160_70, alpha=0.8, lw=1, color=col)
axes[0].plot(DL07tab['U']/8, DL07tab['24']/ DL07tab['8'], 'o', alpha=0.5)
axes[1].plot(DL07tab['U']/8, DL07tab['70']/ DL07tab['24'], 'o', alpha=0.5)
axes[2].plot(DL07tab['U']/8, DL07tab['170']/ DL07tab['70'], 'o', alpha=0.5)
#ax.plot(DL07tab['U']*0.1, DL07tab['70'], 'o', alpha=0.1)
for ax in axes:
ax.set(yscale='log', xscale='log')
axes[0].set(
ylabel=r'24 / 8',
)
axes[1].set(
ylabel=r'70 / 24',
)
axes[2].set(
xlabel='Radiation field: $U$',
ylabel=r'160 / 70',
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)And now the color-color plots
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'grain-color-color.pdf'
sns.set_color_codes("bright")
fig, axes = plt.subplots(2, 1, figsize=(4,8))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth70 = 1.3405e12
Jy = 1e-23
stars = [
["MS10", 0.63, 'Purples_d' ],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d' ],
["BSG", 30.2, 'Greens_d' ]
]
denids = [f"n0{_}" for _ in range(5)]
dnu8, dnu24, dnu70, dnu160 = 1.3886e13, 2.927e12, 1.2491e12, 7.92e11
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
try:
m = CloudyModel(f'models/{prefix}')
except:
continue
Prad = m.data['flux']['F']/light_speed
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
U[eden < 0.5*hden] = np.nan
r24_8 = m.data['emis']['MIPS 24.0000m'] / m.data['emis']['IRAC 8.00000m']
r70_24 = m.data['emis']['PAC1 70.0000m'] / m.data['emis']['MIPS 24.0000m']
r160_70 = m.data['emis']['PAC3 160.000m'] / m.data['emis']['PAC1 70.0000m']
r24_8 *= dnu8/dnu24
r70_24 *= dnu24/dnu70
r160_70 *= dnu160/dnu70
axes[0].plot(r24_8, r70_24, alpha=0.8, lw=1, color=col)
axes[1].plot(r70_24, r160_70, alpha=0.8, lw=1, color=col)
axes[0].plot(DL07tab['24']/ DL07tab['8'], DL07tab['70']/ DL07tab['24'], 'o', alpha=0.5)
axes[1].plot(DL07tab['70']/ DL07tab['24'], DL07tab['170']/ DL07tab['70'], 'o', alpha=0.5)
#ax.plot(DL07tab['U']*0.1, DL07tab['70'], 'o', alpha=0.1)
for ax in axes:
ax.set(yscale='log', xscale='log')
axes[0].set(
xlabel=r'24 / 8',
ylabel=r'70 / 24',
)
axes[1].set(
xlabel=r'70 / 24',
ylabel=r'160 / 70',
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)import glob
folders = glob.glob('DL07-data/U*')
header = [['U', '8', '24', '70', '110', '170'], None]
data = []
for folder in folders:
_, Ustring = folder.split('U')
U = float(Ustring)
datafile = f'U{Ustring}_{Ustring}_MW3.1_00.txt'
with open(folder + '/' + datafile) as f:
lines = f.readlines()
j70, j110, j170 = [float(_.split()[2]) for _ in lines if 'PAC' in _]
j24, _, _ = [float(_.split()[2]) for _ in lines if 'SST-MIPS' in _]
_, _, _, j8 = [float(_.split()[2]) for _ in lines if 'SST-IRAC' in _]
data.append([U, j8, j24, j70, j110, j170])
data.sort(key=lambda x: x[0])
data = header + data| U | 8 | 24 | 70 | 110 | 170 |
|---|---|---|---|---|---|
| 0.1 | 2.124e-18 | 8.856e-18 | 9.302e-17 | 3.393e-16 | 1.139e-15 |
| 0.15 | 3.184e-18 | 1.337e-17 | 1.731e-16 | 6.078e-16 | 1.736e-15 |
| 0.2 | 4.246e-18 | 1.793e-17 | 2.661e-16 | 8.962e-16 | 2.316e-15 |
| 0.3 | 6.375e-18 | 2.714e-17 | 4.845e-16 | 1.523e-15 | 3.432e-15 |
| 0.4 | 8.499e-18 | 3.654e-17 | 7.624e-16 | 2.208e-15 | 4.449e-15 |
| 0.5 | 1.062e-17 | 4.599e-17 | 1.052e-15 | 2.903e-15 | 5.448e-15 |
| 0.7 | 1.488e-17 | 6.538e-17 | 1.737e-15 | 4.346e-15 | 7.253e-15 |
| 0.8 | 1.701e-17 | 7.518e-17 | 2.104e-15 | 5.078e-15 | 8.112e-15 |
| 1.0 | 2.127e-17 | 9.479e-17 | 2.838e-15 | 6.544e-15 | 9.83e-15 |
| 1.2 | 2.555e-17 | 1.154e-16 | 3.72e-15 | 8.014e-15 | 1.125e-14 |
| 1.5 | 3.197e-17 | 1.463e-16 | 5.043e-15 | 1.022e-14 | 1.339e-14 |
| 2.0 | 4.265e-17 | 1.994e-16 | 7.383e-15 | 1.382e-14 | 1.663e-14 |
| 2.5 | 5.33e-17 | 2.545e-16 | 9.895e-15 | 1.734e-14 | 1.948e-14 |
| 3.0 | 6.395e-17 | 3.096e-16 | 1.241e-14 | 2.085e-14 | 2.233e-14 |
| 4.0 | 8.537e-17 | 4.283e-16 | 1.777e-14 | 2.739e-14 | 2.696e-14 |
| 5.0 | 1.068e-16 | 5.486e-16 | 2.32e-14 | 3.383e-14 | 3.137e-14 |
| 7.0 | 1.498e-16 | 8.091e-16 | 3.429e-14 | 4.566e-14 | 3.872e-14 |
| 8.0 | 1.713e-16 | 9.439e-16 | 3.988e-14 | 5.133e-14 | 4.205e-14 |
| 12.0 | 2.58e-16 | 1.524e-15 | 6.211e-14 | 7.231e-14 | 5.361e-14 |
| 15.0 | 3.233e-16 | 1.99e-15 | 7.867e-14 | 8.677e-14 | 6.097e-14 |
| 20.0 | 4.328e-16 | 2.83e-15 | 1.056e-13 | 1.09e-13 | 7.17e-14 |
| 25.0 | 5.431e-16 | 3.75e-15 | 1.318e-13 | 1.288e-13 | 8.053e-14 |
| 100.0 | 2.276e-15 | 2.28e-14 | 4.691e-13 | 3.404e-13 | 1.608e-13 |
| 300.0 | 7.62e-15 | 1.088e-13 | 1.118e-12 | 6.544e-13 | 2.562e-13 |
| 1000.0 | 3.299e-14 | 6.326e-13 | 2.606e-12 | 1.239e-12 | 4.081e-13 |
| 3000.0 | 1.433e-13 | 3.124e-12 | 5.127e-12 | 2.068e-12 | 5.974e-13 |
| 10000.0 | 7.118e-13 | 1.585e-11 | 1.004e-11 | 3.467e-12 | 8.887e-13 |
| 30000.0 | 3.018e-12 | 6.033e-11 | 1.736e-11 | 5.323e-12 | 1.249e-12 |
| 100000.0 | 1.663e-11 | 2.14e-10 | 3.021e-11 | 8.324e-12 | 1.803e-12 |
| 300000.0 | 9.378e-11 | 5.503e-10 | 4.696e-11 | 1.202e-11 | 2.463e-12 |
Write the above table to FITS format
from astropy.table import Table
tab = Table(rows=TAB[1:], names=TAB[0])
tab.write('DL07-data/emissivities.fits', overwrite=True)
- Models are calculated above in Cloudy models with the interstellar radiation field
- We want to compare directly with the Draine:2007a results
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'isrf-j70-vs-U.jpg'
sns.set_color_codes("bright")
fig, ax = plt.subplots(figsize=(4,4))
F_habing = 1.6e-3
F_mathis = 0.0217
light_speed = 2.99792458e10
kB = 1.3806503e-16
bandwidth70 = 1.3405e12
Jy = 1e-23
L = 1e6*3.82e33
# Data from Draine:2007a
DL07tab = Table.read('DL07-data/emissivities.fits')
# Read model
m = CloudyModel(f'models-isrf/dustrad-n00-isrf')
R = m.data['rad']['radius']
F = L / (4*np.pi*R**2)
U = F / F_mathis
hden = m.data['ovr']['hden']
j70 = m.data['emis']['PAC1 70.0000m'] / (4*np.pi*hden*bandwidth70*Jy)
ax.plot(U, j70, alpha=0.8, lw=1, color='r', label="Cloudy")
# Approximately correct for extinction
tau_V = m.data['ovr']['AV(point)'] * np.log(10.0)/2.5
U *= np.exp(tau_V)
ax.plot(U, j70, alpha=0.8, lw=1, color='r', ls='--', label="Extinction corrected")
ax.plot(DL07tab['U'], DL07tab['70'], 'o', alpha=0.5, color='k', label="DL07")
ax.legend()
ax.set(
yscale='log',
xscale='log',
xlabel=r'Radiation field: $U = L_*\, /\, (4 \pi \,R^{2} \, c \, u_{\mathrm{MMP83}})$',
ylabel=r'Emissivity at 70 $\mu$m: $j_\nu$ [Jy cm$^{2}$ sr$^{-1}$ H$^{-1}$]',
#xlim=[1.0, 1.e6],
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile, dpi=300)
fig.savefig(pltfile.replace('.jpg', '.pdf'))
PDF version: isrf-j70-vs-U.pdf
- ISRF SED is given in Appendix C of Mezger:1982a
- OB star SEDs come from Cloudy models
- [2019-04-15 Mon] Write SED to a file that can be used with Cloudy
| 8 | 4.33e-5 |
| 10 | 3.26e-5 |
| 12 | 2.21e-5 |
| 20 | 7.13e-6 |
| 25 | 1.14e-5 |
| 30 | 1.70e-5 |
| 40 | 3.06e-5 |
| 50 | 4.08e-5 |
| 60 | 4.76e-5 |
| 70 | 4.51e-5 |
| 80 | 4.20e-5 |
| 90 | 3.62e-5 |
| 100 | 3.16e-5 |
| 150 | 1.36e-5 |
| 200 | 5.61e-6 |
| 300 | 1.39e-6 |
| 400 | 5.61e-7 |
| 600 | 1.32e-7 |
| 800 | 4.20e-8 |
| 1000 | 2.14e-8 |
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.modeling.blackbody import blackbody_lambda
from astropy import units as u
from astropy.table import Table
from cloudytab import CloudyModel
def BB(wav, T):
"""Blackbody taking care of unit conversions
Input units: `wav` should be in micron, `T` should be in Kelvin
Output units: erg / (cm2 micron s sr)
"""
return blackbody_lambda(
wav*u.micron, T*u.K).to(
u.erg/(u.micron*u.s*u.cm**2*u.sr))
def mmp1983sed(wav):
"""Spectral Energy Distribution for interstellar radiation field
Appendix C of Mezger et al. 1982 with corrections from Appendix A
of Mathis et al. 1983 and additions from App B of same
"""
wav = np.atleast_1d(wav)
rslt = np.empty_like(wav)
# Masks for wavelength ranges
m1 = wav < 0.0912
m2 = (wav >= 0.0912) & (wav < 0.110)
m3 = (wav >= 0.110) & (wav < 0.134)
m4 = (wav >= 0.134) & (wav < 0.246)
m5 = (wav >= 0.246) & (wav < 8.0)
m6 = wav >= 7.2
# Piecewise functional forms for 4 pi J_lambda
rslt[m1] = 0.0
rslt[m2] = 38.57*wav[m2]**3.4172
rslt[m3] = 2.045e-2
rslt[m4] = 7.115e-4*wav[m4]**-1.6678
rslt[m5] = 4*np.pi*(
1.00e-14*BB(wav[m5], 7500) +
1.00e-13*BB(wav[m5], 4000) +
4.00e-13*BB(wav[m5], 3000)
)
# Multiply <8 mic part by wav to get SED
rslt[~m6] *= wav[~m6]
# While for > 8 mic part, we interpolate on the log SED not the F_lam
wtab, ftab = np.hsplit(np.array(TAB).astype('float'), 2)
logwav, logf = np.log10(wtab.ravel()), np.log10(ftab.ravel())
rslt[m6] = 10**np.interp(np.log10(wav[m6]), logwav, logwav + logf)
# Return the SED: lambda F_lambda
return rslt
pltfile = f'sed-comparison.pdf'
fig, ax = plt.subplots(figsize=(6, 3))
star = "MS20"
prefix = f"dustrad-n00-{star}"
m = CloudyModel(f'models/{prefix}')
nu = m.data['cont']['Cont nu']
nuFnu_inc = m.data['cont']['incident']
wav = 0.0912/nu
# Write ISRF SED to file: nu is already in Rydberg, but by default we
# want F_nu, so divide sed by nu (only shape is used by cloudy, so no
# need to worry about units)
isrf_Fnu = mmp1983sed(wav) / nu
mask = isrf_Fnu > 0.0
with open("models-isrf/isrf_mmp1983.sed", "w") as f:
f.write("\n".join(["\t".join(map(str, pair)) for pair in zip(nu[mask], isrf_Fnu[mask])]))
dilution = 1.0/(4*np.pi*u.pc.to(u.cm)**2)
ax.plot(wav, dilution*nuFnu_inc, label="O9.5V star @ 1 pc")
ax.plot(wav, 80*mmp1983sed(wav), label=r"IS radiation $\times 80$")
vmax = dilution*nuFnu_inc.max()
ax.legend()
ax.set(
xlim=[0.01, 1000.],
ylim=[1e-5*vmax, 2.*vmax],
yscale='log',
xscale='log',
xlabel='Wavelength, micron',
ylabel=r'$\lambda F_{\lambda}$, erg/s/cm$^2$',
)
sns.despine()
fig.tight_layout()
fig.savefig(pltfile)Sed table is written to models-isrf/isrf_mmp1983.sed
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'drift-tdust.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots()
F_habing = 1.6e-3
light_speed = 2.99792458e10
kB = 1.3806503e-16
stars = [
["MS10", 0.63, 'Reds_d'],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d'],
["BSG", 30.2, 'Greens_d']
]
denids = [f"n0{_}" for _ in range(5)]
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
m = CloudyModel(f'models/{prefix}')
Prad = m.data['flux']['F']/light_speed
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
for grain in m.data['gdrift'].colnames[1:]:
mm = m.data["gabun"][grain] >= 0.3*m.data["gabun"][grain].max()
ax.plot(m.data["gtemp"][grain][mm], m.data['gdrift'][grain][mm],
alpha=0.7, lw=0.5, color=col)
ax.set(
yscale='log',
xscale='log',
xlabel=r'$T_\mathrm{grain}$',
ylabel=r'$V_\mathrm{drift}$',
xlim=[10.0, 3000.0],
)
fig.savefig(pltfile)import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'drift-xH.pdf'
sns.set_color_codes("bright")
fig, ax = plt.subplots()
F_habing = 1.6e-3
light_speed = 2.99792458e10
kB = 1.3806503e-16
stars = [
["MS10", 0.63, 'Purples_d'],
["MS20", 5.45, 'Oranges_d'],
["MS40", 22.2, 'Blues_d'],
["BSG", 30.2, 'Greens_d']
]
denids = [f"n0{_}" for _ in range(5)]
for star, L4, cmap in stars:
L = 1e4*3.82e33*L4
cols = sns.color_palette(cmap, n_colors=len(denids))
for denid, col in zip(denids, cols):
# Read model
prefix = f"dustrad-{denid}-{star}"
m = CloudyModel(f'models/{prefix}')
Prad = m.data['flux']['F']/light_speed
hden = m.data['ovr']['hden']
eden = m.data['ovr']['eden']
Te = m.data['ovr']['Te']
Pgas = (hden + eden)*kB*Te
Fbol = L / (4*np.pi*m.data['rad']['radius']**2)
G_n = m.data['flux']['F_F']/(F_habing*m.data['ovr']['hden'])
# G_n = Fbol/(F_habing*m.data['ovr']['hden'])
x = m.data['ovr']['HII']
xx = x**2 / (1.0 - x)
for grain in m.data['gdrift'].colnames[1::2]:
mm = m.data["gabun"][grain] >= 0.9*m.data["gabun"][grain].max()
msub = m.data['gdrift'][grain] < 10.0
ax.plot(xx[mm & msub], m.data['gdrift'][grain][mm & msub], alpha=0.7, lw=0.5, color=col)
ax.plot(xx[mm & ~msub], m.data['gdrift'][grain][mm & ~msub], alpha=0.7, lw=0.5, color=col)
ax.set(
yscale='log',
xscale='log',
xlabel='$x^{2} / (1 - x)$',
ylabel=r'$V_\mathrm{drift}$, km/s',
xlim=[1e-3, 2.0e5],
)
fig.savefig(pltfile)from matplotlib import pyplot as plt
import seaborn as sns
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'models/{PREFIX}.png'
m = CloudyModel(f'models/{PREFIX}')
sns.set_palette("Reds_d", n_colors=10)
sns.set_color_codes()
fig, axes = plt.subplots(4, 1, sharex=True, figsize=(6, 9))
radius_pc = (m.data["rad"]["radius"])*u.cm.to(u.pc)
kT_eV = (m.data["ovr"]["Te"]*u.K*k_B).to(u.eV)
for gtype in m.data["gdrift"].colnames[1:]:
if gtype.startswith('sil'):
style = dict(lw=1.0, alpha=0.6, ls='-')
else:
style = dict(lw=0.6, alpha=1.0, ls='--')
axes[-1].plot(radius_pc, m.data["gdrift"][gtype], label=gtype, **style)
axes[-2].plot(radius_pc, m.data["gpot"][gtype], label=gtype, **style)
axes[-3].plot(radius_pc, m.data["gabun"][gtype], label=gtype, **style)
mm = m.data["gabun"][gtype] >= 0.3*m.data["gabun"][gtype].max()
axes[0].plot(radius_pc[mm], m.data["gtemp"][gtype][mm], label=gtype, **style)
axes[0].plot(radius_pc, m.data["ovr"]["Te"], color='g')
axes[-3].plot(radius_pc, m.data["gabun"]["total"], color='g', lw=1.5, label="_nolabel_")
abun_max = m.data["gabun"]["total"].max()
axes[-2].plot(radius_pc, kT_eV)
axes[-2].axhline(0.0, color='0.5', lw=0.5)
axes[-3].legend(ncol=2, loc="lower right", fontsize="xx-small")
axes[0].set(
xscale="log",
yscale="log",
ylabel="Temperature, K",
ylim=[0, None],
)
axes[1].set(
yscale="log",
ylabel="Acceleration",
)
axes[-3].set(
xscale="log",
yscale="log",
ylabel="Grain abundance",
ylim=[0.003*abun_max, 1.5*abun_max]
)
axes[-2].set(
xscale="log",
yscale="symlog",
#yticks=[-1, 0, 1, 10],
ylabel="Grain potential, eV",
)
axes[-1].set(
xscale="log",
yscale="log",
xlabel="Radius, pc",
ylabel="Drift velocity, km/s"
)
fig.suptitle(PREFIX, y=0.99)
sns.despine()
fig.tight_layout(h_pad=0.1)
fig.savefig(pltfile, dpi=300)file:
from matplotlib import pyplot as plt
import seaborn as sns
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'models/multi-dustprops.pdf'
selected_models = "n04-MS40", "n00-MS10", "n00-BSG"
sns.set_palette("Reds_d", n_colors=10)
sns.set_color_codes()
fig, axes = plt.subplots(4, 3, sharex="col", sharey="row", figsize=(12, 9))
for jcol, suffix in enumerate(selected_models):
m = CloudyModel(f'models/dustrad-{suffix}')
radius_pc = (m.data["rad"]["radius"])*u.cm.to(u.pc)
kT_eV = (m.data["ovr"]["Te"]*u.K*k_B).to(u.eV)
for gtype in m.data["gdrift"].colnames[1:]:
if gtype.startswith('sil'):
style = dict(lw=1.0, alpha=0.6, ls='-')
else:
style = dict(lw=0.6, alpha=1.0, ls='--')
axes[-3, jcol].plot(radius_pc, m.data["gabun"][gtype]/m.data["ovr"]["hden"],
label=gtype, **style)
mm = m.data["gabun"][gtype] >= 0.5*m.data["gabun"][gtype].max()
axes[0, jcol].plot(radius_pc[mm], m.data["gtemp"][gtype][mm], label=gtype, **style)
axes[-1, jcol].plot(radius_pc[mm], m.data["gdrift"][gtype][mm], label=gtype, **style)
axes[-2, jcol].plot(radius_pc[mm], m.data["gpot"][gtype][mm], label=gtype, **style)
axes[0, jcol].plot(radius_pc, m.data["ovr"]["Te"], color='g')
axes[-3, jcol].plot(radius_pc, m.data["gabun"]["total"]/m.data["ovr"]["hden"],
color='k', lw=1.5, label="_nolabel_")
abun_max = m.data["gabun"]["total"].max()
axes[-2, jcol].plot(radius_pc, kT_eV, color="g")
axes[-2, jcol].axhline(0.0, color='0.5', lw=0.5)
axes[0, jcol].set_title(suffix)
axes[0, 0].text(0.01, 6500, "Gas", color="g")
axes[0, 0].text(0.003, 200, "Grains", color="r")
axes[-2, 0].text(0.01, 0.4, "$k T$", color="g")
axes[-3, 1].legend(ncol=2, loc="center", fontsize="x-small", title="Grain type")
axes[0, 0].set(
xscale="log",
yscale="log",
ylabel="Temperature, K",
ylim=[15.0, 1.5e4],
)
axes[-3, 0].set(
xscale="log",
yscale="log",
ylabel="Grain abundance",
ylim=[3e-29, 2e-26]
)
axes[-2, 0].set(
xscale="log",
yscale="symlog",
yticks=[-1, 0, 1, 10],
yticklabels=["-1", "0", "1", "10"],
ylabel="Grain potential, eV",
ylim=[-3, 40.0]
)
axes[-1, 0].set(
yscale="log",
ylabel="Drift velocity, km/s"
)
for ax in axes[-1, :]:
ax.set(xscale="log", xlabel="Radius, pc")
sns.despine()
fig.tight_layout(h_pad=0.1)
fig.savefig(pltfile, dpi=300)- Saving every zone should be unnecessary
- we can just look at the incident and transmitted continuum for the last zone
- Since we don’t go into the PDR much, this should be representative of the FUV and EUV attenuation throughout the nebula
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'spectra-{PREFIX}.pdf'
m = CloudyModel(f'models/{PREFIX}')
fig, ax = plt.subplots()
nu = m.data['cont']['Cont nu']
nuFnu_inc = m.data['cont']['incident']
nuFnu_trans = m.data['cont']['trans']
nuFnu_tot = m.data['cont']['total']
vmax = nuFnu_inc.max()
ax.plot(nu, nuFnu_inc)
ax.plot(nu, nuFnu_tot, lw=0.6)
ax.plot(nu, nuFnu_trans, lw=0.3)
L_bol = np.trapz(nuFnu_inc/nu, nu)
mfuv = (nu >= 912.0/2000.0) & (nu <= 1.0)
meuv = (nu > 1.0) & (nu <= 4.0)
L_fuv = np.trapz(nuFnu_inc[mfuv]/nu[mfuv], nu[mfuv])
L_euv = np.trapz(nuFnu_inc[meuv]/nu[meuv], nu[meuv])
ax.axvspan(912.0/2000.0, 1.0, color='0.9')
ax.axvspan(1.0, 4.0, color='0.95')
ax.text(0.6, 2e-4*vmax, f"FUV\n{100*L_fuv/L_bol:.2f}%")
ax.text(1.5, 2e-4*vmax, f"EUV\n{100*L_euv/L_bol:.2f}%")
ax.set(
xlim=[0.05, 4.0],
ylim=[1e-6*vmax, 2.*vmax],
yscale='log',
xscale='log',
xlabel='Photon energy, Rydberg',
ylabel=r'$\nu F_{\nu}$',
)
fig.savefig(pltfile)- This might not be so useful with the current models since they go too deep for comparison with the bow shocks
- But it will be a start at least
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'wav-spectra-{PREFIX}.pdf'
m = CloudyModel(f'models/{PREFIX}')
fig, ax = plt.subplots()
LSUN = 3.82e33
nu = m.data['cont']['Cont nu']
wav = 0.0912/nu
nuFnu_inc = m.data['cont']['incident']/LSUN
nuFnu_trans = m.data['cont']['trans']/LSUN
nuFnu_diffuse = m.data['cont']['DiffOut']/LSUN
vmax = nuFnu_inc.max()
ax.plot(wav, nuFnu_inc, lw=0.2, label='incident')
ax.plot(wav, nuFnu_diffuse, lw=0.6, label='reflected')
ax.plot(wav, nuFnu_trans, lw=1.0, label='transmitted')
L_bol = np.trapz(nuFnu_inc/nu, nu)
mfuv = (nu >= 912.0/2000.0) & (nu <= 1.0)
meuv = (nu > 1.0) & (nu <= 4.0)
L_fuv = np.trapz(nuFnu_inc[mfuv]/nu[mfuv], nu[mfuv])
L_euv = np.trapz(nuFnu_inc[meuv]/nu[meuv], nu[meuv])
ax.legend()
ax.set(
xlim=[0.01, 1000.],
ylim=[1e-4*vmax, 2.*vmax],
yscale='log',
xscale='log',
xlabel='Wavelength, micron',
ylabel=r'$\nu L_{\nu}$, $\mathrm{L}_{\odot}$',
)
fig.savefig(pltfile)- I am testing saving this every zone so I can find FUV and EUV flux vs radius
- This makes an enormous file!
-rw-r--r-- 1 will staff 247M Mar 17 23:58 test-dust-tlusty.cont
- This makes an enormous file!
- Seems to be no way to get Cloudy to save flux in a band as function of depth
- For this example run (MS10 low density), the FUV is not absorbed at all because the ionization parameter is low, so dust τ is negligible
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from astropy.table import Table
from astropy import units as u
from astropy.constants import k_B
from cloudytab import CloudyModel
pltfile = f'fluxes-{PREFIX}.pdf'
m = CloudyModel(f'models/{PREFIX}')
# Split spectrum into a list of tables, one for each spatial zone
nzones = len(m.data['ovr'])
split_tabs = [Table(t) for t in np.split(m.data['cont'].as_array(), nzones)]
fig, ax = plt.subplots()
for tab in split_tabs[::10]:
ax.plot(tab['Cont nu'], tab['trans'])
ax.set(
xlim=[0.1, 10.0],
ylim=[1e25, 1e37],
yscale='log',
xscale='log',
xlabel='Photon energy, Rydberg',
ylabel='Flux',
)
fig.savefig(pltfile)(require 'generic-x) ;; we need this
(define-generic-mode
'cloudy-input-mode ;; name of the mode to create
'("#") ;; comments start with '#'
'("set" "stop" "hden" "table" "blackbody" "title" "element" "constant" "cmb"
"print" "save" "iterate" "time" "end" "cosmic ray" "coronal" "phi(h)"
"abundances" "luminosity" "table") ;; some keywords
'(("\\(#\\($\\| .*\\)\\)" 1 'font-lock-comment-face t)
("\\(//\\($\\| .*\\)\\)" 1 'font-lock-comment-face t)
("=" . 'font-lock-operator-face) ;; '=' is an operator
("\\b\\(scale\\|log\\|linear\\|file\\|units\\)\\b" . 'font-lock-constant-face)
("\\b\\(no\\|end\\|stop\\)\\b" 1 'font-lock-negation-char-face t)
("^title \\(.*\\)$" 1 'font-lock-doc-face t)
("^$" 1 'show-tabs-tab t)
("
\\(.*\\)$" 1 'font-lock-doc-face t)
) ;;
'("\\.in$") ;; files for which to activate this mode
nil ;; other functions to call
"A mode for Cloudy input files" ;; doc string for this mode
)