- docs/092.C-0323.pdf
- Evernote note with finding charts
Object |
---|
HH 514 |
HH 528 |
HH 523 |
HH 529-I |
HH 529-II/III |
HH 202-N |
HH 203 |
HH 204 |
- Discussions between Eduardo and Will during Eduardo’s visit to Morelia
- Gives red-shifted shoulder/wing to many lines
- More prominent in higher ionization (e.g., [O III]) but also seen in [N II]
- See the positive channels of these isovelocity maps:
- figs/multi-panel-isovel-oiii-nii.pdf
- Based on Doi:2004a and Garcia-Diaz:2008a
- This mainly maps the distribution of dust in the neutral/molecular filaments
- Same ones that are traced by fluorescent O I lines
- Secondary influence of orientation effects
- In the most red-shifted channels (Vhel > +50 km/s), face-on illuminated filaments predominate
- Edge-on illuminated filaments are seen at intermediate red-shifted velocities (Vhel ~ +40 km/s). Examples are:
- Bright Bar
- SW Compact Bar
This is from Garcia-Diaz:2007a
- Diagnostics of density/temperature from the O II permitted lines
- See Peimbert & Peimbert 2013ApJ…778…89P
- Also, see Will’s notes on analysis of this for nebula
- [X] Add those notes to project. Done
- Adal spectra
- Manu spectra
- OW Coordinate notation
- See Fig 3 of Doi:2004a
- Orion velocity maps
- Odell:1997b “High Velocity Features in the Orion Nebula”
- http://adsabs.harvard.edu/abs/1997AJ....114.2016O
- FP spectra of Huygens region, HH 203/4
- Doi:2004a “Internal Velocities in the Orion Nebula: Large Radial Velocity Features”
- http://adsabs.harvard.edu/abs/2004AJ....127.3456D
- Multiple longslit echelle spectra
- HH objects, Big Arc, DOH objects
- Henney:2007b “Large-Scale Flows from Orion South”
- Garcia-Diaz:2007a “Velocity Structure in the Orion Nebula. I. Spectral Mapping in Low-Ionization Lines”
- Garcia-Diaz:2008a “Velocity Structure in the Orion Nebula. II. Emission Line Atlas of Partially Ionized to Fully Ionized Gas”
- Odell:1997b “High Velocity Features in the Orion Nebula”
- Orion proper motions and imaging
- Hester:1991a
- Odell:1997a “Herbig Haro Objects in the Orion Nebula.”
- http://adsabs.harvard.edu/abs/1997AJ....114..730O
- HST imaging: HH 202, 202, 204/4
- Odell:2015a “THE NATURE AND FREQUENCY OF OUTFLOWS FROM STARS IN THE CENTRAL ORION NEBULA CLUSTER”
- Spectrophotometry
- Blagrave:2006a “A Photoionized Herbig-Haro Object in the Orion Nebula”
- http://adsabs.harvard.edu/abs/2006ApJ…644.1006B
- HH 529 spectra
- Mesa-Delgado:2009a “Properties of the ionized gas in HH 202. II: Results from echelle spectrophotometry with UVES”
- Espiritu:2017a “DUST DESTRUCTION AND THERMAL INHOMOGENEITIES IN THE ORION NEBULA. RESULTS OF HH 202.”
- Weilbacher:2015a “A MUSE map of the central Orion Nebula (M 42)”
- McLeod:2016a “A nebular analysis of the central Orion nebula with MUSE”
- Blagrave:2006a “A Photoionized Herbig-Haro Object in the Orion Nebula”
- Infrared imaging
- Smith:2005a “Thermal Dust Emission from Proplyds, Unresolved Disks, and Shocks in the Orion Nebula”
- http://adsabs.harvard.edu/abs/2005AJ....130.1763S
- 11.7 micron emission from HH 529 I and III, and jet
- Smith:2005a “Thermal Dust Emission from Proplyds, Unresolved Disks, and Shocks in the Orion Nebula”
- Neutral gas
- van-der-Werf:2013a
- Atomic physics
- Laha:2017a “Emission Line Ratios of FE III as Astrophysical Plasma Diagnostics”
Muy buenas preguntas las dos!
¿por qué 150 km/s?. En un choque con un átomo de H a 100 km/s hay 52.2 eV. Las colisiones de O+ con electrones tienen poco efecto porque la masa es muy pequeña, pero ¿qué pasa con los nucleos?
No es la ionización colisional que nos importa aquí. Sí tienes razón que un choque de 100 km/s produce O++ por colisiones - por eso hay la posibilidad de [O III] de las zonas dos (ionización) y tres (enfriamiento). Pero para el precursor la pregunta es sí la radiación EUV emitida por el gas poschoque es capaz de ionizar O++ en el gas antes del choque. Es eso que requiere velocidades mayores
Si T([OIII]) ya está en equilibrio, ¿porqué es más caliente que T([NII])?
La verdad, eso no entiendo. Tienes razón que es un argumento en contra de lo que propongo… Voy a pensar más sobre esto.
Lástima que la rendija no fue más larga!
Olvidaba mencionar sobre la figura 1 – hice una versión más simple para ir en la introducción si están de acuerdo. Luego la versión con todos los contornos, etc puede ir en la discusión. ¿Les parece bien?
Sí, lo de la T([OIII]) tiene una interpretación compleja. Le he dado unas vueltas. Basado en lo mismo que has encontrado, si en las imágenes 2D del HST vemos que [OIII]/Ha es más alto en una capa periférica delgada de HH204 por la emisión en la zona de enfriamiento esto también debería suceder en zonas equivalentes de la geometría 3D. Entonces deberíamos tener dos capas con [OIII]/Ha altas integradas en la línea de visión de la rendija. Es decir, que la forma paraboloidal de HH204 debería estar recubierta por una manta de gas con [OIII]/Ha alto, por lo que nuestra rendija integraría dos veces esa manta (por delante y por detrás de HH204). El grosor de la manta sería variable -dependiendo de la densidad del gas- en el bowshock sería casi cero, como lo mostraste en la figura 2c del PDF que nos mandaste, mientras que en el jet, con densidades algo menores, dicha capa es más gruesa (y aporta más emisión). Esto podría ser complementario a lo que comentas, que en el espectro de UVES de HH204 hay emisión principalmente del bowshock y del Mach disk. En la zona del bowshock, la zona de enfriamiento no aporta casi nada porque es la zona más densa y es ahí (en una distancia de menos de 5 mpc) donde se observa un aumento importante de la T([OIII]) que está fuera de equilibrio. En la zona del mach disk, ya debería estar enfriada y en equilibrio y a la misma temperatura que T([NII]) pero en esa zona, la capa de la zona de enfriamiento debería ser más grande, porque la densidad es algo menor y tal vez esa aportación es lo que sube la temperatura un poco. ¿Qué piensas sobre eso?
Encontré un artículo de Raymond sobre el Bucle de Cisne que sería util para citar. Aunque hay unas diferencias con nuestro caso (en particular, densidad e ionización) es suficiente similar para usar como analogía. Esta figura muestra el caso de una velocidad un poco más alta que la nuestra en donde se aprecia la diferencia entre la longitud de enfriamiento y la anchura de la zona de emissión [O III]
https://ui.adsabs.harvard.edu/abs/2020ApJ…894..108R
- From Fig 6 of the paper, we have density increasing from 1e4 to 2e4 as one approaches the shock
- For a working surface, we predict a jump of M^2 at the shock, and a total change of M^2 exp(0.5/M^2) up to the CD
M M^2 M^2 exp(0.5/M^2) 0.0 1.00 1.00 0.2 1.00 1.02 0.4 1.00 1.08 0.6 1.00 1.20 0.8 1.00 1.38 1.0 1.00 1.65 1.1 1.21 1.83 1.15 1.32 1.93 1.18 1.39 1.99 1.2 1.44 2.04 1.3 1.69 2.27 1.4 1.96 2.53
- We have already calibrate these images I think
- Yes, they are in hh204/hh204-pc-1994-Sha.fits and similar
- Brightness measurements:
- BG of nose (outside shock): 0.0080 +/- 0.0003
- Peak at nose: 0.022 +/- 0.001
- BG at N wing (outside shock): 0.0076 +/- 0.0002
- Peak at near N wing: 0.0150 +/- 0.0003
- Peak at medium N wing: 0.0112 +/- 0.0004
- Peak at far N wing: 0.0107 +/- 0.0003
- Jet clumps: 0.028
- Jet BG (gaps between clumps): 0.0190 +/- 0.0007
- Bow shock filament in S wing: 0.0108 +/- 0.0002
- Filament BG (still inside bow): 0.0091 +/- 0.0002
- Possible calibration error
- Was the O’Dell calibration for WFC pixels?
- If so, then brightness needs to be multiplied by (0.0996/0.046)**2 = 4.688
- See Finally make ratios 502/656 and 658/656
- It looks like this is indeed the case
- Calculations
- S(Ha) = (α E / 4 pi) n^2 z
- α = 8.6052e-14 cm^3 s^-1 at 10^4 K
- E(Hα) = 3.027e-12 erg
- Some sizes:
- Axial length of “plug”: 7.4 arcsec
- Cylindrical radius of near N wing: 2.9 arcsec
- Thickness of near N wing: 0.3 arcsec
- Cylindrical radius of medium N wing: 4.2 arcsec
- Thickness of medium N wing: 0.5 arcsec
- Cylindrical radius of far N wing: 6.1 arcsec
- Thickness of far N wing: 0.7 arcsec
- Radius of curvature of nose: 3.2 arcsec
- Thickness of nose shell: 0.4 arcsec
- Cylindrical radius of S wing filament: 3.4 arcsec
- Thickness of S wing filament: 0.4 arcsec
- Jet knot diameter: 0.6 arcsec
- Conversion to LOS thickness, dz:
- dz = 2 sqrt(R h)
- Nose: dz = 2 sqrt(3.2 0.4) = 2.3 arcsec
- Near N wing: dz = 2 sqrt(2.9 0.3) = 1.9 arcsec
- Medium N wing: dz = 2 sqrt(4.2 .5) = 2.9 arcsec
- Far N wing: dz = 1
Nose 0.022 0.0080 2.3 2.07e4 Near N wing 0.0150 0.0076 1.9 1.66e4 Med N wing 0.0112 0.0076 2.9 9.35e3 Far N wing 0.0107 0.0076 4.1 7.30e3 S filament 0.0108 0.0091 2.3 7.21e3 Jet knot 0.028 0.0190 0.6 3.25e4 Average plug 0.022 0.0080 5 1.40e4
- Eduardo finds that [O III]/Hβ for the blue component has a U-shaped profile
- Starting at the bottom of the slit (+12 mpc), it is 4959/4861 = 0.08
- Flat minimum of 0.03 to 0.04 from +2 to +10 mpc
- Rises again to 0.06 at +0 mpc
- The initial fall would seem to indicate at ionization gradient due to EUV opacity
- But what causes the increase close to the bow shock?
- And how can we reconcile this with the HST images?
- These imply a roughly constant value - see [[id:A4B5B542-E975-46AC-9A30-93DDC6D50403][Implications for the rim increase in the [O III] / Ha ratio]]
- But they only go to 10 mpc - maybe I should redo it
- Part of the answer may be that Eduardo’s methodology selects only the blue-shifted emission,
- This will predominantly be the jet material, at least towards the bottom of the slit
- At the top, near the bow shock it will include some of the bow shock shell, although some of that will also be at lower velocities
- Eduardo’s fig 8 shows that the O++ abundance also has the same U-shaped profile
pos O+ O++ frac++ +12.2 8.55 7.2 0.045 +10 8.55 6.8 0.018 +5 8.60 6.65 0.011 +0 8.65 7.0 0.022 - This seems to track the [O III]/Ha, which is odd given the T increase between +5 and +0
- Fig 6 shows it goes from 12000 to 17000 K
- This should effect the emissivity ratio quite a lot. From pyneb I find a ratio of 2.5
- This is larger than the increase in 4959/4861, which implies that O++/O keeps going down from +5 mpc to +0 mpc
- Check this with a table
pos T, K n, pcc e4959/4861 I4959/4861 12+log(O++/H+) (1) (2) (3) (4) (5) (6) 0.0 17000 20000 42254 0.060 6.15 2.4 15000 17000 31360 0.035 6.05 4.9 12500 14000 19287 0.030 6.19 7.3 12000 14000 17124 0.035 6.31 9.8 12000 11000 17193 0.045 6.42 12.2 12000 10000 17215 0.080 6.67 - Columns:
- (1), (2), (3) from Figure 6: temperature and density for blue-shifted component as a function of position along slit, in mpc with origin at the shock
- (4) Ratio in per-ion emission coefficients of [O III] 4959 to H beta, calculated with PyNeb using columns (2) and (3).
- (5) Observed intensity ratio [O III] 4959 to H beta from Figure 7
- (6) Calculated ionic abundance 12+log(O++/H+) from columns (4) and (5)
- I4959/4861 = e4959/4861 O++/H+
- => O++/H+ = I4959/4861 / e4959/4861
- Columns:
- We see an ionization gradient in the blue-shifted component
MUSE Line | BG | 1 | 2 | 3 | 4 |
---|---|---|---|---|---|
[O I] 6300 | 10000 | 55000 | 82000 | 40000 | 12000 |
[O II] 7319+20 | 40000 | 170000 | 340000 | 310000 | 200000 |
[O III] 5007 | 1400000 | 1600000 | 1800000 | 2100000 | 2200000 |
6300 / 7319+20 | 0.250 | 0.346 | 0.240 | 0.111 | 0.013 |
7319+20 / 5007 | 0.029 | 0.650 | 0.750 | 0.386 | 0.200 |
import numpy as np
import pyneb as pn
o1 = pn.Atom(atom="O1")
o2 = pn.Atom(atom="O2")
o3 = pn.Atom(atom="O3")
tem = 10000.0
den = 2 .0e4
e6300 = o1.getEmissivity(tem, den, wave=6300)
e7319 = o2.getEmissivity(tem, den, wave=7319)
e7320 = o2.getEmissivity(tem, den, wave=7320)
e5007 = o3.getEmissivity(tem, den, wave=5007)
print(f"e(6300) / e(7319+20) = {e6300 / (e7319 + e7320):.3f}")
print(f"e(7319+20) / e(5007) = {(e7319 + e7320) / e5007:.3f}")
- file:~/Google Drive/Colab Notebooks/HH204 neutral oxygen fraction.ipynb
import numpy as np
import pyneb
REDCORR = pyneb.extinction.red_corr.RedCorr(
law='CCM89 Bal07', R_V=5.5, cHbeta=1.0
)
flam = np.log10(REDCORR.getCorrHb(6563))
print("1 + f_lam =", 1 + flam)
- For HH 204, C(Hb) = 0.37 +/- 0.04 from MUSE map
- ~/Dropbox/dib-scatter-hii/data/orion-muse/C_H_beta-bin001.fits
- So correction for Ha is 10**(0.37 0.78) = 1.94
- So this explains the factor of 1.9 used above in Estimate density from Hα brightness?
- [O III] is in red
- Hα is in cyan
- You can see a very thin [O III] layer at the leading edge of the shock
- But why is the [O III] emission not limb-brightened?
- The [O III] is not limb-brightened, possibly because the emissivity in the cooling zone is not significantly higher than in the shell. Or possibly it is, but over such a small distance that it isn’t even resolved with HST
- Wrong The [O III] is limb-brightened
- It is just very weak, and then superimposed on a stronger background from the equilibrium shell
- Also, the width of the
- This is based on what I did years ago in
- The PC images are described in
- They are the
fov2f{502,547,656,658}.fits
files- Each is 4 HDUs, ones for each WFPC2 chip, starting with PC
import pyds9
acs = pyds9.DS9("ACS")
acs.set("cd /Users/will/Work/BobPC/2005")
<<ds9-setup>>
acs.set("rgb")
for item in "bin|colorbar|slice|smooth".split("|"):
acs.set("rgb lock %s yes" % (item))
for chan, filter_, bmin, bmax in [
["blue", 658, 30, 450],
["green", 656, 50, 400],
["red", 502, 7, 35],
]:
acs.set("rgb %s" % (chan))
acs.set("file fov2f%i.fits" % (filter_))
acs.set("scale limits %i %i" % (bmin, bmax))
Save the WCS files
<<ds9-setup>>
rslts = [["Filter", "Limits"], None]
for rgb, f in [["red", "f502"], ["green", "f656"], ["blue", "f658"]]:
acs.set(f"rgb {rgb}")
acs.set(f"wcs save fov2{f}-align-Robberto.wcs")
limits = acs.get("scale limits")
rslts.append([f, limits])
Filter | Limits |
---|---|
f502 | 7 35 |
f656 | 50 400 |
f658 | 30 450 |
- We use the star NIC 1211
- RA, DEC : 83.8433516 -5.4192935
- F502 x, y : 354.7 439.9
- F656 x, y : 355.0 440.1
- F658 x, y : 355.0 439.5
- F547 x, y : 354.8 439.9
- FITS structure is
0 PRIMARY 1 PrimaryHDU 266 (800, 800, 4) float32 1 s2o_cvt.tab 1 TableHDU 353 4R x 49C
- Where the table seems to be a list of BG brightness levels at various positions
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.table import Table
pcdir = Path("../../Work/BobPC/2005")
outdir = Path("hh204")
for filt in "f502", "f656", "f658", "f547":
hdulist = fits.open(pcdir / f"fov2{filt}.fits")
w2d = WCS(hdulist[0].header).celestial
fits.PrimaryHDU(
header=w2d.to_header(),
data=hdulist[0].data[0, :, :]
).writeto(outdir / f"fov2{filt}-PC.fits", overwrite=True)
# print(hdulist.info())
# print(Table.read(hdulist[1]))
- [2021-02-08 Mon] Change to use my new wcsfile package, which is developed with nbdev
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from reproject import reproject_interp
import wcsfile
datadir = Path("hh204")
hdu502 = fits.open(datadir / "fov2f502-PC.fits")[0]
hdu502.header.update(**wcsfile.read(datadir / "fov2f502-align-Robberto-PC.wcs"))
hdu502.writeto(datadir / "hh204-pc-1994-align-Robberto-f502n.fits", overwrite=True)
hdu656 = fits.open(datadir / "fov2f656-PC.fits")[0]
hdu656.header.update(**wcsfile.read(datadir / "fov2f656-align-Robberto-PC.wcs"))
newdata, _ = reproject_interp(hdu656, hdu502.header)
hdu656.data = newdata
hdu656.header = hdu502.header
hdu656.writeto(datadir / "hh204-pc-1994-align-Robberto-f656n.fits", overwrite=True)
hdu658 = fits.open(datadir / "fov2f658-PC.fits")[0]
hdu658.header.update(**wcsfile.read(datadir / "fov2f658-align-Robberto-PC.wcs"))
newdata, _ = reproject_interp(hdu658, hdu502.header)
hdu658.data = newdata
hdu658.header = hdu502.header
hdu658.writeto(datadir / "hh204-pc-1994-align-Robberto-f658n.fits", overwrite=True)
hdu547 = fits.open(datadir / "fov2f547-PC.fits")[0]
hdu547.header.update(**wcsfile.read(datadir / "fov2f547-align-Robberto-PC.wcs"))
newdata, _ = reproject_interp(hdu547, hdu502.header)
hdu547.data = newdata
hdu547.header = hdu502.header
hdu547.writeto(datadir / "hh204-pc-1994-align-Robberto-f547m.fits", overwrite=True)
- I want to put these in physical units, so I can compare with the MUSE spectra
- ODell:2009b defines various coefficients for each filter
- Primary coefficient K1
- Line contamination coefficient K2
- Continuum contamination coefficient K3
- Surface brightness S = (R / K1) (1 - K2 R_other / R - K2 R_cont / R)
- R, R_other, R_cont are count rates
- ODell:2009b Table 3: Recommended values for K1
F469N 0.51e-10 F487N 0.70e-10 F502N 0.82e-10 F656N 1.62e-10 F658N 1.60e-10 F673N 1.42e-10
- This gives S in units of phot/s/cm^2/sr
- So need to multiply by h c / λ to get in energy units
- 6.62606876e-27 2.99792458e10 / 1e-8 = 1.98644544044e-8
import numpy as np
from astropy.io import fits
h_c_over_AA = 1.98644544044e-8
hdu502 = fits.open("hh204-pc-1994-align-Robberto-f502n.fits")[0]
hdu656 = fits.open("hh204-pc-1994-align-Robberto-f656n.fits")[0]
hdu658 = fits.open("hh204-pc-1994-align-Robberto-f658n.fits")[0]
hdu547 = fits.open("hh204-pc-1994-align-Robberto-f547m.fits")[0]
# convert to counts/s/pixel
raw656 = hdu656.data / 200.0
raw658 = hdu658.data / 500.0
raw502 = hdu502.data / 200.0
raw547 = hdu547.data / 50.0
# Use O'Dell 2009 eq 1
S_ha = raw656 * (1.0 - 0.148 * raw658 / raw656 - 0.051 * raw547 / raw656)
S_ha /= 1.62e-10
S_ha *= h_c_over_AA / 6563
S_nii = raw658 * (1.0 - 0.038 * raw656 / raw658 - 0.062 * raw547 / raw658)
S_nii /= 1.60e-10
S_nii *= h_c_over_AA / 6583
S_oiii = raw502 * (1.0 - 0.028 * raw547 / raw502)
S_oiii /= 0.82e-10
S_oiii *= h_c_over_AA / 5007
# Save ratios
fits.PrimaryHDU(
header=hdu502.header,
data=S_oiii / S_ha,
).writeto("hh204-pc-1994-ratio-Oiii-Ha.fits", overwrite=True)
fits.PrimaryHDU(
header=hdu658.header,
data=S_nii / S_ha,
).writeto("hh204-pc-1994-ratio-Nii-Ha.fits", overwrite=True)
fits.PrimaryHDU(
header=hdu502.header,
data=S_oiii / raw547,
).writeto("hh204-pc-1994-ratio-Oiii-cont.fits", overwrite=True)
# And save calibrated surface brightness in erg/s/cm2/sr
fits.PrimaryHDU(
header=hdu656.header,
data=S_ha,
).writeto("hh204-pc-1994-Sha.fits", overwrite=True)
fits.PrimaryHDU(
header=hdu658.header,
data=S_nii,
).writeto("hh204-pc-1994-Snii.fits", overwrite=True)
fits.PrimaryHDU(
header=hdu502.header,
data=S_oiii,
).writeto("hh204-pc-1994-Soiii.fits", overwrite=True)
Check the answers (make new frame by hand)
import pyds9
acs = pyds9.DS9("ACS")
acs.set("cd /Users/will/Dropbox/Orion-HH/hh204")
acs.set("file hh204-pc-1994-ratio-Oiii-Ha.fits")
- Make a figure of the ratio
import numpy as np
import matplotlib.pyplot as plt
from astropy.wcs import WCS
from astropy.io import fits
import regions as rg
import seaborn as sns
figfile = f"hh204-ratio-oiii-ha-HST-PC-1994.pdf"
hdu = fits.open("hh204/hh204-pc-1994-ratio-Oiii-Ha.fits")[0]
w = WCS(hdu.header)
fig, ax = plt.subplots(figsize=(5, 3.75), subplot_kw=dict(projection=w))
sns.set_color_codes("bright")
m = (hdu.data < 0.0) | (hdu.data > 1.0)
hdu.data[m] = 0.335 # median in box to left of HH object
im = ax.imshow(hdu.data, vmin=0.1, vmax=0.55, cmap="gray")
cb = fig.colorbar(im, ax=ax, label="[O III] 5007 / Hα 6563")
regions = rg.read_ds9(f"hh204-slit.reg")
for skyregion in regions:
pixelregion = skyregion.to_pixel(w)
pixelregion.plot(ax=ax, color="w", lw=1.4, ls="--", alpha=0.7)
ax.set(
xlabel="RA (J2000)",
ylabel="Dec (J2000)",
xlim=[100, 550],
ylim=[80, 530],
)
fig.tight_layout(rect=[0.18, 0.05, 1.0, 1.0])
fig.savefig(figfile)
- Illustrative values of ratio (now corrected to put them properly in energy units)
Original Corrected N wing, outside shock 0.195 +/- 0.02 0.34 +/- 0.03 N wing, shock rim 0.21 +/- 0.01 0.36 +/- 0.02 N wing, inside 0.175 +/- 0.02 0.30 +/- 0.03 Slit, outside shock 0.23 +/- 0.02 0.40 +/- 0.03 Slit, inside shock 0.14 +/- 0.02 0.24 +/- 0.03 Slit, further inside shock 0.11 +/- 0.02 0.19 +/- 0.03
- We are measuring 5007/6563
- In the UVES spectra we have 4959/4861
- So, 4959/4861 = (4959/5007) (6563/4861) (5007/6563)
- Where 4959/5007 = 1/3
- And 6563/4861 can come from the MUSE maps: 3.3 upstream of shock
- conversion factor 4959/4861 = 1.1 (5007/6563)
- We now have agreement between MUSE and HST
- MUSE Slit, outside shock: 5007/4861 = 1.30 => 5007/6563 = 1.3/3.3 = 0.39 - close enough!
- So this implies that 4959/4861 = 1.30 / 3 = 0.433
- This is a bit lower than the 0.6 that I get below from UVES
First we subtract continuum, then we subtract upstream BG
[O III] 4363Region | 4363 cont | 4363 fast | 4363 slow | 4363 inter | fast - cont | slow - cont | inter - cont |
---|---|---|---|---|---|---|---|
BG 1 | 10 | 12 | 32 | 15 | 2 | 22 | 5 |
Shock 2 | 13 | 27 | 35 | 22 | 14 | 22 | 9 |
Shock 3 | 15 | 28 | 37 | 20 | 13 | 22 | 5 |
Shock 4 | 15 | 32 | 37 | 22 | 17 | 22 | 7 |
Region | 4959 cont | 4959 fast | 4959 slow | 4959 inter | fast - cont | slow - cont | inter - cont |
---|---|---|---|---|---|---|---|
BG 1 | 8 | 40 | 2400 | 90 | 32 | 2392 | 82 |
Shock 2 | 10 | 200 | 2400 | 150 | 190 | 2390 | 140 |
Shock 3 | 12 | 400 | 2700 | 250 | 388 | 2688 | 238 |
Shock 4 | 11 | 500 | 2800 | 300 | 489 | 2789 | 289 |
Region | 4861 cont | 4861 fast | 4861 slow | 4861 inter | fast - cont | slow - cont | inter - cont |
---|---|---|---|---|---|---|---|
BG 1 | 12 | 100 | 4000 | 2500 | 88 | 3988 | 2488 |
Shock 2 | 15 | 2500 | 4700 | 3500 | 2485 | 4685 | 3485 |
Shock 3 | 20 | 6500 | 4600 | 4500 | 6480 | 4580 | 4480 |
Shock 4 | 15 | 5500 | 4600 | 4200 | 5485 | 4585 | 4185 |
Region | 4363 | 4959 | 4861 | 4363/4959 | 4959/4861 |
---|---|---|---|---|---|
BG 1 | 22 | 2392 | 3988 | 0.009 | 0.600 |
Sh 2 - BG | 0 | 0 | 700 | 0/0 | 0 |
Sh 3 - BG | 0 | 500 | 600 | 0 | 0.833 |
Sh 4 - BG | 0 | 400 | 600 | 0 | 0.667 |
- The fact that 4959/4861 is similar in the BG and BG-subtracted regions, suggests that this is just variation in the nebula
- It is odd that 4959/4861 is so large, given that we calculate it should be 0.433 in the upstream regions, based on the MUSE and HST photometry
- Maybe the problem is the extra broadening of H beta
Region | 4363 | 4959 | 4861 | 4363/4959 | 4959/4861 | 30 x 4363/4861 |
---|---|---|---|---|---|---|
BG 1 | 2 | 32 | 88 | 0.063 | 0.364 | 0.682 |
Sh 2 - BG | 12 | 158 | 2400 | 0.076 | 0.066 | 0.150 |
Sh 3 - BG | 11 | 356 | 6400 | 0.031 | 0.056 | 0.052 |
Sh 4 - BG | 15 | 457 | 5400 | 0.033 | 0.085 | 0.083 |
- So factor of 2.5 decrease in 4363/4959, going from close behind shock to further away
- And then another factor of 3 going to the nebula component
- Comparing with Fig 7 of paper gives similar results for 30 x 4363/4861
- And for 4959/4861 too
- This time, I do not subtract the upstream BG from the shock regions, since I am interested in the fraction of all the emission in each line that comes from the blue-shifted component
- This is for comparison with the HST images, which cannot separate the two.
- Final 3 columns give the fractions:
- f = B / (R + B)
- For each line, it is the fraction of the total intensity that comes from the blueshifted component
- That is, the fraction emitted by high-velocity gas
Region | 4363B | 4363R | 4959B | 4959R | 4861B | 4861R | f4363 | f4959 | f4861 |
---|---|---|---|---|---|---|---|---|---|
BG 1 | 2 | 22 | 32 | 2392 | 88 | 3988 | 0.083 | 0.013 | 0.022 |
Sh 2 | 12 | 22 | 158 | 2390 | 2400 | 4685 | 0.353 | 0.062 | 0.339 |
Sh 3 | 11 | 22 | 356 | 2688 | 6400 | 4580 | 0.333 | 0.117 | 0.583 |
Sh 4 | 15 | 22 | 457 | 2789 | 5400 | 4585 | 0.405 | 0.141 | 0.541 |
- Results:
- For [O III] 4363, the blue fraction is roughly constant at 0.35 to 0.4 for all of the bow shock region
- For [O III] 4959, the blue fraction rises from 0.06 just behind the shock, up to 0.14 at the end of the slit
- For comparison, the HST image of BG-subtracted Soiii goes from 400 in the upstream region, to 460 just behind the shock: increase of 1.15
- And using a sample box on the slit spectra, it goes from 181 - 10 = 171 upstream, to 208 - 10 = 198 behind shock: increase of 1.16
- **Inconsistency** blue fraction is 0.06 but the previous suggests should be 0.15 Why?
- It turns out that the shock has emission to the red as well. If we add this in, the fraction gets to be nearer to 0.15. Some of this may be scattered, but maybe not
- Note that in 4959, the red component is just as strong as the blue component, whereas in 4959 it is hardly detectable.
- This makes me think it might be related more to the nebular component than to the shock.
- But in that case, why does it increase in strength at the position of the shock?
- It is also very weak, but detectable, in [O II] 3726
- For H beta 4861, the fraction rises and falls slightly in the range 0.33 to 0.58
- The above considerations imply that the bow shock contribution to 4959,5007 is of order 10%
- We can measure this along the upper wing where we see the narrow ridge of increased oiii/ha
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
datadir = Path("hh204")
x, Soiii = np.loadtxt(datadir / f"{ID}-Soiii.dat", unpack=True)
_, Sha = np.loadtxt(datadir / f"{ID}-Sha.dat", unpack=True)
# Convert to mpc
x *= 1.94 * 0.045
bg_oiii = Soiii[:IBG].mean()
bg_ha = Sha[:IBG].mean()
ishock = np.argmax(np.gradient(Sha))
x -= x[ishock]
RR = Soiii / Sha # Raw ratio
R_bg = bg_oiii / bg_ha # Nebular background value
R = (Soiii - bg_oiii) / (Sha - bg_ha) # Nebular-corrected value
figfile = f"{ID}-oiii-ha-ratio.pdf"
fig, [ax2, ax] = plt.subplots(2, 1, sharex=True, figsize=(3, 3.75))
ax.fill_between(x[IBG:], R[IBG:], step="mid", color="0.7")
ax.axhline(R_bg, linestyle="--", color="k", alpha=0.5)
ax.axhline(0.0, linestyle="-", color="k", alpha=1.0, lw=0.5)
ax.axvline(0.0, linestyle="-", color="k", alpha=1.0, lw=0.5)
ax.plot(x, RR[:], color="k")
ax.set(ylim=[-0.05, 0.7], xlabel="Displacement, mpc", ylabel="[O III] / Ha")
ax2.axvline(0.0, linestyle="-", color="k", alpha=1.0, lw=0.5)
ax2.axhline(1.0, linestyle="--", color="k", alpha=0.5)
ax2.plot(x, Sha/bg_ha, label=r"H$\alpha$")
ax2.plot(x, Soiii/bg_oiii, label=r"[O III]")
ax2.set(ylim=[0.0, None], ylabel="Relative brightness")
ax2.legend()
sns.despine()
fig.tight_layout()
fig.savefig(figfile)
import numpy as np
from astropy.io import fits
for night in "dic1", "dic2":
fname = f"hh204/blue_2D_{night}.fits"
hdu, = fits.open(fname)
# avoid the top row
BG = np.mean(hdu.data[-6:-4, :], axis=0, keepdims=True)
hdu.data -= BG
hdu.writeto(fname.replace(".fits", "_bgsub.fits"), overwrite=True)
Description and discussion of these spectra is documented here
- Craft document: HH 204 bow shock UVES profiles
- Public link to Craft document
- Cancelled for being too much detail
- Cancelled because we don’t really need this
- In particular, icaz02040 looks interesting
- This is F373N, so [O II]
- Also, icaz01060, although mainly continuum, does have something else
- It is F336W, so not sure what it has
- Seems to be mainly scattered starlight plus the Balmer continuum
- It is the other side of the Balmer jump
- icaz01050 is F224W and that is nearly all scattered, and with ring artefacts around bright stars
- icaz02030 is F280N (what line is that?) and it is nearly all ringing artefacts. But the rings don’t line up with F224W so we can’t use it to remove them.
- We have the following epochs
- 1996 WFPC2/PC f656n
- 2005 ACS/WFC f658n
- 2015 WFC3/UVIS f656n
- There are lots of intermediate epochs between 1996 and 2005, but they are all with the WFC chip, and so lower resolution.
- We already have the correct WCS for the 1996 and 2015 images
- 1996 in hh204/
- fov2f656-PC.fits
- fov2f656-align-Robberto-PC.wcs
- 2015 in ~/Dropbox/Orion-HH-data/ACS/
- icaz01040_drz.fits
- trap-south-align-robberto.wcs
- Turns out that was not good enough - redo with star NIC 1211 as above in Align the images to Robberto 2005
- Result hh204/trap-south-hh204-align-robberto.wcs
- 2005 also in ~/Dropbox/Orion-HH-data/ACS/
- hlsp_orion_hst_acs_strip0l_f658n_v1_drz.fits
- This is a very large file - we should extract a sub-region to be more manageable
- I already did this in [[id:904807B0-9F3B-4299-840E-096A910331E8][Extract just a window of the ACS 2005 Ha+[N II] image]]
- 1996 in hh204/
- x, y, w, h = 7292.1998,2884.8059,2078.195,1347.475
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("../Orion-HH-data/ACS")
# center (x0, y0) and width, height of window
x0, y0, w, h = 7292, 2885, 2000, 2000
# Corners of window to extract from image
xmin, xmax, ymin, ymax = x0 - w//2, x0 + w//2, y0 - w//2, y0 + w//2
rect = slice(ymin, ymax), slice(xmin, xmax)
# Big file from Robberto 2005 program
hdu0 = fits.open(datadir / "hlsp_orion_hst_acs_strip0l_f658n_v1_drz.fits")["SCI"]
# Fix stupid SIP keywords that shouldn't be there
bad_kwds = [_ for _ in hdu0.header if _[:2] in ("A_", "B_")]
for kwd in bad_kwds:
hdu0.header.remove(kwd)
# Now we are ready to create the WCS
w = WCS(hdu0.header)
outdir = Path("hh204")
fits.PrimaryHDU(
data=hdu0.data[rect],
header=w.slice(rect).to_header(),
).writeto(
outdir / "hh204_acs_f658n_2005.fits",
overwrite=True,
)
- This had very annoying errors to begin with because I was trying to update the headers, rather than simply replacing them
- It turns out that is very difficult to get right because of the byzantine preferences of CD and PC, etc. So I just don’t bother now
- This means that we lose info on provenance, but at least the images align properly in DS9
from pathlib import Path import numpy as np from astropy.io import fits from astropy.wcs import WCS from reproject import reproject_interp import wcsfile datadir = Path("hh204") bigdatadir = Path("../Orion-HH-data/ACS") # Orig file from Bally 2015 program hdu2 = fits.open(bigdatadir / "icaz01040_drz.fits")["SCI"] # Remove the CD keywords because they cause problems bad_kwds = [_ for _ in hdu2.header if _.startswith("CD")] for kwd in bad_kwds: hdu2.header.remove(kwd) wcs_aligned = WCS(wcsfile.read(datadir / "trap-south-hh204-align-robberto.wcs")) hdu2.header.update(wcs_aligned.to_header()) # Cropped file from Robberto 2005 program hdu1, = fits.open(datadir / "hh204_acs_f658n_2005.fits") # Reproject 2->1 newdata, _ = reproject_interp(hdu2, hdu1.header) newheader = hdu1.header fits.PrimaryHDU( data=newdata, header=newheader, ).writeto( datadir / "hh204_wfc3_f656n_2015.fits", overwrite=True, ) # Orig file from Bally 1996 program hdu3 = fits.open(datadir / "fov2f656-PC.fits")[0] # Remove the CD keywords because they cause problems bad_kwds = [_ for _ in hdu3.header if _.startswith("CD")] for kwd in bad_kwds: hdu3.header.remove(kwd) wcs_aligned = WCS(wcsfile.read(datadir / "fov2f656-align-Robberto-PC.wcs")) hdu3.header.update(wcs_aligned.to_header()) # Reproject 3->1 newdata, _ = reproject_interp(hdu3, hdu1.header) newheader = hdu1.header fits.PrimaryHDU( data=newdata, header=newheader, ).writeto( datadir / "hh204_wfpc2_f656n_1996.fits", overwrite=True, )
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("hh204")
hdu1, = fits.open(datadir / "hh204_wfpc2_f656n_1996.fits")
hdu2, = fits.open(datadir / "hh204_acs_f658n_2005.fits")
hdu3, = fits.open(datadir / "hh204_wfc3_f656n_2015.fits")
fits.PrimaryHDU(
header=hdu3.header,
data=hdu3.data/hdu2.data,
).writeto(
datadir / f"hh204_f656n_ratio_2015_2005.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu2.header,
data=hdu2.data/hdu1.data,
).writeto(
datadir / f"hh204_f656n_ratio_2005_1996.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu3.header,
data=hdu3.data/hdu1.data,
).writeto(
datadir / f"hh204_f656n_ratio_2015_1996.fits",
overwrite=True,
)
import sys
sys.path.append("../../Source/pyflct")
import os
from pathlib import Path
import numpy as np
import pyflct
from astropy.io import fits
from astropy.wcs import WCS
os.environ["PATH"] = ":".join((os.environ["PATH"], "/Users/will/.local/bin"))
datadir = Path("hh204")
h1, = fits.open(datadir / "hh204_wfpc2_f656n_1996.fits")
h2, = fits.open(datadir / "hh204_acs_f658n_2005.fits")
data_scale = h1.data.max()
sigma = 10.0
vx, vy, vm = pyflct.flct(h1.data.T, h2.data.T,
deltat=1.0, deltas=1.0, sigma=sigma,
thresh=0.1/data_scale,
)
vx[vm==0] = np.nan
vy[vm==0] = np.nan
fits.PrimaryHDU(
header=h1.header,
data=vx.T,
).writeto(datadir / f"hh204_f656n_1996_2005_vx_sig{int(sigma):02d}.fits", overwrite=True)
fits.PrimaryHDU(
header=h1.header,
data=vy.T,
).writeto(datadir / f"hh204_f656n_1996_2005_vy_sig{int(sigma):02d}.fits", overwrite=True)
import sys
sys.path.append("../../Source/pyflct")
import os
from pathlib import Path
import numpy as np
import pyflct
from astropy.io import fits
from astropy.wcs import WCS
os.environ["PATH"] = ":".join((os.environ["PATH"], "/Users/will/.local/bin"))
datadir = Path("hh204")
h1, = fits.open(datadir / "hh204_acs_f658n_2005.fits")
h2, = fits.open(datadir / "hh204_wfc3_f656n_2015.fits")
data_scale = h1.data.max()
sigma = 10.0
vx, vy, vm = pyflct.flct(h1.data.T, h2.data.T,
deltat=1.0, deltas=1.0, sigma=sigma,
thresh=0.1/data_scale,
)
vx[vm==0] = np.nan
vy[vm==0] = np.nan
fits.PrimaryHDU(
header=h1.header,
data=vx.T,
).writeto(datadir / f"hh204_f656n_2005_2015_vx_sig{int(sigma):02d}.fits", overwrite=True)
fits.PrimaryHDU(
header=h1.header,
data=vy.T,
).writeto(datadir / f"hh204_f656n_2005_2015_vy_sig{int(sigma):02d}.fits", overwrite=True)
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("hh204")
hdu, = fits.open(datadir / "hh204_wfc3_f656n_2015.fits")
s1, s2 = "", "_sig10"
s_interval = "f656n_1996_2005"
vx1 = fits.open(datadir / f"hh204_{s_interval}{s1}_vx{s2}.fits")[0].data
vy1 = fits.open(datadir / f"hh204_{s_interval}{s1}_vy{s2}.fits")[0].data
m = vx1 == -20.0
vx1[m] = np.nan
vy1[m] = np.nan
s_interval = "f656n_2005_2015"
vx2 = fits.open(datadir / f"hh204_{s_interval}{s1}_vx{s2}.fits")[0].data
vy2 = fits.open(datadir / f"hh204_{s_interval}{s1}_vy{s2}.fits")[0].data
vxstack = np.stack([vx1, vx2])
vystack = np.stack([vy1, vy2])
vxm = np.nanmean(vxstack, axis=0)
vym = np.nanmean(vystack, axis=0)
vxs = np.nanstd(vxstack, axis=0)
vys = np.nanstd(vystack, axis=0)
vv = np.hypot(vxm, vym)
pa = np.rad2deg(np.arctan2(-vxm, vym)) % 360
sigma = np.hypot(vxs, vys) / vv
s_interval = "f656n_1996_2015_2stage"
fits.PrimaryHDU(
header=hdu.header,
data=vxm,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vx_mean{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vym,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vy_mean{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vxs,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vx_sigma{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vys,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vy_sigma{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=sigma,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vv_sigrel{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vv,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vv{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=pa,
).writeto(
datadir / f"hh204_{s_interval}{s1}_pa{s2}.fits",
overwrite=True,
)
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
figdir = Path(".")
datadir = Path("hh204")
s0 = "f656n_1996_2015_2stage"
figfile = figdir / f"hh204_{s0}_proper_motions.jpg"
hdu, = fits.open(datadir / "hh204_acs_f658n_2005.fits")
w = WCS(hdu)
im = hdu.data
#s1, s2 = "_sharp_16", "_sig20"
s1, s2 = "", "_sig10"
vx = fits.open(datadir / f"hh204_{s0}{s1}_vx_mean{s2}.fits")[0].data
vy = fits.open(datadir / f"hh204_{s0}{s1}_vy_mean{s2}.fits")[0].data
sig = fits.open(datadir / f"hh204_{s0}{s1}_vv_sigrel{s2}.fits")[0].data
vv = np.hypot(vx, vy)
pa = np.rad2deg(np.arctan2(-vx, vy)) % 360
vx[abs(vv) > 7.0] = np.nan
vx[vv < 0.1] = np.nan
vx[(im < 20.0)] = np.nan
vx[(im > 90.0)] = np.nan
vx[(sig > 0.5)] = np.nan
ny, nx = im.shape
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
fig, ax = plt.subplots(figsize=(4.5, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=10, vmax=90, origin="lower", cmap="gray_r")
step = slice(3, None, 7), slice(4, None, 7)
# Miss out the regions that are not in all three epochs
#
# Q = ax.quiver(X[step], Y[step], vx[step], vy[step], pivot="middle",
# headwidth=4, headlength=6, minlength=0.01,
# units="xy", scale=0.125, width=1.5, minshaft=2.0, color="orange", alpha=0.7)
# Q.set_path_effects(
# [
# path_effects.PathPatchEffect(
# edgecolor='white', linewidth=0.1, facecolor='orange', capstyle="projecting", alpha=0.7
# ),
# ]
# )
vx[(sig == 0.0)] = np.nan
Q = ax.quiver(X[step], Y[step], vx[step], vy[step], pivot="middle",
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.3, width=1.5, minshaft=2.0, color="r", alpha=1.0)
Q.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
QK = ax.quiverkey(Q, 0.1, 0.2, 7.0,
"70 km/s", labelcolor="red",
fontproperties=dict(size=6))
QK.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='red', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
ax.set(
xlabel="RA",
ylabel="Dec",
xlim=[370, 1050],
ylim=[470, 1150],
)
fig.tight_layout(rect=(0.2, 0.1, 1.0, 1.0))
fig.savefig(figfile, dpi=1200)
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.table import Table
import regions as rg
regiondir = Path("hh204")
datadir = Path("hh204")
regions = rg.read_ds9(regiondir / "HH204-samples.reg")
# Split into BG and source regions
source_regions = []
bg_regions = []
for region in regions:
label = region.meta["label"]
if "BG" in label:
bg_regions.append(region)
else:
source_regions.append(region)
hdu2015, = fits.open(datadir / "hh204_wfc3_f656n_2015.fits")
w = WCS(hdu2015.header)
badmask = ~np.isfinite(hdu2015.data)
hdu2015.data[badmask] = 0.0
data_images = {
"VV": fits.open(datadir / "hh204_f656n_1996_2015_2stage_vv_sig10.fits")[0].data,
"PA": fits.open(datadir / "hh204_f656n_1996_2015_2stage_pa_sig10.fits")[0].data,
}
# Find mean brightness BG for each group
weights = np.zeros_like(hdu2015.data)
for r in bg_regions:
weights += r.to_pixel(w).to_mask().to_image(hdu2015.data.shape)
bg_brightness = np.average(hdu2015.data, weights=weights)
# Now do measurements on each source region
all_results = []
for region in source_regions:
label = region.meta["label"]
group = label.split()[2]
bg = bg_brightness
results = {"region": label}
mask = region.to_pixel(w).to_mask()
mw = mask.to_image(hdu2015.data.shape).astype(bool)
weights = hdu2015.data[mw] - bg
results["Bright"] = np.round(np.mean(weights), 2)
results["Bright sig"] = np.round(np.std(weights), 2)
results["Bright BG"] = np.round(bg, 2)
for data_label, data in data_images.items():
m = mask.to_image(data.shape).astype(bool)
mean = np.average(data[m], weights=weights)
sigma = np.sqrt(np.average((data[m] - mean)**2, weights=weights))
results[data_label] = np.round(mean, 2)
results[f"{data_label} sig"] = np.round(sigma, 2)
all_results.append(results)
#
# List of dicts to astropy.table.Table
#
# It would be cleaner to just write Table(rows=all_results) but for
# some reason that does alphabetical sorting on column names, which we
# don't want. This way preserves the column order.
table = Table(
names=all_results[0].keys(),
rows=[_.values() for _ in all_results],
)
#
# astropy.table.Table to a list of lists, which babel will format as
# an org table
#
display = [list(table.colnames), None] + [list(row) for row in table]
region | Bright | Bright sig | Bright BG | VV | VV sig | PA | PA sig |
---|---|---|---|---|---|---|---|
HH 204 BS a | 6.09 | 1.3 | 4.01 | 4.97 | 0.21 | 127.08 | 4.27 |
HH 204 MD a | 8.76 | 1.27 | 4.01 | 4.87 | 0.47 | 129.9 | 8.97 |
HH 204 MD b | 9.47 | 0.96 | 4.01 | 4.76 | 0.7 | 131.87 | 6.89 |
HH 204 MD c | 9.07 | 0.97 | 4.01 | 5.18 | 1.1 | 144.15 | 39.09 |
HH 204 MD d | 7.64 | 1.2 | 4.01 | 5.37 | 0.88 | 140.11 | 30.44 |
HH 204 MD e | 8.24 | 1.72 | 4.01 | 4.86 | 1.02 | 152.97 | 52.87 |
HH 204 BS b | 4.37 | 0.8 | 4.01 | 5.18 | 0.57 | 134.61 | 5.02 |
HH 204 BS c | 2.96 | 0.63 | 4.01 | 3.37 | 0.17 | 96.48 | 11.26 |
HH 204 BS d | 3.12 | 0.3 | 4.01 | 2.22 | 0.18 | 97.12 | 6.84 |
HH 204 Knot | 2.21 | 0.63 | 4.01 | 3.35 | 0.16 | 183.61 | 2.39 |
- For an average motion, I take the
BS a
,MD a
,MD b
values to get VV = 4.9, PA = 130- 1 pixel in 9.91 yrs is 10.0 km/s within the uncertainties
- So, we can move the CRPIX values for the 1996 and 2015 epochs
- x : -4.9 sin(130) = -3.75
- y : 4.9 cos(130) = -3.15
- This turns out to be not quite right for the nose
- What we end up with is (x, y)
- 1996
- (-5.5, -6) = 8
- 2015
- (-4.5, -4.3) = 6.2
- average
- (-5.0, -5.2)
- These are saved in
- What we end up with is (x, y)
- And the implications for the velocity are:
- 1996 - 2005: 80 km/s at PA=137
- 2005 - 2015: 62 km/s at PA=134
- Average: 71 +/- 9 km/s at 136 +/- 2 deg
- So now we can find the inclination to the plane of the sky
- The radial velocity of the HH 204 knot is -45 km/s
- So total velocity is V = hypot(71, 45) = 84 +/- 9 km/s
- Inclination is -32 +/- 6 deg
- Previous determination was 27 deg, so this is consistent
- Note that for HH 203, we get a different PM
- roughly 12 km/s at PA=330 in HH 204 frame:
- x = -1.2 sin(330) = 0.6
- y = 1.2 cos(330) = 1.0
- So in cluster frame: (-4.4, -4.2)
- 6.1 = 61 km/s
- PA = 134 deg
- And the radial velocity is higher: -73 km/s
- HH 203 total velocity: 95 km/s
- HH 203 inclination is -50 +/- 6 deg
- roughly 12 km/s at PA=330 in HH 204 frame:
- Use results from previous calculation and save resampled images that are aligned in the nose frame
from pathlib import Path import numpy as np from astropy.io import fits from astropy.wcs import WCS from reproject import reproject_interp import wcsfile datadir = Path("hh204") # Orig file from Bally 2015 program hdu2, = fits.open(datadir / "hh204_wfc3_f656n_2015.fits") # Align on bowshock nose wcs_aligned = WCS(wcsfile.read(datadir / "hh204-2015-nose-frame.wcs")) hdu2.header.update(wcs_aligned.to_header()) # Cropped file from Robberto 2005 program hdu1, = fits.open(datadir / "hh204_acs_f658n_2005.fits") # Reproject 2->1 newdata, _ = reproject_interp(hdu2, hdu1.header) newheader = hdu1.header fits.PrimaryHDU( data=newdata, header=newheader, ).writeto( datadir / "hh204_wfc3_f656n_2015_nose_align.fits", overwrite=True, ) # Orig file from Bally 1996 program hdu3, = fits.open(datadir / "hh204_wfpc2_f656n_1996.fits") wcs_aligned = WCS(wcsfile.read(datadir / "hh204-1996-nose-frame.wcs")) hdu3.header.update(wcs_aligned.to_header()) # Reproject 3->1 newdata, _ = reproject_interp(hdu3, hdu1.header) newheader = hdu1.header fits.PrimaryHDU( data=newdata, header=newheader, ).writeto( datadir / "hh204_wfpc2_f656n_1996_nose_align.fits", overwrite=True, )
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("hh204")
hdu1, = fits.open(datadir / "hh204_wfpc2_f656n_1996_nose_align.fits")
hdu2, = fits.open(datadir / "hh204_acs_f658n_2005.fits")
hdu3, = fits.open(datadir / "hh204_wfc3_f656n_2015_nose_align.fits")
fits.PrimaryHDU(
header=hdu3.header,
data=hdu3.data/hdu2.data,
).writeto(
datadir / f"hh204_f656n_ratio_2015_2005_nose_align.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu2.header,
data=hdu2.data/hdu1.data,
).writeto(
datadir / f"hh204_f656n_ratio_2005_1996_nose_align.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu3.header,
data=hdu3.data/hdu1.data,
).writeto(
datadir / f"hh204_f656n_ratio_2015_1996_nose_align.fits",
overwrite=True,
)
for f in wfc3_f656n_2015 acs_f658n_2005 wfpc2_f656n_1996; do
python /Users/will/Dropbox/Teresa-Turtle/scripts/smooth-image.py hh204_$f.fits -w 16
done
for f in wfc3_f656n_2015 wfpc2_f656n_1996; do
python /Users/will/Dropbox/Teresa-Turtle/scripts/smooth-image.py hh204_${f}_nose_align.fits -w 16
done
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
figdir = Path(".")
datadir = Path("hh204")
for s in "wfpc2_f656n_1996_nose_align", "acs_f658n_2005", "wfc3_f656n_2015_nose_align":
camera, filter_, year = s.split("_")[:3]
figfile = figdir / f"hh204_{year}_nose_align.jpg"
hdu, = fits.open(datadir / f"hh204_{s}_sharp_16.fits")
w = WCS(hdu)
im = hdu.data
fig, ax = plt.subplots(figsize=(4, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=0.5, vmax=2.0, origin="lower", cmap="gray")
ax.set(
xlim=[370, 1050],
ylim=[470, 1150],
)
lon, lat = ax.coords
lon.set_ticklabel_visible(False)
lat.set_ticklabel_visible(False)
fig.tight_layout()
fig.savefig(figfile, dpi=600)
open hh204_*_nose_align.jpg
New [2021-06-21 Mon] - do the same for the star-aligned images. And not spatially filtered
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
figdir = Path("press-release-hh204")
datadir = Path("hh204")
vlim = {
"1996": [40, 400],
"2005": [6, 100],
"2015": [2, 20],
}
for s in "wfpc2_f656n_1996", "acs_f658n_2005", "wfc3_f656n_2015":
camera, filter_, year = s.split("_")[:3]
figfile = figdir / f"hh204_{year}.jpg"
hdu, = fits.open(datadir / f"hh204_{s}.fits")
w = WCS(hdu)
im = hdu.data
im_med = np.median(hdu.data)
fig, ax = plt.subplots(figsize=(4, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=vlim[year][0], vmax=vlim[year][1], origin="lower", cmap="gray")
ax.set(
xlim=[370, 1050],
ylim=[470, 1150],
)
lon, lat = ax.coords
lon.set_ticklabel_visible(False)
lat.set_ticklabel_visible(False)
fig.tight_layout()
fig.savefig(figfile, dpi=600)
ORDER="1996 2005 2015"
convert $(for a in "$ORDER"; do printf -- "-delay 50 %s.jpg " $a; done; ) HH204-bowshock-frame.gif
open -a Safari hh204_multi_epoch_nose_align_EXPORT/HH204-bowshock-frame.gif
ORDER="1996 2005 2015"
convert $(for a in "$ORDER"; do printf -- "-delay 50 hh204_%s.jpg " $a; done; ) HH204-cluster-frame.gif
open -a Safari press-release-hh204/HH204-cluster-frame.gif
Write a general program to do the FLCT method
import sys
sys.path.append("../../Source/pyflct")
import os
from pathlib import Path
import numpy as np
import pyflct
from astropy.io import fits
from astropy.wcs import WCS
os.environ["PATH"] = ":".join((os.environ["PATH"], "/Users/will/.local/bin"))
try:
DATADIR, FN1, FN2, OUT_PREFIX = sys.argv[1:5]
except:
sys.exit(f"Usage: {sys.argv[0]} DATADIR FN1 FN2 OUT_PREFIX [SIGMA]")
try:
SIGMA = float(sys.argv[5])
except:
# Optional input argument
SIGMA = 10.0
data_path = Path(DATADIR)
h1, = fits.open(data_path / f"{FN1}.fits")
h2, = fits.open(data_path / f"{FN2}.fits")
# Remove NaNs in the input data
data1 = h1.data.T
data2 = h2.data.T
goodpixels = np.isfinite(data1*data2)
data1[~goodpixels] = 0.0
data2[~goodpixels] = 0.0
print("Number of goodpixels =", goodpixels.sum())
# Designed for images that have been high-pass filtered and normalized
data_scale = 2.0
sigma = SIGMA
vx, vy, vm = pyflct.flct(data1, data2,
deltat=1.0, deltas=1.0, sigma=sigma,
#thresh=0.1/data_scale,
)
#vx[vm==0] = np.nan
#vy[vm==0] = np.nan
fits.PrimaryHDU(
header=h1.header,
data=vx.T,
).writeto(data_path / f"{OUT_PREFIX}_vx_sig{int(sigma):02d}.fits", overwrite=True)
fits.PrimaryHDU(
header=h1.header,
data=vy.T,
).writeto(data_path / f"{OUT_PREFIX}_vy_sig{int(sigma):02d}.fits", overwrite=True)
Recommended to run these in a shell buffer since they take a few minutes each
python make_pm_map.py hh204\
hh204_wfpc2_f656n_1996_nose_align_sharp_16\
hh204_acs_f658n_2005_sharp_16\
hh204_nose_align_f656n_1996_2005\
> hh204_nose_align_f656n_1996_2005_make_pm_map.out
python make_pm_map.py hh204 hh204_acs_f658n_2005_sharp_16 hh204_wfc3_f656n_2015_nose_align_sharp_16 hh204_nose_align_f656n_2005_2015 > hh204_nose_align_f656n_2005_2015_make_pm_map.out
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("hh204")
hdu, = fits.open(datadir / "hh204_wfc3_f656n_2015_nose_align_sharp_16.fits")
s1, s2 = "", f"_sig{SIG}"
s_interval = "nose_align_f656n_1996_2005"
vx1 = fits.open(datadir / f"hh204_{s_interval}{s1}_vx{s2}.fits")[0].data
vy1 = fits.open(datadir / f"hh204_{s_interval}{s1}_vy{s2}.fits")[0].data
m = vx1 == -20.0
vx1[m] = np.nan
vy1[m] = np.nan
s_interval = "nose_align_f656n_2005_2015"
vx2 = fits.open(datadir / f"hh204_{s_interval}{s1}_vx{s2}.fits")[0].data
vy2 = fits.open(datadir / f"hh204_{s_interval}{s1}_vy{s2}.fits")[0].data
vxstack = np.stack([vx1, vx2])
vystack = np.stack([vy1, vy2])
vxm = np.nanmean(vxstack, axis=0)
vym = np.nanmean(vystack, axis=0)
vxs = np.nanstd(vxstack, axis=0)
vys = np.nanstd(vystack, axis=0)
vv = np.hypot(vxm, vym)
pa = np.rad2deg(np.arctan2(-vxm, vym)) % 360
sigma = np.hypot(vxs, vys) / vv
s_interval = "nose_align_f656n_1996_2015_2stage"
fits.PrimaryHDU(
header=hdu.header,
data=vxm,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vx_mean{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vym,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vy_mean{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vxs,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vx_sigma{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vys,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vy_sigma{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=sigma,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vv_sigrel{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vv,
).writeto(
datadir / f"hh204_{s_interval}{s1}_vv{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=pa,
).writeto(
datadir / f"hh204_{s_interval}{s1}_pa{s2}.fits",
overwrite=True,
)
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.patches import Rectangle
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
import seaborn as sns
sns.set_color_codes()
figdir = Path(".")
datadir = Path("hh204")
s0 = "nose_align_f656n_1996_2015_2stage"
figfile = figdir / f"hh204_{s0}_proper_motions.jpg"
hdu, = fits.open(datadir / "hh204_acs_f658n_2005.fits")
im_full = hdu.data
hdu, = fits.open(datadir / "hh204_acs_f658n_2005_sharp_16.fits")
w = WCS(hdu)
im = hdu.data
#s1, s2 = "_sharp_16", "_sig20"
s1, s2 = "", "_sig16"
vx = fits.open(datadir / f"hh204_{s0}{s1}_vx_mean{s2}.fits")[0].data
vy = fits.open(datadir / f"hh204_{s0}{s1}_vy_mean{s2}.fits")[0].data
sig = fits.open(datadir / f"hh204_{s0}{s1}_vv_sigrel{s2}.fits")[0].data
vv = np.hypot(vx, vy)
pa = np.rad2deg(np.arctan2(-vx, vy)) % 360
bad_pixels = (
(abs(vv) > 8.0) | (vv < 0.01)
| (im < 0.95) | (im > 4.0)
| (sig > 1.0) | (sig == 0.0)
)
bad_pixels = bad_pixels | ~np.isfinite(im)
ny, nx = im.shape
x, y = np.meshgrid(np.arange(nx), np.arange(ny))
# Flat copies of the arrays with only the desired pixels
X = x[~bad_pixels]
Y = y[~bad_pixels]
VX = vx[~bad_pixels]
VY = vy[~bad_pixels]
VV = vv[~bad_pixels]
IM = im[~bad_pixels]
IM_full = im_full[~bad_pixels]
# divide into two parts:
# vv < 2 that we plot densely with arrows
# vv > 2 that we plot more sparsely
# vv > 4 that we plot even more sparsely
VX_slow = VX.copy()
VX_fast = VX.copy()
VX_slow[VV > 2.0] = np.nan
VX_slow[IM_full < 20.0] = np.nan
VX[VV < 2.0] = np.nan
VX[VV > 4.0] = np.nan
VX[IM_full < 20.0] = np.nan
VX_fast[VV < 4.0] = np.nan
VX_fast[IM < 1.02] = np.nan
fig, ax = plt.subplots(figsize=(4.5, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=0.5, vmax=1.5, origin="lower", cmap="gray")
n = len(X)
#fraction = 0.2
fraction = 0.1
rng = np.random.default_rng()
nsamples = int(fraction * n)
weight = IM # **2
samples = rng.choice(np.arange(n), nsamples, replace=False, p=weight/weight.sum())
# Slow parts
Qslow = ax.quiver(X[samples], Y[samples], VX_slow[samples], VY[samples], pivot="middle",
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.2, width=1.5, minshaft=2.0, color="r", alpha=1.0)
Qslow.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='w', linewidth=0.1, facecolor='r', capstyle="projecting",
),
]
)
# Medium parts
step = 4
Q = ax.quiver(X[samples[::step]], Y[samples[::step]],
VX[samples[::step]], VY[samples[::step]],
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.2, width=1.5, minshaft=2.0,
color="r",
alpha=1.0)
Q.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='w', linewidth=0.1, facecolor='r', capstyle="projecting",
),
]
)
# Fast parts
step = 8
Qfast = ax.quiver(X[samples[::step]], Y[samples[::step]],
VX_fast[samples[::step]], VY[samples[::step]],
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.2, width=1.5, minshaft=2.0, color="r", alpha=0.5)
Qfast.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='w', linewidth=0.1, facecolor='r', capstyle="projecting", alpha=0.5,
),
]
)
ax.add_patch(
Rectangle(
(0.015, 0.48),
width=0.2, height=0.18,
facecolor="white", alpha=0.5,
transform=ax.transAxes)
)
y0 = 0.52
for u0 in [2.0, 4.0, 6.0]:
QK = ax.quiverkey(Qslow, 0.07, y0, u0,
f"{int(10*u0):d} km/s", labelcolor="r", labelpos="E",
fontproperties=dict(size=6))
y0 += 0.05
ax.set(
xlabel="RA",
ylabel="Dec",
xlim=[370, 1050],
ylim=[470, 1150],
)
fig.tight_layout(rect=(0.2, 0.1, 1.0, 1.0))
fig.savefig(figfile, dpi=1200)
- We see a solid lump of emission that is all moving at the bow shock velocity
- This is almost stationary in the movie
- It seems to be the inner part of the working surface
- Part of it is the clumpy kite-shaped mass at the head
- And there are also faint tendrils feeding into this from behind
- Some of these may have low velocities towards the working surface, which is evidence for feeding from a dense jet
- Then there is a sheath of backwards moving material in the bow shock frame
- In principle we can use this to try and estimate the shock velocity (and hence whether the upstream gas is moving)
- We could try and find the expansion law for the bow shock wings, or at least fit a theoretical model of how the motion in the bow shock frame should be (conserving velocity parallel to the shock)
- We can try and colorize the image so that the stationary gas is blue amd the moving gas is red, for instance
- Possible sources of the PA140 flow
- Proper motions in the same direction in the Huygens region from Table 6 of ODell:2015a
- 178-328 has V_t = 161 km/s @ PA132 in [O III] only
- HH204 is at PA146 from this source
- 178-328 has V_t = 161 km/s @ PA132 in [O III] only
- Proper motions in the same direction in the Huygens region from Table 6 of ODell:2015a
- The PA108 flow looks promising
- There is evidence for a long filament in this direction from the following lines
- [Fe III] 4658 integrated line
- [Fe III] 5270 integrated line (slightly weaker)
- Si II 6347 blue flank
- Si II 6371 blue flank
- Unfortunately, the Si II lines at 5958, 5979 do not show it
- He I 6678.15 blue flank (specifically, 6676.65 channel)
- Velocity = 3e5 (6676.65 - 6678.15) / 6678.15 = -70 in some frame or other
- Also evidence from the KPNo vcube in Ha, [O III], and [N II]
- But it is difficult to see because of ovelap with HH 528
- There are several small-scale filaments that cross the HH 528 flow with approximately the right orientaion
- They are in the general vicinity of 166-442
- Need to check their proper motions
- Also closer to bow aroud 205-507
- Evidence from H_2 filaments
- Subaru map shows a filament that aligns with this flow, just to NW of the HH 528 head
- And we can also see filaments that maybe feed HH 203
- And possibly the Big Arc flow too
- In fact, the strange partial bubble that I once associated with HH 528 now turns out to be a combination of several unconnected filaments
- /Users/will/Dropbox/dib-scatter-hii/data/orion-hawki/Orion-Bar-H2-minus-045-Brg.fits
- https://www.irya.unam.mx/gente/w.henney/WebWiki/OrionImages.html
- Possible source:
- V2235 Ori, JW401, COUP 534
- This is the most luminous star vaguely along the PA 108 line
- But it doesn’t work for passing through the H_2 filament
- 05 35 13.2006958872 -05 24 55.493104899
- COUP 484 is more likely
- This does line up and gives PA=105, which is close enough
- 05 35 12.4954 -05 24 38.148
- 125-438
- V2235 Ori, JW401, COUP 534
- Although none of these PA seem compatible with the observed proper motion of HH 204
- Or are they?
- There is evidence for a long filament in this direction from the following lines
Fe, Ni, Cr abundance higher by factor of about 10 in shocked gas
- Gets high ADF for C2+ but
- might be fluorescence if it is from the 723X lines
- No, seems mainly based on 4267 which should be unaffected
- Also it is based on 3rd party [C III] NUV observations
- might be fluorescence if it is from the 723X lines
- Born and brought up here in Morelia, Michoacán
- Went to UNAM Mexico City to do an undergraduate degree in physics 2012-2017
- Undergraduate thesis with Manuel Peimbert in IA-UNAM on Temperature fluctuations in Planetary Nebulae
- In 2014 he received Premio Estatal al mérito juvenil del estado de Michoacán
- Then did Master’s degree in Universidad de la Laguna and Instituto de Astrofísica de Canarias in Tenerife
- Research on Abundance Gradient of Helium in Milky Way under direction of César Esteban
- Finally, he is currently do a PhD at the IAC where he is also a resident astronomer
- Advisors are César Esteban and Jorge García Rojas
- Thesis is on Photoionized Herbig-Haro Objects in the Orion Nebula
- which is topic he will talk about today
Nací y viví en Morelia hasta mudarme a la CDMX para estudiar la licenciatura en Física en la UNAM del 2012 al 2017. Ahí hice mi tesis de licenciatura sobre fluctuaciones de temperatura en nebulosas planetarias bajo la dirección del Dr. Manuel Peimbert, quien me introdujo en la Astronomía. Por mi desempeño académico me dieron el Premio Estatal al mérito juvenil del estado de Michoacán 2014. Al terminar la licenciatura, obtuve la beca de la Fundación Carolina para Latinoamérica para estudiar el Máster en Astrofísica en la Universidad de La Laguna/IAC. Hice mi trabajo de fin de Máster bajo la dirección de César Esteban sobre el gradiente de Helio en la Vía Láctea. Luego gané una beca de Astrofísico Residente del IAC donde hago mi doctorado. Mi proyecto de doctorado es sobre Objetos Herbig-Haro fotoionizados en la Nebulosa de Orión, bajo la dirección de César Esteban y Jorge García Rojas.
Dear Dr. Mendez-Delgado, I don’t know that we have met either in my one visit to IAC years ago or at any of the Orion Nebula conferences held since then. In any event, I am Bob O’Dell (C. R. O’Dell on papers), presently a part-time Distinguished Research Professor at Vanderbilt University in Nashville, Tennessee USA. You’ve generously cited my relevant papers in your own publications. I’ve read with great interest your papers “Photoionized Herbig-Haro objects in the Orion Nebula through deep high-spectral resolution spectroscopy I: HH 529 II and III (published earlier this year in MNRAS) and …. II: HH 204 (submitted to the ApJ but available through arXiv). A tremendous amount of effort must have gone into these two very thorough studies, congratulations. As someone who has worked on the Orion Nebula for 50+ years, I hope that you will welcome a few comments from me. At age 84 I can see that I won’t be able to continue my work on Orion much longer and these comments may encourage you to address some of the open questions. PAPER I: First a mea culpa. You point out at the end of Section 9.2 that your results for the tangential motions are in agreement with O’Dell and Henney 2008, but in disagreement with O’Dell et al. (2015). I have looked through my electronic notes used in the 2015 study, and cannot justify what I wrote there. Examining the original data, it is obvious that I made an error somewhere. Unfortunately, I don’t have access to my paper notes of that period (because of the pandemic) and cannot say anything further. Since you have already published your own results and pointed out the agreement with O & H 2008 and disagreement with O et al. 2015, I don’t think that I need to publish an erratum. Earlier in Section 9.2 you discuss the line broadening. You include the sigma-nt term. You probably know that I first noted that such a term was necessary in explaining Orion’s FWHM, concluding that this component carries a significant fraction of the total energy of the ionized gas. Although I’ve speculated on its cause, I think that your observation of the background gas emission may allow characterizing it in terms of ionization state and atomic mass in a superior fashion to what my collaborators and I have done before, perhaps leading to its explanation. I hope that you will look into this. If you or one of your collaborators cannot, please let me know and I’ll pursue it, using your published data. As you probably know, for several decades my collaborators and I have noted that there is usually a red wing on emission-lines with the best S/N ratio. I see this in your Fig.4, where it is expected (at twice the velocity difference of the emitting layer and the backscattering dust in the underlying PDR). Using V_PDR=28 and V_[OIII]=16, this leads to the expectation of the backscattered component to be at V=40, about where the red bump appears in the figure. Any backscattering of the shocked material will be much more redshifted and of course, weaker. Furthermore, the shocked gas is weak and probably not close to the PDR, so it is not surprising that its backscattered component is not seen. Of course we don’t know the wavelength dependence of the backscattering component or its dependence on the separation of the emitting and backscattering layers, so I cannot comment on the effect of this red component on derived characteristics such as abundances and physical conditions. PAPER-II: Fig.20. In this case the HH 204 component is strongest and the backscattered component may be present. It would be informative to deconvolve this spectrum into more than the three apparent components. I can see that there is the red wing on the Main Ionization Front component (labeled Orion Nebula). It is worth looking for the backscattered component of the “Blue Layer” (unlikely because of its weakness) and the HH 204 component. The latter should be highly redshifted and well defined. If it is not there, then this confirms that the HH 204 shock feature is not close to the underlying PDR. If you can’t do this, send me a velocity calibrated copy of the .fits image and I’ll do it for you. Section 8. I don’t see the PA 140 component on any of my images (the same as you use). Can you describe it more completely? End of paragraph 1. The values of 124 and 127 that I give are what I judged to be the axes of symmetry of HH 203 and HH 204. I see your component HH 203c, which appears to be the apex of a bowshock at 5:35:22.36-5:25:07.7 (I advise you give the coordinates). HH 203a (5:35:22.19:-5:25:05.5) and HH 203b (5:35:21.73-5:24:56.6). HH 203c is peculiar in that there are a pair of small shocks at 5:35:22.23-5:25:07.4 aimed about about PA 80 that overlap it. Personally, I despair of identifying all the driving jets and components of HH 203-204. The better one looks, the more one sees. For example the broad structure at 5:35:21.76-5:25:17.4 - does it fit in? `Section 9.1. I think that you have underestimated the amount of work done on the DBL, perhaps being led by co-author Will Henney to emphasize this feature in the vicinity of Theta2OriA. Will really deserves the credit for identifying the feature in the SE portion of the nebula as he attempted to interpret it (Deharveng’s study was a one-off presentation of some observations). My colleagues and I have written a series of papers which establish that this feature extends across much of the Huygens Region (O’Dell et al. 2020, ApJ, 891, 46) and have attempted to model its position (Abel et al 2019 ApJ, 881, 130). The most recent of these papers overlapped with the finalizing of your Paper-II and I’m not surprised that you do not mention them. The DBL is brighter and easier to study in the core of the Huygens Region, and in some areas becomes the strongest velocity component. I can’t say that we really understand it and we have another positional study underway. By this point we have an idea of the approximate position, certainly its position with respect to the multiple other foreground layers of gas. This leads us to strongly prefer the name “Nearer Ionized Layer” rather than the less descriptive “Blue Diffuse Layer”. Fig. 22, I see the Backscattered components of the Orion Nebula’s emission in the panels or [SIII] and [OIII]. I hope that you find these comments to be useful. At the minimum it means that there is at least one astronomer out in the world who has studied your excellent papers in detail. Best Regards, Bob O'Dell
- Disagreement between OH2008 and OD2015
Vt 08 Vro 08 HH 529-I 54 -31 HH 529-II 92 -31 HH 529-III 40 -30 - So recalculate the proper motions
- Used in OD2015
- GO 5085 - WFPC2 - O’Dell - 1995-01
- GO 5469 - WFPC2 - Bally - 1995-03
- Turns out I already aligned these to Robberto using Astrodrizzle
- See ~/Dropbox/OrionHST-2012/Combine/
- I have done all the PC frames in F502N but only HST3 frame (PC+WFC) in F656N
- But this is enough for HH 529
- GO 11038 - WFPC2 - 2007-11 - Biretta
- GO 12543 - WFC3 - 2012-01
- Newer and unused programs
- GO 13826 - ACS - 2015-10 - Robberto - F775W only
- GO 13419 - WFC3 - 2015-01 - Bally
- I already have these in ~/Dropbox/OrionBally-2016/data-2016-01-11/
- See also Proper motion of helical jet
- GO 9825 - ACS - 2004-01 - Bally - Ha
- These are longer exposures than Robberto
- Which is good for faint stuff
- But bad for saturated stars
- Also, in the HH 529 field they have the spikes at an annoying angle
- And the alignment isn’t quite right either
- These are longer exposures than Robberto
- GO 10246 - ACS - 2005-04 - Robberto - Ha and continuum only
- OH2008
- GO 10921 - WFPC2 - 2007-12 - O’Dell
- 5193, 5469 as above
- Also 8121, 9460 for HH 502 (way to S) and BN/KL
- Also
- Megeath Cycle 25 Program 15141
- Uses a novel drift and shift technique to survey much of the Orion molecular filament at 1.6 micron
- See Momcheva:2017a
- Epoch 2019-05 to 2019-08
- Might be useful, but data reduction could be difficult
- We have several options that use the same instrument:
- GO 12543 → 13419 is 3 years of WFC3 0.045 arcsec
- GO 10246 → 13826 is 10 years of ACS 0.05 arcsec, but F775W only
- Try this first
- 5085 → 11038 is 12 years of WFPC2 0.1 arcsec
- This is better because no stars
- All the files are in ~/Dropbox/Orion-HH-data/ACS/
- On the other hand, the filters aren’t quite the same
- But we do have the WFPC2 images from the Robberto program too
- Although unfortunately they are not aligned very well
- We also have both Bally and Robberto for 2005 ACS
- But Bally needs to be shifted a little bit
- In the F775W version below, I am now taking a very large window
- But here I want to concentrate on HH 529 for now
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("../Orion-HH-data/ACS")
# Corners of window to extract from image
xmin, xmax, ymin, ymax = 8120, 9670, 3770, 5710
rect = slice(ymin, ymax), slice(xmin, xmax)
# Big file from Robberto 2005 program
hdu0 = fits.open(datadir / "hlsp_orion_hst_acs_strip0l_f658n_v1_drz.fits")["SCI"]
# Fix stupid SIP keywords that shouldn't be there
bad_kwds = [_ for _ in hdu0.header if _[:2] in ("A_", "B_")]
for kwd in bad_kwds:
hdu0.header.remove(kwd)
# Now we are ready to create the WCS
w = WCS(hdu0.header)
outdir = Path("proper-motions/data")
fits.PrimaryHDU(
data=hdu0.data[rect],
header=w.slice(rect).to_header(),
).writeto(
outdir / "hh529_acs_f658n_2005.fits",
overwrite=True,
)
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from reproject import reproject_interp
datadir = Path("proper-motions/data")
bigdatadir = Path("../Orion-HH-data/ACS")
# Orig file from Bally 2015 program
hdu2 = fits.open(bigdatadir / "icaz01040_drz.fits")["SCI"]
# Use the alignment from trap-south-align-robberto.wcs
hdu2.header["CRPIX1"] = 2060
hdu2.header["CRPIX2"] = 2185.5
hdu2.header["CRVAL1"] = 83.8261121432188
hdu2.header["CRVAL2"] = -5.40646485004188
# Cropped file from Robberto 2005 program
hdu1, = fits.open(datadir / "hh529_acs_f658n_2005.fits")
# Reproject 2->1
newdata, _ = reproject_interp(hdu2, hdu1.header)
newheader = hdu2.header.copy()
newheader.update(WCS(hdu1.header).to_header())
fits.PrimaryHDU(
data=newdata,
header=newheader,
).writeto(
datadir / "hh529_wfc3_f656n_2015.fits",
overwrite=True,
)
# Orig file from Bally 1996 program
hdu3 = fits.open(bigdatadir / "F656N-HST3-align-rob_drz_sci.fits")[0]
# Alignment tweak
hdu3.header["CRVAL1"] += 0.05/3600
hdu3.header["CRVAL2"] -= 0.05/3600
# Reproject 2->1
newdata, _ = reproject_interp(hdu3, hdu1.header)
newheader = hdu3.header.copy()
newheader.update(WCS(hdu1.header).to_header())
fits.PrimaryHDU(
data=newdata,
header=newheader,
).writeto(
datadir / "hh529_wfpc2_f656n_1996.fits",
overwrite=True,
)
import sys
sys.path.append("../../Source/pyflct")
import os
from pathlib import Path
import numpy as np
import pyflct
from astropy.io import fits
from astropy.wcs import WCS
os.environ["PATH"] = ":".join((os.environ["PATH"], "/Users/will/.local/bin"))
datadir = Path("proper-motions/data")
h1, = fits.open(datadir / "hh529_wfpc2_f656n_1996.fits")
h2, = fits.open(datadir / "hh529_acs_f658n_2005.fits")
data_scale = h1.data.max()
sigma = 10.0
vx, vy, vm = pyflct.flct(h1.data.T, h2.data.T,
deltat=1.0, deltas=1.0, sigma=sigma,
thresh=0.1/data_scale,
)
vx[vm==0] = np.nan
vy[vm==0] = np.nan
fits.PrimaryHDU(
header=h1.header,
data=vx.T,
).writeto(datadir / f"hh529_f656n_1996_2005_vx_sig{int(sigma):02d}.fits", overwrite=True)
fits.PrimaryHDU(
header=h1.header,
data=vy.T,
).writeto(datadir / f"hh529_f656n_1996_2005_vy_sig{int(sigma):02d}.fits", overwrite=True)
import sys
sys.path.append("../../Source/pyflct")
import os
from pathlib import Path
import numpy as np
import pyflct
from astropy.io import fits
from astropy.wcs import WCS
os.environ["PATH"] = ":".join((os.environ["PATH"], "/Users/will/.local/bin"))
datadir = Path("proper-motions/data")
h1, = fits.open(datadir / "hh529_acs_f658n_2005.fits")
h2, = fits.open(datadir / "hh529_wfc3_f656n_2015.fits")
data_scale = h2.data.max()
sigma = 10.0
vx, vy, vm = pyflct.flct(h1.data.T, h2.data.T,
deltat=1.0, deltas=1.0, sigma=sigma,
thresh=0.1/data_scale,
)
vx[vm==0] = np.nan
vy[vm==0] = np.nan
fits.PrimaryHDU(
header=h1.header,
data=vx.T,
).writeto(datadir / f"hh529_f656n_2005_2015_vx_sig{int(sigma):02d}.fits", overwrite=True)
fits.PrimaryHDU(
header=h1.header,
data=vy.T,
).writeto(datadir / f"hh529_f656n_2005_2015_vy_sig{int(sigma):02d}.fits", overwrite=True)
import sys
sys.path.append("../../Source/pyflct")
import os
from pathlib import Path
import numpy as np
import pyflct
from astropy.io import fits
from astropy.wcs import WCS
os.environ["PATH"] = ":".join((os.environ["PATH"], "/Users/will/.local/bin"))
datadir = Path("proper-motions/data")
h1, = fits.open(datadir / "hh529_wfpc2_f656n_1996.fits")
h2, = fits.open(datadir / "hh529_wfc3_f656n_2015.fits")
data_scale = h2.data.max()
sigma = 5.0
vx, vy, vm = pyflct.flct(h1.data.T, h2.data.T,
deltat=1.0, deltas=1.0, sigma=sigma,
thresh=0.1/data_scale,
)
vx[vm==0] = np.nan
vy[vm==0] = np.nan
fits.PrimaryHDU(
header=h1.header,
data=vx.T,
).writeto(datadir / f"hh529_f656n_1996_2015_vx_sig{int(sigma):02d}.fits", overwrite=True)
fits.PrimaryHDU(
header=h1.header,
data=vy.T,
).writeto(datadir / f"hh529_f656n_1996_2015_vy_sig{int(sigma):02d}.fits", overwrite=True)
python /Users/will/Dropbox/Teresa-Turtle/scripts/smooth-image.py hh529_wfc3_f656n_2015.fits -w 16
python /Users/will/Dropbox/Teresa-Turtle/scripts/smooth-image.py hh529_acs_f658n_2005.fits -w 16
- First using the entire 20 year interval
- The problem with this is that things move to far so we need to use a wide kernel (15 pix) and even then it doesn’t pick up all the motions properly
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
figdir = Path("proper-motions")
datadir = figdir / "data"
s0 = "f656n_1996_2015"
figfile = figdir / f"hh529_{s0}_proper_motions.jpg"
hdu, = fits.open(datadir / "hh529_wfc3_f656n_2015_sharp_16.fits")
w = WCS(hdu)
im = hdu.data
#s1, s2 = "_sharp_16", "_sig20"
s1, s2 = "", "_sig15"
vx = fits.open(datadir / f"hh529_{s0}{s1}_vx{s2}.fits")[0].data
vy = fits.open(datadir / f"hh529_{s0}{s1}_vy{s2}.fits")[0].data
vv = np.hypot(vx, vy)
pa = np.rad2deg(np.arctan2(-vx, vy)) % 360
fits.PrimaryHDU(
header=hdu.header,
data=vv,
).writeto(
datadir / f"hh529_{s0}{s1}_vv{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=pa,
).writeto(
datadir / f"hh529_{s0}{s1}_pa{s2}.fits",
overwrite=True,
)
vx[abs(vv) > 12.0] = np.nan
vx[vv < 0.2] = np.nan
vx[(im < 1.0)] = np.nan
ny, nx = im.shape
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
fig, ax = plt.subplots(figsize=(10, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=0.8, vmax=1.5, origin="lower", cmap="gray_r")
step = slice(None, None, 4), slice(None, None, 4)
Q = ax.quiver(X[step], Y[step], vx[step], vy[step], pivot="middle",
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.5, width=1.5, minshaft=2.0, color="r", alpha=1.0)
Q.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
QK = ax.quiverkey(Q, 0.05, 0.2, 8.0,
"40 km/s", labelcolor="red",
fontproperties=dict(size=8))
QK.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
#ax.text(0.05, 0.95, "[O III] proper motions", color="w", transform=ax.transAxes)
ax.set(
xlabel="RA",
ylabel="Dec",
# xlim=[150, 1000],
# ylim=[200, 550],
)
fig.tight_layout(rect=(0.1, 0.1, 1.0, 1.0))
fig.savefig(figfile, dpi=1200)
- Next using 10 year intervals
- The problem this time is that the 2005 epoch is ACS, which lets in some [N II] emission
- This does not effect the high-ionization shocks much
- But it does produce some spurious apparent motions in the ionization fronts
- This can be spotted as opposite motions between 1996-2005 and 2005-2015
- The problem this time is that the 2005 epoch is ACS, which lets in some [N II] emission
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
figdir = Path("proper-motions")
datadir = figdir / "data"
s0 = "f656n_2005_2015"
figfile = figdir / f"hh529_{s0}_proper_motions.jpg"
hdu, = fits.open(datadir / "hh529_wfc3_f656n_2015_sharp_16.fits")
w = WCS(hdu)
im = hdu.data
#s1, s2 = "_sharp_16", "_sig20"
s1, s2 = "", "_sig10"
vx = fits.open(datadir / f"hh529_{s0}{s1}_vx{s2}.fits")[0].data
vy = fits.open(datadir / f"hh529_{s0}{s1}_vy{s2}.fits")[0].data
vv = np.hypot(vx, vy)
pa = np.rad2deg(np.arctan2(-vx, vy)) % 360
fits.PrimaryHDU(
header=hdu.header,
data=vv,
).writeto(
datadir / f"hh529_{s0}{s1}_vv{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=pa,
).writeto(
datadir / f"hh529_{s0}{s1}_pa{s2}.fits",
overwrite=True,
)
vx[abs(vv) > 6.0] = np.nan
vx[vv < 0.1] = np.nan
vx[(im < 1.0)] = np.nan
ny, nx = im.shape
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
fig, ax = plt.subplots(figsize=(10, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=0.8, vmax=1.5, origin="lower", cmap="gray_r")
step = slice(None, None, 4), slice(None, None, 4)
Q = ax.quiver(X[step], Y[step], vx[step], vy[step], pivot="middle",
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.25, width=1.5, minshaft=2.0, color="r", alpha=1.0)
Q.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
QK = ax.quiverkey(Q, 0.05, 0.2, 8.0,
"40 km/s", labelcolor="red",
fontproperties=dict(size=8))
QK.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
#ax.text(0.05, 0.95, "[O III] proper motions", color="w", transform=ax.transAxes)
ax.set(
xlabel="RA",
ylabel="Dec",
# xlim=[150, 1000],
# ylim=[200, 550],
)
fig.tight_layout(rect=(0.1, 0.1, 1.0, 1.0))
fig.savefig(figfile, dpi=1200)
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
figdir = Path("proper-motions")
datadir = figdir / "data"
s0 = "f656n_1996_2005"
figfile = figdir / f"hh529_{s0}_proper_motions.jpg"
hdu, = fits.open(datadir / "hh529_acs_f658n_2005_sharp_16.fits")
w = WCS(hdu)
im = hdu.data
#s1, s2 = "_sharp_16", "_sig20"
s1, s2 = "", "_sig10"
vx = fits.open(datadir / f"hh529_{s0}{s1}_vx{s2}.fits")[0].data
vy = fits.open(datadir / f"hh529_{s0}{s1}_vy{s2}.fits")[0].data
vv = np.hypot(vx, vy)
pa = np.rad2deg(np.arctan2(-vx, vy)) % 360
fits.PrimaryHDU(
header=hdu.header,
data=vv,
).writeto(
datadir / f"hh529_{s0}{s1}_vv{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=pa,
).writeto(
datadir / f"hh529_{s0}{s1}_pa{s2}.fits",
overwrite=True,
)
vx[abs(vv) > 6.0] = np.nan
vx[vv < 0.1] = np.nan
vx[(im < 1.0)] = np.nan
ny, nx = im.shape
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
fig, ax = plt.subplots(figsize=(10, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=0.8, vmax=1.5, origin="lower", cmap="gray_r")
step = slice(None, None, 4), slice(None, None, 4)
Q = ax.quiver(X[step], Y[step], vx[step], vy[step], pivot="middle",
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.25, width=1.5, minshaft=2.0, color="r", alpha=1.0)
Q.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
QK = ax.quiverkey(Q, 0.05, 0.2, 8.0,
"40 km/s", labelcolor="red",
fontproperties=dict(size=8))
QK.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
#ax.text(0.05, 0.95, "[O III] proper motions", color="w", transform=ax.transAxes)
ax.set(
xlabel="RA",
ylabel="Dec",
# xlim=[150, 1000],
# ylim=[200, 550],
)
fig.tight_layout(rect=(0.1, 0.1, 1.0, 1.0))
fig.savefig(figfile, dpi=1200)
- [X] Next thing to do: take an average of the two 10-year intervals
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("proper-motions/data")
hdu, = fits.open(datadir / "hh529_wfc3_f656n_2015_sharp_16.fits")
s1, s2 = "", "_sig10"
s_interval = "f656n_1996_2005"
vx1 = fits.open(datadir / f"hh529_{s_interval}{s1}_vx{s2}.fits")[0].data
vy1 = fits.open(datadir / f"hh529_{s_interval}{s1}_vy{s2}.fits")[0].data
m = vx1 == -20.0
vx1[m] = np.nan
vy1[m] = np.nan
s_interval = "f656n_2005_2015"
vx2 = fits.open(datadir / f"hh529_{s_interval}{s1}_vx{s2}.fits")[0].data
vy2 = fits.open(datadir / f"hh529_{s_interval}{s1}_vy{s2}.fits")[0].data
vxstack = np.stack([vx1, vx2])
vystack = np.stack([vy1, vy2])
vxm = np.nanmean(vxstack, axis=0)
vym = np.nanmean(vystack, axis=0)
vxs = np.nanstd(vxstack, axis=0)
vys = np.nanstd(vystack, axis=0)
vv = np.hypot(vxm, vym)
pa = np.rad2deg(np.arctan2(-vxm, vym)) % 360
sigma = np.hypot(vxs, vys) / vv
s_interval = "f656n_1996_2015_2stage"
fits.PrimaryHDU(
header=hdu.header,
data=vxm,
).writeto(
datadir / f"hh529_{s_interval}{s1}_vx_mean{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vym,
).writeto(
datadir / f"hh529_{s_interval}{s1}_vy_mean{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vxs,
).writeto(
datadir / f"hh529_{s_interval}{s1}_vx_sigma{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vys,
).writeto(
datadir / f"hh529_{s_interval}{s1}_vy_sigma{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=sigma,
).writeto(
datadir / f"hh529_{s_interval}{s1}_vv_sigrel{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=vv,
).writeto(
datadir / f"hh529_{s_interval}{s1}_vv{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=pa,
).writeto(
datadir / f"hh529_{s_interval}{s1}_pa{s2}.fits",
overwrite=True,
)
- Now do a better figure with the arrows.
- We can take into account the relative sigma
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
figdir = Path("proper-motions")
datadir = figdir / "data"
s0 = "f656n_1996_2015_2stage"
figfile = figdir / f"hh529_{s0}_proper_motions.jpg"
hdu, = fits.open(datadir / "hh529_wfc3_f656n_2015_sharp_16.fits")
w = WCS(hdu)
im = hdu.data
#s1, s2 = "_sharp_16", "_sig20"
s1, s2 = "", "_sig10"
vx = fits.open(datadir / f"hh529_{s0}{s1}_vx_mean{s2}.fits")[0].data
vy = fits.open(datadir / f"hh529_{s0}{s1}_vy_mean{s2}.fits")[0].data
sig = fits.open(datadir / f"hh529_{s0}{s1}_vv_sigrel{s2}.fits")[0].data
vv = np.hypot(vx, vy)
pa = np.rad2deg(np.arctan2(-vx, vy)) % 360
vx[abs(vv) > 6.0] = np.nan
vx[vv < 0.1] = np.nan
vx[(im < 1.0)] = np.nan
vx[(sig > 0.75)] = np.nan
ny, nx = im.shape
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
fig, ax = plt.subplots(figsize=(10, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=0.8, vmax=1.5, origin="lower", cmap="gray_r")
step = slice(2, None, 5), slice(3, None, 5)
Q = ax.quiver(X[step], Y[step], vx[step], vy[step], pivot="middle",
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.125, width=1.5, minshaft=2.0, color="orange", alpha=0.7)
Q.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='orange', capstyle="projecting", alpha=0.7
),
]
)
vx[(sig == 0.0)] = np.nan
Q = ax.quiver(X[step], Y[step], vx[step], vy[step], pivot="middle",
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.125, width=1.5, minshaft=2.0, color="r", alpha=1.0)
Q.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
QK = ax.quiverkey(Q, 0.05, 0.2, 4.0,
"40 km/s", labelcolor="red",
fontproperties=dict(size=8))
QK.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
#ax.text(0.05, 0.95, "[O III] proper motions", color="w", transform=ax.transAxes)
ax.set(
xlabel="RA",
ylabel="Dec",
xlim=[200, 1500],
ylim=[100, 900],
)
fig.tight_layout(rect=(0.1, 0.1, 1.0, 1.0))
fig.savefig(figfile, dpi=1200)
This is a lot easier than what we have been doing, and it might be useful
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("proper-motions/data")
hdu1, = fits.open(datadir / "hh529_wfpc2_f656n_1996.fits")
hdu2, = fits.open(datadir / "hh529_acs_f658n_2005.fits")
hdu3, = fits.open(datadir / "hh529_wfc3_f656n_2015.fits")
fits.PrimaryHDU(
header=hdu3.header,
data=hdu3.data/hdu2.data,
).writeto(
datadir / f"hh529_f656n_ratio_2015_2005.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu2.header,
data=hdu2.data/hdu1.data,
).writeto(
datadir / f"hh529_f656n_ratio_2005_1996.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu3.header,
data=hdu3.data/hdu1.data,
).writeto(
datadir / f"hh529_f656n_ratio_2015_1996.fits",
overwrite=True,
)
- Define regions in DS9
- proper-motions/HH529-samples.reg
- Ellipses are regions we measure PM
- Dashed boxes are for background brightness measurements
- [ ] Eventually I should repeat for HH 529-I, HH 529-East, and HH 1149, etc
- Plan is to calculate weighted average of V and PA for each region
- Weighting by BG-subtracted brightness
- We can do it separately for the 1996-2005 and 2005-2015 intervals, where available
- This only works because the WCS is the same for all the images
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.table import Table
import regions as rg
regiondir = Path("proper-motions")
datadir = regiondir / "data"
regions = rg.read_ds9(regiondir / "HH529-samples.reg")
# Split into BG and source regions
source_regions = []
bg_regions = {}
for region in regions:
label = region.meta["label"]
if "BG" in label:
# BG samples are grouped into II, III, etc
bg_group = label.split()[0]
if bg_group not in bg_regions:
bg_regions[bg_group] = []
bg_regions[bg_group].append(region)
else:
source_regions.append(region)
hdu2015, = fits.open(datadir / "hh529_wfc3_f656n_2015.fits")
w = WCS(hdu2015.header)
badmask = ~np.isfinite(hdu2015.data)
hdu2015.data[badmask] = 0.0
data_images = {
# "VV 96->05": fits.open(datadir / "hh529_f656n_1996_2005_vv_sig10.fits")[0].data,
# "PA 96->05": fits.open(datadir / "hh529_f656n_1996_2005_pa_sig10.fits")[0].data,
# "VV 05->15": fits.open(datadir / "hh529_f656n_2005_2015_vv_sig10.fits")[0].data,
# "PA 05->15": fits.open(datadir / "hh529_f656n_2005_2015_pa_sig10.fits")[0].data,
"VV": fits.open(datadir / "hh529_f656n_1996_2015_2stage_vv_sig10.fits")[0].data,
"PA": fits.open(datadir / "hh529_f656n_1996_2015_2stage_pa_sig10.fits")[0].data,
}
# Find mean brightness BG for each group
bg_brightness = {}
for bg_group, bgrs in bg_regions.items():
weights = np.zeros_like(hdu2015.data)
for r in bgrs:
weights += r.to_pixel(w).to_mask().to_image(hdu2015.data.shape)
bg_brightness[bg_group] = np.average(hdu2015.data, weights=weights)
# Now do measurements on each source region
all_results = []
for region in source_regions:
label = region.meta["label"]
group = label.split()[0]
bg = bg_brightness[group]
results = {"region": label}
mask = region.to_pixel(w).to_mask()
mw = mask.to_image(hdu2015.data.shape).astype(bool)
weights = hdu2015.data[mw] - bg
results["Bright"] = np.round(np.mean(weights), 2)
results["Bright sig"] = np.round(np.std(weights), 2)
results["Bright BG"] = np.round(bg, 2)
for data_label, data in data_images.items():
m = mask.to_image(data.shape).astype(bool)
mean = np.average(data[m], weights=weights)
sigma = np.sqrt(np.average((data[m] - mean)**2, weights=weights))
results[data_label] = np.round(mean, 2)
results[f"{data_label} sig"] = np.round(sigma, 2)
all_results.append(results)
#
# List of dicts to astropy.table.Table
#
# It would be cleaner to just write Table(rows=all_results) but for
# some reason that does alphabetical sorting on column names, which we
# don't want. This way preserves the column order.
table = Table(
names=all_results[0].keys(),
rows=[_.values() for _ in all_results],
)
#
# astropy.table.Table to a list of lists, which babel will format as
# an org table
#
display = [list(table.colnames), None] + [list(row) for row in table]
region | Bright | Bright sig | Bright BG | VV | VV sig | PA | PA sig |
---|---|---|---|---|---|---|---|
III knot | 13.47 | 1.21 | 11.85 | 2.67 | 0.03 | 129.84 | 0.9 |
III knot N | 11.2 | 1.06 | 11.85 | 2.98 | 0.17 | 104.75 | 3.3 |
III knot S | 7.57 | 0.82 | 11.85 | 3.03 | 0.23 | 124.9 | 4.24 |
III apex | 9.7 | 1.3 | 11.85 | 3.56 | 0.1 | 106.79 | 6.79 |
III apex N | 5.49 | 1.69 | 11.85 | 3.33 | 0.34 | 90.18 | 3.18 |
III apex S | 6.18 | 1.61 | 11.85 | 3.16 | 0.14 | 117.18 | 3.12 |
III lower wing | 4.8 | 0.88 | 11.85 | 3.92 | 0.19 | 136.86 | 4.63 |
III lower wing S | 3.16 | 0.87 | 11.85 | 3.15 | 0.12 | 151.74 | 10.23 |
II a | 3.3 | 0.46 | 13.2 | 2.09 | 0.93 | 117.04 | 58.35 |
II c | 2.91 | 0.42 | 13.2 | 3.52 | 0.86 | 151.16 | 87.38 |
II b | 4.62 | 1.05 | 13.2 | 2.58 | 0.51 | 107.14 | 4.44 |
PA | PA + 90 | ||||||
129.84 | 219.84 | ||||||
104.75 | 194.75 | ||||||
124.9 | 214.90 | ||||||
106.79 | 196.79 | ||||||
90.18 | 180.18 | ||||||
117.18 | 207.18 | ||||||
136.86 | 226.86 | ||||||
151.74 | 241.74 | ||||||
117.04 | 207.04 | ||||||
151.16 | 241.16 | ||||||
107.14 | 197.14 |
I have now renamed these features for the figure.
region | new name | Bright | Bright sig | Bright BG | VV | VV sig | PA | PA sig | Rel Bright |
---|---|---|---|---|---|---|---|---|---|
III knot | HH 529 III b2 | 13.47 | 1.21 | 11.85 | 2.67 | 0.03 | 129.84 | 0.9 | 1.14 +/- 0.10 |
III knot N | HH 529 III b1 | 11.2 | 1.06 | 11.85 | 2.98 | 0.17 | 104.75 | 3.3 | 0.95 +/- 0.09 |
III knot S | HH 529 III b3 | 7.57 | 0.82 | 11.85 | 3.03 | 0.23 | 124.9 | 4.24 | 0.64 +/- 0.07 |
III apex | HH 529 III a2 | 9.7 | 1.3 | 11.85 | 3.56 | 0.1 | 106.79 | 6.79 | 0.82 +/- 0.11 |
III apex N | HH 529 III a1 | 5.49 | 1.69 | 11.85 | 3.33 | 0.34 | 90.18 | 3.18 | 0.46 +/- 0.14 |
III apex S | HH 529 III a3 | 6.18 | 1.61 | 11.85 | 3.16 | 0.14 | 117.18 | 3.12 | 0.52 +/- 0.14 |
III lower wing | HH 529 III a4 | 4.8 | 0.88 | 11.85 | 3.92 | 0.19 | 136.86 | 4.63 | 0.41 +/- 0.07 |
III lower wing S | HH 529 III a5 | 3.16 | 0.87 | 11.85 | 3.15 | 0.12 | 151.74 | 10.23 | 0.27 +/- 0.07 |
II a | HH 529 II a | 3.3 | 0.46 | 13.2 | 2.09 | 0.93 | 117.04 | 58.35 | 0.25 +/- 0.03 |
II c | HH 529 II c | 2.91 | 0.42 | 13.2 | 3.52 | 0.86 | 151.16 | 87.38 | 0.22 +/- 0.03 |
II b | HH 529 II b | 4.62 | 1.05 | 13.2 | 2.58 | 0.51 | 107.14 | 4.44 | 0.35 +/- 0.08 |
And a version for exporting to latex
Region | VV | PA | Rel Bright |
---|---|---|---|
HH 529 III b1 | 30 ± 2 | 105 ± 3 | 0.95 ± 0.09 |
HH 529 III b2 | 27 ± 1 | 130 ± 1 | 1.14 ± 0.10 |
HH 529 III b3 | 30 ± 2 | 125 ± 4 | 0.64 ± 0.07 |
HH 529 III a1 | 33 ± 3 | 90 ± 3 | 0.46 ± 0.14 |
HH 529 III a2 | 36 ± 1 | 107 ± 7 | 0.82 ± 0.11 |
HH 529 III a3 | 32 ± 1 | 117 ± 3 | 0.52 ± 0.14 |
HH 529 III a4 | 39 ± 2 | 137 ± 5 | 0.41 ± 0.07 |
HH 529 III a5 | 32 ± 1 | 152 ± 10 | 0.27 ± 0.07 |
HH 529 II a | 21 ± 9 | 117 ± 58 | 0.25 ± 0.03 |
HH 529 II b | 26 ± 5 | 107 ± 4 | 0.35 ± 0.08 |
HH 529 II c | 35 ± 9 | 151 ± 87 | 0.22 ± 0.03 |
- HH 529 III a2
- V = 36 +/- 1 km/s
- PA = 107 +/- 7 deg
- contrast = 0.82 +/- 0.11
- HH 529 III b1
- V = 30 +/- 2 km/s
- PA = 105 +/- 3 deg
- contrast = 0.95 +/- 0.09
- HH 529 II a
- V = 21 +/- 9 km/s
- PA = 120 +/- 50 deg
- contrast = 0.25 +/- 0.03
- HH 529 II b
- V = 26 +/- 5 km/s
- PA = 107 +/- 4 km/s
- contrast = 0.35 +/- 0.08
- HH 529 II is the most straightforward, since the components are narrow
- Suggesting either bullets are internal working surfaces
- II a has an extent on the [Fe III] PV array of Fig 2e of [-40, -20] km/s
- Thermal broadening should be negligible (3 km/s) and instrumental width is 6.5 km/s
- So the corrected FWZI is sqrt(20**2 - 6.5**2 - 3**2) = 18.5 km/s, which should be the working surface speed relative to the upstream gas
- The average Vhel is -30 km/s, which gives Vr = -57 km/s
- Proper motion gives Vt = 21 +/- 9 km/s, so V = 61 km/s and inc = arctan2(-57, 21 +/- 9) = -70 +/- 8 deg
- II b has Vhel in range [-25, -45]
- V_WS - V_us = 18.5
- Average Vhel = -35 => Vr = -62
- Vt = 26 +/- 5 => V = -67 km/s and inc = -67 +/- 4
- HH 529 III is a combination of a2 and b1 in the spectrum
- a2 is part of the bow shock
- whereas b1 may be something different
- To start with, analyse them as a single entity
- the fact that it goes down to near-systemic velocities suggests it is a terminal bow
- highest velocity is -45 km/s => Vr = -72 km/s
- Average Vt = 33 km/s
- If we used
- Vmax = V 0.5 (1 + sin i)
- Vt = V cos i
- Vmax/Vt = r = 0.5 (seci + tan i)
- 2r - tan i = sec i
- 4r^2 + tan^2 i - 4 r tan i = 1 + tan^2 i
- tan i = (4r^2 - 1)/4r = r - 1/4r
- r = 72/(33 +/- 3) = 2.18 +/- 0.20 => r - 1/4r = 2.07 +/- 0.20 => inc = 64 +/- 2 deg
- V = 72 / 0.5 (1 + sin(64 +/- 2)) = 75.8 +/- 0.6 km/s
- Alternative V = (33 +/- 3) / cos(64 +/- 2) = 75.3 +/- 0.7 - close enough
- But this means that we should get a most red-shifted velocity of Vr = V 0.5 (sin i - 1) = -0.05 V = 4 km/s
- This would mean Vhel = +30 km/s
- This would be hard to see - maybe there is evidence in [O I]
- Alternatively, we could suppose that the most redshifted velocity is actually Vhel = 0
- Suggested at least for III b1 from the spectrum
- This would give V_WS - V_us = 45 km/s
- Average Vhel = -23 km/s => Vr = -50 km/s
- Vt = 30 +/- 2 km/s => V = 58 +/- 3 km/s
- Implying upstream V = 13 km/s
- inc = arctan2(-50, 30) = -59 deg
- Going back,
- 0.5 (1 + sin(i)) V = 0.93 45 = 42 km/s
- 0.5 (-1 + sin(i)) V = -0.07 45 = -3 km/s
- This implies that the rest frame of upstream is Vr = +3 km/s, which is -30 km/s wrt systemic
- I don’t quite understand this, but lets carry on
- If that was at the same angle (60 deg), then upstream V = 30 / sin(60) = 35 km/s
- So Vt = 35 cos(60) = 17.5 km/s
- This would be a prediction for III a2 but we actually have 36 km/s!
- Center on the star that is superimposed on HH 529-III
- CRPIX1, CRPIX2 = 4060.407, 4230.3817
- CRVAL1, CRVAL2 = 83.820644, -5.3998766
- Extract a region around HH529 and save both epochs with the same image pixels
- Use
reproject
for this
- Use
- First extract a more manageable window of the 2005 image
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
datadir = Path("../Orion-HH-data/ACS")
# Corners of window to extract from image
xmin, xmax, ymin, ymax = 7324, 10724, 2787, 5717
rect = slice(ymin, ymax), slice(xmin, xmax)
# Big file from Robberto 2005 program
hdu0 = fits.open(datadir / "hlsp_orion_hst_acs_strip0l_f775w_v1_drz.fits")["SCI"]
# Fix stupid SIP keywords that shouldn't be there
bad_kwds = [_ for _ in hdu0.header if _[:2] in ("A_", "B_")]
for kwd in bad_kwds:
hdu0.header.remove(kwd)
# Now we are ready to create the WCS
w = WCS(hdu0.header)
outdir = Path("proper-motions/data")
fits.PrimaryHDU(
data=hdu0.data[rect],
header=w.slice(rect).to_header(),
).writeto(
outdir / "hh529_acs_f775w_2005.fits",
overwrite=True,
)
- And use FLCT to calculate proper motions
import sys
sys.path.append("../../Source/pyflct")
import os
from pathlib import Path
import numpy as np
import pyflct
from astropy.io import fits
from astropy.wcs import WCS
os.environ["PATH"] = ":".join((os.environ["PATH"], "/Users/will/.local/bin"))
datadir = Path("proper-motions/data")
h1, = fits.open(datadir / "hh529_acs_f775w_2005.fits")
h2, = fits.open(datadir / "hh529_acs_f775w_2015.fits")
# Mask out pixels that are star/diffraction/bleed
thresh = 80.0
m = (h1.data > thresh) | (h2.data > thresh)
# Set first epoch equal to second there
h1.data[m] = h2.data[m]
data_scale = h1.data.max()
sigma = 10.0
vx, vy, vm = pyflct.flct(h1.data.T, h2.data.T,
deltat=1.0, deltas=1.0, sigma=sigma,
thresh=0.1/data_scale,
)
vx[vm==0] = np.nan
vy[vm==0] = np.nan
fits.PrimaryHDU(
header=h1.header,
data=vx.T,
).writeto(datadir / f"hh529_acs_f775w_vx_sig{int(sigma):02d}.fits", overwrite=True)
fits.PrimaryHDU(
header=h1.header,
data=vy.T,
).writeto(datadir / f"hh529_acs_f775w_vy_sig{int(sigma):02d}.fits", overwrite=True)
- Make an high-pass filtered version
python /Users/will/Dropbox/Teresa-Turtle/scripts/smooth-image.py hh529_acs_f775w_2015.fits -w 16
python /Users/will/Dropbox/Teresa-Turtle/scripts/smooth-image.py hh529_acs_f775w_2005.fits -w 16
- And now do FLCT on the high-pass filtered images
import sys
sys.path.append("../../Source/pyflct")
import os
from pathlib import Path
import numpy as np
import pyflct
from astropy.io import fits
from astropy.wcs import WCS
os.environ["PATH"] = ":".join((os.environ["PATH"], "/Users/will/.local/bin"))
datadir = Path("proper-motions/data")
h1, = fits.open(datadir / "hh529_acs_f775w_2005_sharp_16.fits")
h2, = fits.open(datadir / "hh529_acs_f775w_2015_sharp_16.fits")
# Mask out pixels that are star/diffraction/bleed
thresh = 1.6
m = (h1.data > thresh) | (h2.data > thresh)
# Set first epoch equal to second there
h1.data[m] = h2.data[m]
data_scale = h1.data.max()
sigma = 20.0
vx, vy, vm = pyflct.flct(h1.data.T, h2.data.T,
deltat=1.0, deltas=1.0, sigma=sigma,
thresh=0.1/data_scale,
)
vx[vm==0] = np.nan
vy[vm==0] = np.nan
fits.PrimaryHDU(
header=h1.header,
data=vx.T,
).writeto(datadir / f"hh529_acs_f775w_sharp_16_vx_sig{int(sigma):02d}.fits", overwrite=True)
fits.PrimaryHDU(
header=h1.header,
data=vy.T,
).writeto(datadir / f"hh529_acs_f775w_sharp_16_vy_sig{int(sigma):02d}.fits", overwrite=True)
- Plot the proper motions
- Convert to physical velocities
- ACS pixel scale is 1.38888888888896E-05 deg = 0.05 arcsec = 50 mas
- Time span
- MJD = 57081.4732 = 2015-02-28 = 2015.16
- 2005-04 = 2005.25
- dt = 2015.16 - 2005.25 = 9.91 years
- 1 pixel in 9.91 yrs is 50 / 9.91 = 5.045 mas/yr
- At a distance of 417 pc, this is 5.045 417 au / 1000 yr km = 9.973 km/s
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.patheffects as path_effects
import seaborn as sns
from astropy.io import fits
from astropy.wcs import WCS
figdir = Path("proper-motions")
datadir = figdir / "data"
figfile = figdir / "hh529_acs_f775w_proper_motions.jpg"
hdu, = fits.open(datadir / "hh529_acs_f775w_2015_sharp_16.fits")
w = WCS(hdu)
im = hdu.data
#s1, s2 = "_sharp_16", "_sig20"
s1, s2 = "", "_sig10"
vx = fits.open(datadir / f"hh529_acs_f775w{s1}_vx{s2}.fits")[0].data
vy = fits.open(datadir / f"hh529_acs_f775w{s1}_vy{s2}.fits")[0].data
vv = np.hypot(vx, vy)
pa = np.rad2deg(np.arctan2(-vx, vy)) % 360
fits.PrimaryHDU(
header=hdu.header,
data=vv,
).writeto(
datadir / f"hh529_acs_f775w{s1}_vv{s2}.fits",
overwrite=True,
)
fits.PrimaryHDU(
header=hdu.header,
data=pa,
).writeto(
datadir / f"hh529_acs_f775w{s1}_pa{s2}.fits",
overwrite=True,
)
vx[abs(vv) > 5.0] = np.nan
vx[vv < 0.1] = np.nan
vx[(im < 1.0)] = np.nan
ny, nx = im.shape
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
fig, ax = plt.subplots(figsize=(10, 4), subplot_kw=dict(projection=w))
ax.imshow(im, vmin=0.8, vmax=1.5, origin="lower", cmap="gray_r")
# A faint version for the faint emission
step = slice(None, None, 7), slice(None, None, 7)
Q = ax.quiver(X[step], Y[step], vx[step], vy[step], pivot="middle",
headwidth=4, headlength=6, minlength=0.01,
units="xy", scale=0.2, width=1.5, minshaft=2.0, color="r", alpha=1.0)
Q.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
QK = ax.quiverkey(Q, 0.05, 0.2, 4.0,
"40 km/s", labelcolor="red",
fontproperties=dict(size=8))
QK.set_path_effects(
[
path_effects.PathPatchEffect(
edgecolor='white', linewidth=0.1, facecolor='red', capstyle="projecting",
),
]
)
#ax.text(0.05, 0.95, "[O III] proper motions", color="w", transform=ax.transAxes)
ax.set(
xlabel="RA",
ylabel="Dec",
xlim=[150, 1000],
ylim=[200, 550],
)
fig.tight_layout(rect=(0.1, 0.1, 1.0, 1.0))
fig.savefig(figfile, dpi=1200)
ORDER="HH-529-ACS-GO10246-2005-04.png HH-529-ACS-GO13826-2015-10.png"
convert $(for a in "$ORDER"; do printf -- "-delay 25 %s " $a; done; ) HH529-motions.gif
ORDER="HH529-III-bowshock-frame-2005.png HH529-III-bowshock-frame-2015.png"
convert $(for a in "$ORDER"; do printf -- "-delay 25 %s " $a; done; ) HH529-III-bowshock-frame.gif
- HH 529-III
- PA = 98 +/- 10 deg
- V_t = 35 +/- 2 km/s
- HH 529-II
- PA = 105 +/- 13
- V_t = 28 +/- 5
- HH 529-I
- PA = 83 +/- 50
- V_t = 25 +/- 15
- HH 529-east
- V_t = 28 +/- 6
- PA = 95 +/- 13
- HH 529-1.5
- V_t = 13 +/- 3
- PA = 61 +/- 20
- HH 529-0
- V_t = 25 +/- 7
- PA = 96 +/- 60
- East shock
- PA = 170 +/- 8
- V_t = 9 +/- 2 km/s
- HH 529-Ram
- V_t = 25 +/- 10
- PA = 130 +/- 70
- HH 998
- PA = 224 +/- 24
- V_t = 54 +/- 27
- Move the CRPIX according to the motion
- -35 sin(98) / 10 = -3.466 pix
- +35 cos(98) / 10 = -0.487 pix
- [676.0 - 3.466, -37.0 - 0.487] = [672.534, -37.487]
- Turns out that the bowshock frame is a bit different
- CRPIX = [671.534, -37.987]
- So shift is [671.534, -37.987] - [676.0, -37.0] = [-4.466, -0.987]
- 4.5738 pixels = 46 km/s
- PA = arctan2(4.466, -0.987) = 102.4 deg
- Why is this different from the results from FLCT?
- Because there is sub-structure to HH 529-III with different velocities
- The bow shock is moving the fastest
- There is a bright knot, with its own half mini-bow to north, which is moving slower
- Objects that I name here (not all are exactly new)
- HH 529-1.5
- This is what ODell:2015a calls “several west-moving intervening high-ionization knots” and “a faint series of moving knots centered 8.3 arcsec east of COUP 769. The axis of these knots is about 82°” (sec 3.3.1.1)
- They are aligned with HH 529-east
- HH 529-0
- This is just above the E end of the “East shock”, about 2 arcsec W of the HH 998 knot
- HH 529-1.5
- Epoch: 2012
- Published in Forbrich:2016a
- Data is ~/Dropbox/RadioProplyds/JanData/concat_nouvra_inner4k.fits
- In principle, these should be better absolute coordinates
- They don’t quite align with Robberto 2005
- I had to subtract (1.0, -0.5) pixels from the CRPIX of the radio data
- But a better solution would be to shift the optical coordinates
This is now in a separate file because it was too big:
- Remove the labels
- [X] Query about position
- Corrected in email from Eduardo
- Filters to use
- I had thought about [O I]/[O II]/[O III]
- But the [O II] is from 2015 while [O I] is 1996, so there will be motions
- Could use all from 1996
- Or all from 2015
- I had thought about [O I]/[O II]/[O III]
- Dry run with DS9
- I am comparing two different sets (the Bally ones have too weird a set of filters)
xpaset -p RGB rgb red xpaget RGB file xpaget RGB scale limits xpaset -p RGB rgb green xpaget RGB file xpaget RGB scale limits xpaset -p RGB rgb blue xpaget RGB file xpaget RGB scale limits
The WFPC2 1996 set
/Users/will/Work/BobPC/2002/final502-radec.fits 130 750 /Users/will/Work/BobPC/2002/final658-radec.fits 200 2200 /Users/will/Work/BobPC/2002/final673-radec.fits 70 500
The ACS 2005 set:
/Users/will/Work/OrionTreasury/acs/hlsp_orion_hst_acs_strip0l_f658n_v1_drz.fits[SCI] 20 140 /Users/will/Work/OrionTreasury/acs/hlsp_orion_hst_acs_strip0l_f555w_v1_drz.fits[SCI] 20 130 /Users/will/Work/OrionTreasury/acs/hlsp_orion_hst_acs_strip0l_f435w_v1_drz.fits[SCI] 7 35
- I am comparing two different sets (the Bally ones have too weird a set of filters)
final502-radec | 130 | 750 |
final658-radec | 200 | 2200 |
final673-radec | 70 | 500 |
1354 | 1369 | 760 | 600 |
from pathlib import Path
import matplotlib.pyplot as plt
from astropy.visualization import make_lupton_rgb
from astropy.io import fits
from astropy.wcs import WCS
import regions as rg
datadir = Path(DIR)
# Unpack the channel info from the table
[rf, r1, r2], [gf, g1, g2], [bf, b1, b2] = TAB
def load_and_scale_image(fn, v1, v2, ihdu=HDU):
data = fits.open(datadir / f"{fn}.fits")[ihdu].data
return (data - v1) / (v2 - v1)
w = WCS(fits.open(datadir / f"{rf}.fits")[HDU].header)
image_r = load_and_scale_image(rf, r1, r2)
image_g = load_and_scale_image(gf, g1, g2)
image_b = load_and_scale_image(bf, b1, b2)
image = make_lupton_rgb(image_r, image_g, image_b, stretch=STRETCH, Q=Q)
[x0, y0, dx, dy], = WINDOW
figfile = f"rgb-lupton-hh{HH}-{SUFFIX}.jpg"
fig, ax = plt.subplots(figsize=(6.2, 5), subplot_kw=dict(projection=w))
ax.imshow(image)
regions = rg.read_ds9(f"hh{HH}-slit.reg")
for skyregion in regions:
pixelregion = skyregion.to_pixel(w)
pixelregion.plot(ax=ax, color="w", lw=0.7, ls="--")
ax.set(
xlabel="RA (J2000)",
ylabel="Dec (J2000)",
xlim=[x0-dx/2, x0+dx/2],
ylim=[y0-dy/2, y0+dy/2],
)
fig.tight_layout(rect=[0.15, 0.07, 1.0, 1.0])
fig.savefig(figfile, dpi=300)
hlsp_orion_hst_acs_strip0l_f658n_v1_drz | 20 | 140 |
hlsp_orion_hst_acs_strip0l_f555w_v1_drz | 20 | 130 |
hlsp_orion_hst_acs_strip0l_f435w_v1_drz | 7 | 35 |
8899 | 4312 | 1500 | 1165 |
final502-radec | 130 | 750 |
final658-radec | 200 | 2200 |
final656-radec | 600 | 2500 |
1354 | 1369 | 1140 | 900 |
- Coordinates of slit
- RA(J2000)=05h35𝑚22𝑠.72,
- DEC(J2000)=−05 25 20.42
- Coordinates of zero-point of section 5
- RA(J2000)=05 35 22.63, DEC(J2000)=−05 25 21.86
- Correction 05 35 22.82 −05 25 21.76
- Second correction from Eduardo [2021-01-12 Tue]
- RA: 05 35 22.81, DEC:−05 25 21.86
- 8 pixels SE of the slit center
- 0.246 arcsec/pixel => 2 arcsec from slit center, so 3 arcsec from slit end.
- Original observing proposal
- 05 35 22.74
- -05 25 19.47
final502-radec | 30 | 250 |
final658-radec | 100 | 2000 |
final656-radec | 200 | 2000 |
400 | 630 | 800 | 630 |
- 2005 from Robberto is the reference frame
- 2015-01 Bally is aligned more or less by
trap-south-align-robberto.wcs
- Still seems off by a fraction of a pixel, but not important
- For the 1993 mosaic, we need to make quite a big correction
- Use the unsaturated star in between HH 203 and 204
- 5:35:22.4015 -5:25:09.449
- 83.8433551 -5.4192964
- Pixels in 1993 mosaics: 388.39887 567.14165
- This alignment is saved in mosaic-1996-HH204-align-robberto.wcs
- Rather than re-use the HH 529 program, I should just write another one
- I can use the WCS adjustments described above
- It looks like the
.wcs
file can be interpreted as TOML - Except that TOML has decided that “0.” is not a valid float, so we need to fix that
import toml def read_wcs(wcsfile: str) -> dict: "Read WCS (key,value) pairs from DS9's *.wcs files" with open(wcsfile) as stream: return toml.loads(_fixups(stream.read())) def _fixups(s: str) -> str: s = s.replace("0.\n", "0.0\n") return s if __name__ == "__main__": data = read_wcs("mosaic-1996-HH204-align-robberto.wcs") print(data)
Test out the parser:
python read_wcs.py
- It looks like the
- Here is the new program. It is almost identical to the previous one, except for applying the WCS correction.
- We also have changed the slit length to be
- 38 pixels times 0.246 arcsec = 9.348 arcsec
- And cut off the LHS with 244-440 since it is not on the large-scale image
570 | 680 | 693 | 693 |
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import make_lupton_rgb
from astropy.io import fits
from astropy.wcs import WCS
import regions as rg
from read_wcs import read_wcs
import seaborn as sns
datadir = Path(DIR)
# Unpack the channel info from the table
[rf, r1, r2], [gf, g1, g2], [bf, b1, b2] = TAB
def load_and_scale_image(fn, v1, v2, ihdu=HDU):
data = fits.open(datadir / f"{fn}.fits")[ihdu].data
return (data - v1) / (v2 - v1)
hdr = fits.open(datadir / f"{rf}.fits")[HDU].header
hdr.update(**read_wcs("mosaic-1996-HH204-align-robberto.wcs"))
w = WCS(hdr)
image_r = load_and_scale_image(rf, r1, r2)
image_g = load_and_scale_image(gf, g1, g2)
image_b = load_and_scale_image(bf, b1, b2)
image = make_lupton_rgb(image_r, image_g, image_b, stretch=STRETCH, Q=Q)
[x0, y0, dx, dy], = WINDOW
figfile = f"rgb-lupton-hh{HH}-{SUFFIX}.jpg"
fig, ax = plt.subplots(figsize=(5, 5), subplot_kw=dict(projection=w))
sns.set_color_codes("bright")
ax.imshow(image)
# Add contours of the high-velocity h alpha
kpnodir = Path("../../Work/BobKPNO/2004")
oiii_hdu = fits.open(kpnodir / "vcube.oiii.fits")[0]
oiii_hdr = oiii_hdu.header
oiii_hdr.update(**read_wcs(kpnodir / "../DOH-radec-vel-oiii.wcs"))
nii_hdu = fits.open(kpnodir / "vcube.nii.fits")[0]
nii_hdr = nii_hdu.header
nii_hdr.update(**read_wcs(kpnodir / "../DOH-radec-vel-oiii.wcs"))
wkpno = WCS(oiii_hdr).celestial
# centered on fits channel 13 of oiii
vmap_blue = np.mean(oiii_hdu.data[12:13, :, :], axis=0)
# centered on fits channel 17 of nii
vmap_red = np.mean(nii_hdu.data[16:17, :, :], axis=0)
# subtract continuum
vmap_blue -= np.mean(oiii_hdu.data[0:5, :, :], axis=0)
vmap_red -= np.mean(nii_hdu.data[0:5, :, :], axis=0)
mstar = np.mean(oiii_hdu.data[0:5, :, :], axis=0) > 5.0
vmap_blue[mstar] = 0.0
vmap_red[mstar] = 0.0
# scaling to match contour levels
vmap_blue *= 3
vmap_red *= 2
ax.contour(vmap_blue,
transform=ax.get_transform(wkpno),
levels=[40, 160],
colors="c",
linewidths=[1.0, 2.0],
)
ax.contour(vmap_red,
transform=ax.get_transform(wkpno),
levels=[250, 1000],
# cmap="BuPu",
colors="y",
linewidths=[1.0, 2.0],
)
# centered on fits channel 17 of nii
vmap_knot = np.mean(nii_hdu.data[20:23, :, :], axis=0)
vmap_knot -= np.mean(nii_hdu.data[0:5, :, :], axis=0)
m = oiii_hdu.data[32, :, :] < 500.0
vmap_knot[m] = 0.0
vmap_knot[mstar] = 0.0
vmap_knot *= 0.9
ax.contour(vmap_knot,
transform=ax.get_transform(wkpno),
levels=[130, 160],
# cmap="BuPu",
colors="w",
linewidths=[1.0, 2.0],
)
regions = rg.read_ds9(f"hh{HH}-slit.reg")
for skyregion in regions:
pixelregion = skyregion.to_pixel(w)
pixelregion.plot(ax=ax, color="w", lw=0.7, ls="--")
ax.set(
xlabel="RA (J2000)",
ylabel="Dec (J2000)",
xlim=[x0-dx/2, x0+dx/2],
ylim=[y0-dy/2, y0+dy/2],
)
fig.tight_layout(rect=[0.18, 0.05, 1.0, 1.0])
fig.savefig(figfile, dpi=300)
480 | 610 | 420 | 420 |
final502-radec | 30 | 180 |
final658-radec | 100 | 2000 |
final656-radec | 200 | 2000 |
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import make_lupton_rgb
from astropy.io import fits
from astropy.wcs import WCS
import regions as rg
from read_wcs import read_wcs
import seaborn as sns
datadir = Path(DIR)
# Unpack the channel info from the table
[rf, r1, r2], [gf, g1, g2], [bf, b1, b2] = TAB
def load_and_scale_image(fn, v1, v2, ihdu=HDU):
data = fits.open(datadir / f"{fn}.fits")[ihdu].data
return (data - v1) / (v2 - v1)
hdr = fits.open(datadir / f"{rf}.fits")[HDU].header
hdr.update(**read_wcs("mosaic-1996-HH204-align-robberto.wcs"))
w = WCS(hdr)
image_r = load_and_scale_image(rf, r1, r2)
image_g = load_and_scale_image(gf, g1, g2)
image_b = load_and_scale_image(bf, b1, b2)
image = make_lupton_rgb(image_r, image_g, image_b, stretch=STRETCH, Q=Q)
[x0, y0, dx, dy], = WINDOW
figfile = f"rgb-lupton-hh{HH}-{SUFFIX}.jpg"
fig, ax = plt.subplots(figsize=(5, 5), subplot_kw=dict(projection=w))
sns.set_color_codes("bright")
ax.imshow(image)
regions = rg.read_ds9(f"hh{HH}-slit.reg")
for skyregion in regions:
pixelregion = skyregion.to_pixel(w)
pixelregion.plot(ax=ax, color="w", lw=0.7, ls="--")
ax.set(
xlabel="RA (J2000)",
ylabel="Dec (J2000)",
xlim=[x0-dx/2, x0+dx/2],
ylim=[y0-dy/2, y0+dy/2],
)
fig.tight_layout(rect=[0.18, 0.05, 1.0, 1.0])
fig.savefig(figfile, dpi=300)
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import make_lupton_rgb
from astropy.io import fits
from astropy.wcs import WCS
import regions as rg
from read_wcs import read_wcs
import seaborn as sns
kpnodir = Path("../../Work/BobKPNO/2004")
oiii_hdu = fits.open(kpnodir / "vcube.oiii.fits")[0]
oiii_hdr = oiii_hdu.header
oiii_hdr.update(**read_wcs(kpnodir / "../DOH-radec-vel-oiii.wcs"))
nii_hdu = fits.open(kpnodir / "vcube.nii.fits")[0]
nii_hdr = nii_hdu.header
nii_hdr.update(**read_wcs(kpnodir / "../DOH-radec-vel-oiii.wcs"))
wkpno = WCS(oiii_hdr).celestial
image_r = np.mean(oiii_hdu.data[22:25, :, :], axis=0)
image_r /= 200.0
image_g = np.mean(oiii_hdu.data[17:23, :, :], axis=0)
image_g /= 100.0
image_b = np.mean(oiii_hdu.data[11:18, :, :], axis=0)
image_b /= 80.0
image1 = make_lupton_rgb(image_r, image_g, image_b, stretch=STRETCH, Q=Q)
image_r = np.mean(nii_hdu.data[22:26, :, :], axis=0)
image_r /= 1500.0
image_g = np.mean(nii_hdu.data[14:23, :, :], axis=0)
image_g /= 2000.0
image_b = np.mean(nii_hdu.data[10:15, :, :], axis=0)
image_b /= 3000.0
image2 = make_lupton_rgb(image_r, image_g, image_b, stretch=STRETCH, Q=Q)
figfile = f"rgb-lupton-channels.jpg"
fig, [ax1, ax2] = plt.subplots(1, 2,
figsize=(8, 5),
subplot_kw=dict(projection=wkpno))
sns.set_color_codes("bright")
ax1.imshow(image1)
ax2.imshow(image2)
ax1.set(
xlabel="RA (J2000)",
ylabel="Dec (J2000)",
)
ax2.set(
xlabel="RA (J2000)",
ylabel=" ",
)
#fig.tight_layout(rect=[0.18, 0.05, 1.0, 1.0])
fig.savefig(figfile, dpi=300)
- Velocities for the RGB images:
line chan k1 k2 kmean adjust V_hel V_ONC dV [O III] B 11 18 14 0 -46 -71 16 [O III] G 17 23 19.5 0 -24. -49. 14 [O III] R 22 25 23 0 -10 -35 8 [N II] B 10 15 12 4 -50 -75 12 [N II] G 14 23 18 4 -26 -51 20 [N II] R 22 26 23.5 4 -4. -29. 10 [N II] neb B 23 28 25 4 2 -23 12 [N II] neb G 27 32 29 4 18 -7 12 [N II] neb R 32 37 34 4 38 13 12
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
import regions as rg
import seaborn as sns
datadir = Path("../dib-scatter-hii/data/orion-hawki")
h2_hdu = fits.open(datadir / "Orion-Bar-H2-minus-045-Brg.fits")["CHIP1.INT1"]
w = WCS(h2_hdu.header)
fig, ax = plt.subplots(
figsize=(4.5, 5),
subplot_kw=dict(projection=w),
)
m = (h2_hdu.data > 2400.0) | (h2_hdu.data < 750)
h2_hdu.data[m] = 1314.0 # median value
ax.imshow(h2_hdu.data, vmin=1000, vmax=2400, cmap="gray_r")
x0, y0 = 1740, 870
dx, dy = 550, 800
figfile = f"h2-jet-hh204.jpg"
ax.set(
xlabel="RA (J2000)",
ylabel="Dec (J2000)",
xlim=[x0 - dx/2, x0 + dx/2],
ylim=[y0 - dy/2, y0 + dy/2],
)
fig.savefig(figfile, dpi=300)
Doing it in Affinity Designer
- I suspect that HH204 jet is the PA=153 jet that looks like a sheath on the HH203 jet
- The HH203 jet is at PA=117 although could be slightly higher
- This PA is from putting a line through the knot peaks
- But really, the only good PA will come from proper motions
- And proper motions of jet knots at that
- I use the vcube.oiii.fits and vcube.nii.fits and vcube.ha.fits
- Last one is best s/n
- Contours of Ha
- cyan: channel 16, 17, 18 = -38 +/- 6 km/s
- blue: channel 9, 10, 11 = -66 +/- 6 km/s
- Note that velocities are wrt Ha line center, which is about 18 km/s heliocentric
- Use [O III] at -54 km/s for the PA120 jet
- Slightly to the blue of peak for better isolation
- Use [N II] at -38 for the PA150 jet
- Slightly to the red of peak for better isolation
- All columns except V_ONC are heliocentric
[N II] | Ha | [O III] | BEST | V_ONC | |
---|---|---|---|---|---|
PA150 jet | -42 | -40 | -40 | -40 | -65 |
PA120 jet LO | -52 | -50 | -46 | -50 | -75 |
PA120 jet HI | -52 | -48 | -48 | -48 | -73 |
HH 203 head | -50 | -48 | -46 | -48 | -73 |
Splash | -48 | -44 | -45 | -44 | -69 |
HH 204 head | -23 | -18 | -20 | -20 | -45 |
HH 204 pre | -20 | - | -18 | -18 | -43 |
[N II] knot | -20 | -20 | - | -20 | -45 |
HH 204 N edge | -2 | -4 | -2 | -2 | -27 |
HH 204 S edge | -2 | -2 | 2 | 0 | -25 |
HH 203 S head | -14 | -10 | -14 | -14 | -39 |
- This may be related to the HH 204 outflow
- Take Robberto 2005 ACS as ground truth, as always
- Align on the proplyd 197-427
- Celestial coordinates from 2005: 83.8318880 -5.4073798
- Pixel coordinates:
- 2005 ACS: 1190.0 1879.0
- 2015 WFC3: 2097.5 1656.0 2097.5236 1657.0835
- 1996 WFPC2: 799.7 998.9
- Make special WCS files for the 1996 and 2015 images
- Rotation seems off in 1996
- 0.1 / 3600 = 2.77777777778e-5
- Currently CD1_2 / CD1_1 = 3.804223E-08 / -2.76664E-05 = -1.37503361478e-3
- -0.078783573163 deg
- CD2_1 / CD2_2 = 4.913188E-09 / 2.764793E-05 = 1.77705455707e-4
- 0.0101817725013 deg
- CD1_1 = -2.76664E-05
- arccos(2.76664E-05 / 2.77777777778e-5) = arccos(0.996) = 5.133 deg
- That is enormous! Must mean that
- Rotation seems off in 1996
ORDER="hh204-rgb-1993.png hh204-rgb-2005.png hh204-rgb-2015.png"
convert $(for a in "$ORDER"; do printf -- "-delay 100 %s " $a; done; ) HH204-motions.gif
open -a Safari hh204/HH204-motions.gif
files="mean-Fe_III-4658-patfixx.fits
mean-Fe_II-8617-patfixx.fits
mean-Ar_IV-4740-patfixx.fits
mean-Cl_IV-8046-patfixx.fits
mean-H_I-6563-patfixx.fits
mean-H_I-4861-patfixx.fits
mean-He_I-5876-patfixx.fits
mean-He_I-6678-patfixx.fits
mean-N_II-6548-patfixx.fits
mean-N_II-6583-patfixx.fits
mean-O_I-6300-patfixx.fits
mean-O_III-4959-patfixx.fits
mean-O_III-5007-patfixx.fits
mean-S_II-6716-patfixx.fits
mean-S_II-6731-patfixx.fits
mean-S_III-9069-patfixx.fits"
for f in $files;do
cp -pv /Volumes/SSD-2TB/OrionMuse/LineMaps/$f .
done
files="linesum-Fe_III-4881.fits
linesum-Fe_III-4881-bin004.fits
linesum-Fe_III-4658-bin004.fits
linesum-Fe_III-5270-bin004.fits
linesum-Fe_III-4702-bin004.fits
ratio-4702-4658-multibin-SN0010.fits
ratio-4658-4861-bin004.fits
ratio-9069-9229.fits
ratio-6583-6563.fits
linesum-Fe_III-4702-multibin-SN0005.fits"
for f in $files;do
cp -pv /Volumes/SSD-2TB/OrionMuse/LineMaps/$f .
done
Original batch
files="linesum-Fe_III-5270-multibin-SN0005.fits
linesum-Fe_III-4658-multibin-SN0005.fits
linesum-Fe_II-8617-multibin-SN0005.fits
linesum-Ni_II-7378-multibin-SN0005.fits
linesum-Si_III-5740-bin032.fits
linesum-S_III-9069.fits
linesum-S_III-6312.fits
linesum-Ar_III-7136.fits
linesum-O_III-5007.fits
mean-Ar_III-7136-patfixx-bin001.fits
mean-O_II-7318-patfixx-bin001.fits
mean-O_I-8446-bin001.fits
sigma-O_I-8446-bin001.fits
mean-Ni_II-7378-bin008.fits
muse-derived-Ne-iii-multibin-SN0003.fits
muse-derived-Ne-multibin-SN0010.fits
muse-derived-Te-multibin-SN0030.fits
muse-derived-Te-iii-multibin-SN0003.fits
linesum-O_III-5007.fits
linesum-O_II-7318.fits
linesum-O_I-8446.fits
linesum-N_II-6583.fits
linesum-Fe_III-5270-multibin-SN0005.fits
linesum-O_II-7330.fits
linesum-O_I-8446-multibin-SN0005.fits
linesum-Ni_II-7378-multibin-SN0005.fits
linesum-Fe_II-8617-multibin-SN0005.fits
linesum-S_III-9069.fits
linesum-Ar_III-7136.fits
linesum-O_III-5007.fits
ratio-5007-7318.fits"
for f in $files;do
cp -pv /Volumes/SSD-2TB/OrionMuse/LineMaps/$f .
done
import numpy as np
from astropy.io import fits
hdu = {
id_: fits.open(f"linesum-{id_}.fits")[0]
for id_ in ["O_III-5007", "O_II-7330", "O_II-7318", "O_I-8446"]
}
pairs = [
["O_I-8446", "O_II-7330"],
["O_II-7318", "O_II-7330"],
["O_II-7318", "O_III-5007"],
]
for wav1, wav2 in pairs:
ratio = hdu[wav1].data/hdu[wav2].data
fits.PrimaryHDU(
header=hdu[wav1].header,
data=ratio,
).writeto(f"ratio-{wav1}-{wav2}.fits", overwrite=True)
- 4881/4658 varies from 0.3 @ 1000 pcc to 0.45 @ 1e4, then falls again at higher densities
- 5271/4658 is mainly T sensitive when n = 1e3 → 1e5: 0.6 @ 6000 K, falling to 0.4 @ 14,000 K
- 4702/4658 is density-sensitive (and monotonic this time), from 0.26 @ 1000 pcc to 0.33 @ 1e4 pcc to 0.46 @ 1e5 pcc
import numpy as np
from astropy.io import fits
hdu = {
wav: fits.open(f"linesum-Fe_III-{wav}-bin004.fits")["SCALED"]
for wav in ["4881", "5270", "4702", "4658"]
}
pairs = [
["4881", "4658"],
["5270", "4658"],
["4702", "4658"],
]
mask = (hdu["4658"].data > 15000.0)
for wav1, wav2 in pairs:
ratio = hdu[wav1].data/hdu[wav2].data
ratio[~mask] = np.nan
fits.PrimaryHDU(
header=hdu[wav1].header,
data=ratio,
).writeto(f"ratio-Fe_III-{wav1}-{wav2}-bin004.fits", overwrite=True)
Measurements of features in [Fe III] lines
Feature | 4658 | 4702/4658 | 4881/4658 | 5270/4658 | N(S II)/1000 |
---|---|---|---|---|---|
Outside Bar | 0.26 +/- 0.06 | 0.37 +/- 0.05 | 0.56 +/- 0.06 | ||
Bar NE | 0.29 +/- 0.04 | 0.42 +/- 0.04 | 0.57 +/- 0.05 | ||
HH 204 | 0.30 +/- 0.05 | 0.45 +/- 0.03 | 0.64 +/- 0.06 | ||
HH 203 | 0.31 +/- 0.07 | 0.43 +/- 0.03 | 0.61 +/- 0.04 | ||
Bar SW | 0.31 +/- 0.06 | 0.44 +/- 0.06 | 0.58 +/- 0.06 | ||
HH 529-III/II | 0.31 +/- 0.07 | 0.46 +/- 0.07 | 0.59 +/- 0.08 | ||
knots | 0.32 +/- 0.07 | 0.48 +/- 0.06 | 0.56 +/- 0.13 | ||
shock | 0.32 +/- 0.05 | 0.47 +/- 0.05 | 0.66 +/- 0.06 | ||
minibar | 0.32 +/- 0.04 | 0.43 +/- 0.02 | 0.60 +/- 0.06 | 8.4 +/- 1.0 | |
HH 202 | 0.33 +/- 0.04 | 0.48 +/- 0.05 | 0.59 +/- 0.04 | 5.0 +/- 1.2 | |
HH 529-0 | 0.35 +/- 0.03 | 0.45 +/= 0.06 | 0.68 +/- 0.05 | ||
HH 269 base | 0.37 +/- 0.06 | 0.32 +/- 0.02 | 0.64 +/- 0.06 |
- Slit size: 1 x 12 arcsec in red; 1 x 10 arcsec in blue
- Slit center: RA(J2000) 05:35:16.80, DEC(J2000) -05:23:57.48
- For instance ~/Dropbox/SPMFEB13/OrionS/spm-vstack-oiii.fits
- Zoom factor needs to be 4.5x10 to match 1.68 on MUSE images
- Bowshock III and bowshock-east show up best in the Vhel = +4 channel
- actually, +5.6 - see below
- II and the III knot show up best around -28
- 0 and Ram show peak atound -46
- strange velocity gradient
- Diagonally across II/III from NW → SE
- III bowshock contrast (wrt interpolation from either side)
- -4.4 : 1500 = 300% - starting to see knot too
- -2.4 : 2000 = 400%
- -0.4 : 2400 = 400%
- +1.6 : 2500 = 250%
- +3.6 : 2700 = 180%
- +5.6 : 2800 = 107%
- +7.6 : 2700 = 67%
- +9.6 : 2500 = 42%
- +11.6 : 2000 = 25%
- +13.6 : 1750 = 20%
- +15.6: 1500 = 16%
- +17.6: 1000 = 11%
- So, peak is at +5.6
- Also, half-maximum points are at [-4.4, +15.6], with mean at +5.6 too
- This is strongest in S wing of III
- There are some bright knots just below HH 529-III
- ~/Dropbox/Orion-HH-data/MUSE/fe-ii-knots.reg
- They turn out to be HH 1149-CW bows
- I had thought this might be feeding HH 203
- But now it doesn’t look like it
- Some interesting low-ionization knots near HST 10
- ~/Dropbox/Orion-HH-data/MUSE/oi-sii-knots.reg
- Possibly an extension of the HH 1149-CW flow
- Not moving very fast
- At least that is what I get by comparing by eye F658N 1996 with F373N 2015
final658-radec.fits
- Shifted 5 pix to CRPIX1=1360
icaz02040_drz.fits
- Also possibly visible in 2005 ACS F658N
- At least that is what I get by comparing by eye F658N 1996 with F373N 2015
- Bow shock with fat shell driven by p159-350 (?) moving to NE at 15 km/s
- Fainter, thinner bow shock moving to E at 25 km/s
- Fed by jet from p155-337
- These look like these are a super-slow jet that is feeding HH 1149-CW
- But actually they are the wings of the HH 529 bow shocks passing over a stationary filament
- Star is AC Ori, also known as:
- [HC2000] 202 - near-IR
- COUP 768 - X-ray
- JW503
- [FRM2016] 266 - Radio
- This is projected onto the HH 529-III shock
- It has a shadow disk, surrounded by faint, circularly symmetric emission
- This is similar to 171-340, but even more circular
- Looking in catalogues
- Coordinates are 5:35:16.9570 -5:23:59.567
- Searching in SIMBAD finds MLLA 264
- Muench:2002a - The Luminosity and Mass Function of the Trapezium Cluster: From B Stars to the Deuterium-burning Limit
- Source 264 in catalog
- Coordinates 05 35 16.95, -05 23 59.6
- Agree within 0.1 arcsec
- Full record for this source
- Only photometry is from NTT
- J = 15.26 +/- 0.16
- H = 14.62 +/- 0.06
- K = 13.90 +/- 0.01
- Cross-ref to [HC2000] 182
- This is Hillenbrand:2000a
- Same coords (within 0.1 arcsec)
- H=14.264, K=13.237 mags are a bit brighter
- Other matches
- [FRM2016] 329
- Forbrich:2016a - Radio 4cm and 6cm survey
- Flux density 0.1964 +/- 0.0026 mJy
- Eisner:2018a
- Calls it HC182
- F(850 mic) = 2.0 +/- 0.2 mJy
- F(dust) = 2.0 +/- 0.2 mJy
- M(dust) = 3.5 +/- 0.3 Mearth = 1e-5 Msun
- R(disk) = 14.8 +/- 1.9 AU = 0.03 arcsec
- This puts it at the faint and small end of the distribution of detected sources
- In the HST images I measure a disk shadow diameter of 0.13 arcsec
- This gives a radius of 0.065 arcsec = 27 AU
- [FRM2016] 329
- 5:35:17.5401 -5:23:55.023
- Previously classified as proplyd in OW94
- Very similar morphology
- Circular rim of Ha
- Inner disk extinction
- Both excitation mechanisms seem to be operating for these lines
- 5158.7924 Å is roughly equal parts collisional and fluorescent in the Keck spectrum of 244-440 adjacent nebula
- We clearly see the velocity splitting of about 10 km/s between the two components
- 5262 is mainly fluorescent, but weaker
- Both lines are fluorescent in the proplyd itself - maybe collisional deexcitation of collisional line
- Can also see differences in the MUSE maps, which cover broader range of wavelengths
- 7155 is mainly collisional
- 4815 is mainly fluorescent
- This line is also in the Keck data, but I seem to have never reduced it
- Actually I did
- It is slightly weaker than 5262 and has another high-ionization line next to it at 4816
- This line is also in the Keck data, but I seem to have never reduced it
- 8617 is the strongest and is a combination of both
- 5158.7924 Å is roughly equal parts collisional and fluorescent in the Keck spectrum of 244-440 adjacent nebula
- There is plenty of discussion in Keck notes
- [[id:4F674E2D-0B0C-4584-9333-40C6131BD3F2][[Fe II] lines]]
- These are also combinations of at least two mechanisms
- Weilbacher:2015a - Fig 28 shows maps of velocity difference nii-ha and sii-nii
- But these are mean velocities, so they are heavily effected by high-velocity components, so pretty useless for studying the main ionization front except for in restricted regions of the nebula where we can be sure there are no jets, etc
- Garcia-Diaz:2008a - Fig 18 shows distribution of peak velocities for different pairs of lines
- This should be relatively unaffected by high-velocity flows except in the restricted regions where they dominate the brightness
- Table 2 shows the average differences
- nii - oiii = 5.6 km/s
- sii - siii = 2.6 km/s
- oi - sii = 4.5 km/s
- Henney:2005b appendix A gives a simple model of plane parallel weak-D front
- This includes neutral and singly-ionized line only
- Does not include a neutral fluorescent line, but that should be represented approximately by zero velocity
[2020-10-12 Mon]
This paper provides an abundance analysis for two photoionized HH objects in the Orion nebula. A main conclusion from the work is that iron appears depleted, but some is released into the gas phase by the HH shock. There is also quite a bit of discussion about fluctuations, and abundance discrepancy factors between collisionally-excited lines and recombination lines, as well as some proper motion and velocity measurements. My overall reactions reading through this were (a) wow, this is a lot of work, it takes effort to sift through hundreds of line ratios, and (b) the data look lovely and the analysis seems very thorough - meticulous error bars on everything, nuances of each line understood, even tidbits associated with 'ghost' lines. The writing style is clear. It's an impressive work that certainly deserves publication in MNRAS. I do have a couple of comments, and then a few minor ones. (*) One consequence of being so thorough is that the paper does drag, and is rather long. I'll leave it to the authors to consider, but I saw little use in all the details of Table 13 and 14 (or even Tables 10-12), or the long discussions of each line ratio in section 8. I'm not sure how useful the whole t^2 formalism is for a shocked filamentary object. I don't know what we learn from that, really. Some of these tables and discussion might live more happily in an Appendix. It's up to the authors, but they should take a step back from this tome of data and be sure that the primary takeaways stand out from the details, taking care to only include tables that are relevant to the main messages. (*) An important assumption that underlies this entire analysis is that a single temperature and density, with some fluctuations, characterizes the emitting gas. This assumption is only true if the cooling zone behind the shock does not emit enough to alter the spectrum, so the shock merely serves to create a dense blob that can be photoionized by the ambient radiation. However, this point is not really made until the first paragraph in section 11 on page 24, with a quick reference to Henney 2002. There is then a reference to compression ~ M^2, but this is for the non-magnetic case, and even a small seed field changes that in HH shocks (see 2015 ApJ 811, 12). I feel that this issue warrants discussion much earlier. I would mention it in the introduction at the level it is mentioned in section 11, and then make a new subsection that would go in front of section 4.1, and be a bit more pedagogical. A radiative shock wave is not amenable to standard Te/Ne analysis because there is no single Te, Ne. That's why so little information exists on abundances in HH objects. As radiation increases the first effect is to ionize the preshock gas, which can increase flux in some lines (e.g. HH 47D has strong [O II] because zeta Pup and gamma Vel ionize O to O+, so a weak shock shows [O II]). The next level up in ionization is that ionization enhances forbidden lines by kicking out electrons in the cooling zone. Then as ionization increases, the sheaths of jets get photoionized, like in Carina, and finally you photoionize the whole thing so that nothing goes below about 8000K. After that there is some point where you just don't see the cooling zone at all. The paper should argue based on some rough numbers that the radiation field is high enough in Orion that you are in the last category. Are we really sure that the enhanced OIII abundance in HH 529 doesn't arise in some way from the shocked cooling zone? At some point if you move HH 529 away from the photoionizing sources it's going to look more like a preionized shock with a cooling zone, like HH 47D. Note, a shock can be 'fully ionized' in the sense that it has no OI, yet still be affected by the cooling zone behind the shock in that it will enhance OIII relative to OII. I have to wonder if this also affects the ADF in HH 529. Some minor comments Fig 1 caption: The labeling is unclear. It took a while for me to realize that I, II, and III were roman numerals and not |, ||, ||| indicative of slit positions or something. The authors need to move the labels a bit away from the objects and then mark them with arrows, e.g. II ----> (points to object) p 2: add a reference to Fig 1 inside the parentheses for the O'Dell/Henney reference: \citep[][; Fig~1]{odellhenney08} so the reader knows the objects are labeled in the figure. p13: Table 10 should note these are logarithmic abundances with n(H) defined to be 12. p25: The compression really gives the magnetosonic Mach number. Some inclusion of magnetic field discussion is needed here. Finally, some discussion of Fe destruction in HH shocks is needed. A paper just out sees evidence for it in the Mach disk component of HH 32 (2020 AJ 160, 165), and references in that paper will guide to some other attempts to measure this.
- proper-motions/hh529-finding-chart.pdf
- Reply to referee
We have made the suggested changes Fig 1 in order to make the labelling clearer. In addition, we have changed the font to make it more obvious that I, II, and III are roman numerals.
- Energy flux through shock is 0.5 m n_0 V^3 = 0.5 m n_0 M^3 c_0^3
- Units erg/cm^2/s
- Equilibrium cooling from shell is Λ n_2^2 h
- Same units:
- n_2 = (V/c_0)^2 n_0 = M^2 n_0
- Λ = 2e-24 erg cm^3 / s at 1e4 K
- Ratio is
- 0.5 m n_0 M^3 c_0^3 / Λ M^4 n_0^2 h
- = (0.5 m c_0^3 / Λ) 1/ M n_0 h
- = 5e17 cm^-2 / M n_0 h
- h = R / M^2
- So Ratio = 5e17 M / R n_0
- n_0 = 100 pcc
- R = 10 arcsec = 6e17 cm
- M = 8 or so
- => Ratio = 0.067
- So not that small!
- Star radial velocity is 5.74 LSR (Boyden+ 2020), so 23.84 km/s heliocentric
- Very close to the average cluster velocity of 25
- Claimed errors are very small: less than 0.1 km/s
- Knot and jet velocity are +150 km/s heliocentric
- So relative velocity of jet wrt star is 150 - 24 = +126 km/s
- Theissen 2022ApJ…926..141T
- Observes our star according to Table 1 (HC252)
- But not included in any following tables
- [2022-04-01 Fri 09:26] I sent an email to the corresponding author
- Source IDs
- [HC2000] 252
- [OW94] 170-337
- 2019AJ....157..109K
- 170-337 is source 118
- Vizier record of source 118
- Proper motion
- RA: 2.97 +/- 0.37 mas/yr
- Dec: 2.67 +/- 0.68 mas/yr
- Total: 3.994 +/- 0.5314
- If distance is 410 pc, then parallax = 2.439 mas
- V = 4.74 mu / P
- V = [5.772 +/- 0.7191, 5.189 +/- 1.322] km/s
- Total (8 +/- 1) km/s at PA = (48 +/- 8) deg
- This out to be the ideal reference
- HST ASTROMETRY IN THE ORION NEBULA CLUSTER: CENSUS OF LOW-MASS RUNAWAYS
- But 170-337 is not included in their catalog!
- O’Dell & Henney 2008
- 48 km/s at PA 357, but assuming D = 440 pc
- Correcting to 410 pc
- RA = (410/440) 48 sin(357) = -2.341 km/s
- Dec = (410/440) 48 cos(357) = 44.67 km/s
- Subtract off motion of star
- RA = -2.341 - 5.772 = -8.113
- Dec = 44.67 - 5.189 = 39.481
- V = hypot(-8.113, 39.481) = 40.31 km/s
- PA = arctan(-8.113 / 39.481) = 348 deg
- V_r = (+126 +/- 10) km/s
- V_t = (40 +/- 6) km/s
- V3d = 132.2 +/- 9.7
- Inclination to POS is 72.39 +/- 2.804
- So inclination to LOS is 18 +/- 3 deg
- As opposed to 17 deg last time
- Disk mass
- Dust: 23 +/- 5 Mearth = 6.9e-5 Msun
- Gas:
- assume gas/dust = 100
- 0.007 Msun
- I load up the Muse color images of different ionization zones in Preview, as different tabs of the same window
- I match the image sizes, for instance with
cmd 0
- I
ctl tab
to switch between the tabs - Optionally I use
`
to activate the magnification window - Everything that is an ionization front will show movement
- the less dense it is, the further it will move
- whereas things that stay still may be various things:
- extinction features if they happen to be dark
- or they could be so dense that we can’t see the ionization gradient
- or they could be jets, where the ionization is collisional, rather than radiative
- using this technique, I have identified the jet driving HH 203/4 I think
- actually, not so fast - not really clear whether it is or not
- it is definitely a jet, but it might be an unrelated flow
- it is pretty low ionization
- there is also the faint [O I] filament that we mentioned in the Tsamis 2013 paper
- and there is the evidence from the high-velocity [O III] and where it disappears
- look into this more when I have more time