UVES observations of HH objects in Orion

Original Proposal

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 2019 August

Discussions between Eduardo and Will during Eduardo’s visit to Morelia

The back-scattered component

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

Finding chart for nebular features

This is from Garcia-Diaz:2007a

Other ideas for Eduardo to pursue

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

Background reading

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”
  - http://adsabs.harvard.edu/abs/2007AJ....133.2192H
- Garcia-Diaz:2007a “Velocity Structure in the Orion Nebula. I. Spectral Mapping in Low-Ionization Lines”
  - http://adsabs.harvard.edu/abs/2007AJ....133..952G
- Garcia-Diaz:2008a “Velocity Structure in the Orion Nebula. II. Emission Line Atlas of Partially Ionized to Fully Ionized Gas”
  - http://adsabs.harvard.edu/abs/2008RMxAA..44..181G
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”
  - http://adsabs.harvard.edu/abs/2015AJ....150..108O
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”
  - http://adsabs.harvard.edu/abs/2009MNRAS.395..855M
- Espiritu:2017a “DUST DESTRUCTION AND THERMAL INHOMOGENEITIES IN THE ORION NEBULA. RESULTS OF HH 202.”
  - http://adsabs.harvard.edu/abs/2017cfe..confE..10E
- Weilbacher:2015a “A MUSE map of the central Orion Nebula (M 42)”
  - http://adsabs.harvard.edu/abs/2015A%26A…582A.114W
- McLeod:2016a “A nebular analysis of the central Orion nebula with MUSE”
  - http://adsabs.harvard.edu/abs/2016MNRAS.455.4057M
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
Neutral gas
- van-der-Werf:2013a
Atomic physics
- Laha:2017a “Emission Line Ratios of FE III as Astrophysical Plasma Diagnostics”
  - http://adsabs.harvard.edu/abs/2017ApJ…841....3L

Large-scale velocity maps

Combined channel maps of [O III] and [N II]

figs/multi-panel-isovel-oiii-nii.pdf

Velocity-colorized [O III] 5007 image

Velocity-colorized [N II] 6583 image

Shock emission in HH 204

Final stuff on HH 204

Abstract should not exceed 250 words

Answer Eduardo’s earlier queries

Muy buenas preguntas las dos!

Anchored Note, page 3

¿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

Anchored Note, page 3

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?

Cygnus loop analogy

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

Density jump in the jet

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

Estimate density from Hα brightness

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

Explanation for the profiles along the slit

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

(1)	(2)	(3)	(4)	(5)	(6)
pos	T, K	n, pcc	e4959/4861	I4959/4861	12+log(O++/H+)
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

Trapping of the ionization front in the jet

We see an ionization gradient in the blue-shifted component

Oxygen ionization ratios

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

Notebook in Google Colab

file:~/Google Drive/Colab Notebooks/HH204 neutral oxygen fraction.ipynb

Extinction correction from C(Hb)

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] versus Hα emission

[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?

Return to this issue for the HH 204 paper

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

pyds9 commands to open the files

This is based on what I did years ago in
- ~/Dropbox/Org/hh-orion.org
The PC images are described in
- ~/Work/BobPC/2005/README
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

Plan to get the [O III] rim analysis done

Align the images to Robberto 2005

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

Correct for continuum using F547M

This is now done in Finally make ratios 502/656 and 658/656 using the calibration constants from ODell:2009b

Make ratio images of [O III] / Ha

First make them 2D images - currently they are 3d

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

Reproject all to the F502 corrected WCS

[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

Comparison with UVES values

First we subtract continuum, then we subtract upstream BG

[O III] 4363

Region	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

[O III] 4959

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

Hb 4861

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

Red component (nebula)

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

Blue component

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

Comparison between red (R) and blue (B) components

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

Implications for the rim increase in the [O III] / Ha ratio

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)

Remove upstream nebula from the 2D spectra

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

Do plots of profiles through the shock

Cancelled for being too much detail

Look at the 2015 WFC3 images

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.

Do proper motions for HH 204

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.

Reproject HH 204 images onto the 2005 coordinate grid

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]]

Extract ACS image around HH 204

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

Put them all on the same grid

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

Make ratio images

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

Do flct to get 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("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)

Combine the two PM intervals to get variance

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

Make a figure of HH 204 proper motions

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)

Measure per-feature motions for HH 204

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
  - hh204/hh204-1996-nose-frame.wcs
  - hh204/hh204-2015-nose-frame.wcs
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

Align in frame of reference of the bow shock nose

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

Make high-pass filtered versions

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

Make GIFs

Save images of the FITS files

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)

GIFs in bow shock nose frame

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

Redo proper motions in bow shock frame

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

Redo the combining of the two intervals

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

Make a figure of the proper motions in the bow shock frame

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)

Conclusions from proper motions and new directions

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

Document with Will’s additions to paper on HH 204

will-hh204-figs.tex

Looking for evidence of the jet in HH 204

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
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
- Although none of these PA seem compatible with the observed proper motion of HH 204
  - Or are they?

Results from Adal’s previous UVES observations of HH202-S

Dust destruction

Fe, Ni, Cr abundance higher by factor of about 10 in shocked gas

ADFs

Gets high ADF for C²⁺ 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

First look at data

HH 529 dataset

Eduardo talk [2021-06-17 Thu]

Will presentation of Eduardo

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

Eduardo’s text

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.

Bob comments [2021-07-13 Tue]

HH 204 jets

Original message

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

Proper motions

HH 529

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

Relevant programs

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

Plan of action

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

WFPC2 1996 → ACS 2005 → WFC3 2015 F656N

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

Extract just a window of the ACS 2005 Ha+[N II] image

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

Reproject the 1996 and 2015 epochs onto 2005

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

Use FLCT to calculate proper motions 1996 → 2005 → 2015

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)

High-pass filtered version of the 2015 epoch WFC3 Hα filter

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

Make vector graphs of the proper motions.

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

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

Average and dispersion 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)

Make some ratio images

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

Measure the per-feature motions in HH 529-III and II

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

Which proper motion features are in the UVIS slit

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

Comparison with radial velocities

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!

ACS 10-year interval

Fix the WCS of epoch 2

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

Make some GIFs

Frame of reference of star cluster

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

Frame of reference of HH 529-III bow shock

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

Measurements of HH 529 PM

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

Align on the HH529-III motion

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

Figure showing the ACS proper motions for HH 529

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

Radio observations

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

Permitted lines

This is now in a separate file because it was too big:

muse-oii-permitted.org

Make a new version of Fig 1

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

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

Make lupton versions

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

Repeat for HH 204 - Fig 1 of Paper II

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

Aligning the epochs for HH 204

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

Re-do the rgb image for HH 204

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

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

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

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)

Make RGB channel maps of nii and ha

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

Image of the H_2 emission from the possible HH 203 jet

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)

Now make the large-scale map that shows the high-velocity flow

Doing it in Affinity Designer

Comments on the axes

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

Add contours of the jets

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

New strategy

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

Measurements of Radial velocities of different features

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

The low-ionization knot above the Bar

This may be related to the HH 204 outflow

Aligning the HST image

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

GIFs of HH 204 1993, 2005, 2015

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

Compare with MUSE spectral maps

Copy the files that we want from Volumes/SSD-2TB/OrionMuse/LineMaps

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

Extra line ratios

O I / [O II]

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)

[Fe III] density and temperature sensitive

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

Region for the slit

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

Compare with SPM horizontal slits

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

Analysis of Bowshock III profile

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

Yet another component at -12 km/s

This is strongest in S wing of III

Implication for 3D kinematics

Flow from HH 1149-CW down to Bright Bar

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

Interesting proper motion projects in the field of HH 529

Colliding shocks - 162-341

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

Overlapping filaments - 152-359 and 153-360

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

Moving dark filament

Multipolar jet from 160-353

Star is AC Ori, also known as:
- [HC2000] 202 - near-IR
- COUP 768 - X-ray
- JW503
- [FRM2016] 266 - Radio

New proplyd? 170-360

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

Comparison with 175-355

5:35:17.5401 -5:23:55.023
Previously classified as proplyd in OW94
Very similar morphology
- Circular rim of Ha
- Inner disk extinction

Acceleration in the main ionization front

Data from the Keck proplyd study

[Fe II] lines - collisions and fluorescence

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
  - 8617 is the strongest and is a combination of both
There is plenty of discussion in Keck notes
- [[id:4F674E2D-0B0C-4584-9333-40C6131BD3F2][[Fe II] lines]]

Si II lines - recombination and fluorescence

These are also combinations of at least two mechanisms

Previous empirical studies

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

Previous theory

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

Dealing with the referee report

[2020-10-12 Mon]

Referee report

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.

New version of Figure 1

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.

Shock versus shell emission

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!

HH 514 jet from proplyd 170-337

Figure 1

figs/hh514-finding.pdf

Radial velocity of jet and star

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

Proper motions

Proper motions of star

Source IDs
- [HC2000] 252
- [OW94] 170-337

Kim et al 2019

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

Platais et al 2020

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!

Proper motion of jet knot

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

3D motion of jet

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

Masses and mass loss rates

Disk mass
- Dust: 23 +/- 5 Mearth = 6.9e-5 Msun
- Gas:
  - assume gas/dust = 100
  - 0.007 Msun

A new way of identifying HH jets

Basic method

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

Name		Name	Last commit message	Last commit date
Latest commit History 61 Commits
__pycache__		__pycache__
auto		auto
docs		docs
eduardo-2019-09		eduardo-2019-09
eduardo-2020-06		eduardo-2020-06
figs		figs
hh204		hh204
hh204_multi_epoch_nose_align_EXPORT		hh204_multi_epoch_nose_align_EXPORT
line-shifts		line-shifts
muse		muse
notebooks		notebooks
press-release-hh204		press-release-hh204
proper-motions		proper-motions
shock-vs-shell		shock-vs-shell
.gitignore		.gitignore
.log		.log
HH204-bowshock-frame.gif		HH204-bowshock-frame.gif
README.org		README.org
h2-jet-hh204.jpg		h2-jet-hh204.jpg
hh204-N-wing-oiii-ha-ratio.pdf		hh204-N-wing-oiii-ha-ratio.pdf
hh204-finding-chart-lores.pdf		hh204-finding-chart-lores.pdf
hh204-finding-chart-simple.afdesign		hh204-finding-chart-simple.afdesign
hh204-finding-chart-simple.pdf		hh204-finding-chart-simple.pdf
hh204-finding-chart.afdesign		hh204-finding-chart.afdesign
hh204-finding-chart.jpg		hh204-finding-chart.jpg
hh204-finding-chart.pdf		hh204-finding-chart.pdf
hh204-medium-slit-oiii-ha-ratio.pdf		hh204-medium-slit-oiii-ha-ratio.pdf
hh204-ratio-oiii-ha-HST-PC-1994.pdf		hh204-ratio-oiii-ha-HST-PC-1994.pdf
hh204-ratio-oiii-ha-annotated.afdesign		hh204-ratio-oiii-ha-annotated.afdesign
hh204-ratio-oiii-ha-annotated.pdf		hh204-ratio-oiii-ha-annotated.pdf
hh204-slit.reg		hh204-slit.reg
hh204-uves-slit-oiii-ha-ratio.pdf		hh204-uves-slit-oiii-ha-ratio.pdf
hh204-wide-slit-oiii-ha-ratio.pdf		hh204-wide-slit-oiii-ha-ratio.pdf
hh204_1996_nose_align.jpg		hh204_1996_nose_align.jpg
hh204_2005_nose_align.jpg		hh204_2005_nose_align.jpg
hh204_2015_nose_align.jpg		hh204_2015_nose_align.jpg
hh204_f656n_1996_2005_make_pm_map.out		hh204_f656n_1996_2005_make_pm_map.out
hh204_f656n_1996_2015_2stage_proper_motions.jpg		hh204_f656n_1996_2015_2stage_proper_motions.jpg
hh204_multi_epoch_nose_align.afdesign		hh204_multi_epoch_nose_align.afdesign
hh204_nose_align_f656n_1996_2005_make_pm_map.out		hh204_nose_align_f656n_1996_2005_make_pm_map.out
hh204_nose_align_f656n_1996_2015_2stage_proper_motions.jpg		hh204_nose_align_f656n_1996_2015_2stage_proper_motions.jpg
hh204_nose_align_f656n_1996_2015_2stage_proper_motions.pdf		hh204_nose_align_f656n_1996_2015_2stage_proper_motions.pdf
hh204_nose_align_f656n_2005_2015_make_pm_map.out		hh204_nose_align_f656n_2005_2015_make_pm_map.out
hh529-iii-slit.reg		hh529-iii-slit.reg
hh529-knots.tex		hh529-knots.tex
make_pm_map.py		make_pm_map.py
mosaic-1996-HH204-align-robberto.wcs		mosaic-1996-HH204-align-robberto.wcs
muse-oii-permitted.org		muse-oii-permitted.org
oii-fluorescence.org		oii-fluorescence.org
oiii-4959-BG-spec.dat		oiii-4959-BG-spec.dat
read_wcs.py		read_wcs.py
rgb-lupton-channels-oiii.jpg		rgb-lupton-channels-oiii.jpg
rgb-lupton-channels.jpg		rgb-lupton-channels.jpg
rgb-lupton-hh204-1996-contours.jpg		rgb-lupton-hh204-1996-contours.jpg
rgb-lupton-hh204-1996-simple.jpg		rgb-lupton-hh204-1996-simple.jpg
rgb-lupton-hh204-1996.jpg		rgb-lupton-hh204-1996.jpg
rgb-lupton-hh204-1996fix.jpg		rgb-lupton-hh204-1996fix.jpg
rgb-lupton-hh529-1996.jpg		rgb-lupton-hh529-1996.jpg
rgb-lupton-hh529-1996b.jpg		rgb-lupton-hh529-1996b.jpg
rgb-lupton-hh529-1996b.pdf		rgb-lupton-hh529-1996b.pdf
rgb-lupton-hh529-2005.jpg		rgb-lupton-hh529-2005.jpg
rgb-lupton-hh529-hh204.jpg		rgb-lupton-hh529-hh204.jpg
rgb-lupton-hh529-iii-1996.jpg		rgb-lupton-hh529-iii-1996.jpg
rgb-lupton-hh529.py		rgb-lupton-hh529.py
testin.dat		testin.dat
testout.dat		testout.dat
two-temp-TB-30000.pdf		two-temp-TB-30000.pdf
two-temp-Toiii-Tsiii.pdf		two-temp-Toiii-Tsiii.pdf
will-hh204-figs (1).tex		will-hh204-figs (1).tex
will-hh204-figs.log		will-hh204-figs.log
will-hh204-figs.out		will-hh204-figs.out
will-hh204-figs.pdf		will-hh204-figs.pdf
will-hh204-figs.tex		will-hh204-figs.tex
will-hh204-figs_comments.pdf		will-hh204-figs_comments.pdf
will-hh204-figs_comments.txt		will-hh204-figs_comments.txt
will-hh204-refs.bib		will-hh204-refs.bib
will-hh204-two-zone.pdf		will-hh204-two-zone.pdf
working-surface-hh204-single.pdf		working-surface-hh204-single.pdf
working-surface-hh204.jpg		working-surface-hh204.jpg
working-surface-hh204.pdf		working-surface-hh204.pdf

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

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

F469N	0.51e-10
F487N	0.70e-10
F502N	0.82e-10
F656N	1.62e-10
F658N	1.60e-10
F673N	1.42e-10

	Vt 08	Vro 08
HH 529-I	54	-31
HH 529-II	92	-31
HH 529-III	40	-30

div-B-equals-0/orion-hh

Folders and files

Latest commit

History

Repository files navigation

UVES observations of HH objects in Orion

Original Proposal

Discussions 2019 August

The back-scattered component

Finding chart for nebular features

Other ideas for Eduardo to pursue

Background reading

Large-scale velocity maps

Combined channel maps of [O III] and [N II]

Velocity-colorized [O III] 5007 image

Velocity-colorized [N II] 6583 image

Shock emission in HH 204

Final stuff on HH 204

Abstract should not exceed 250 words

Answer Eduardo’s earlier queries

Anchored Note, page 3

Anchored Note, page 3

Cygnus loop analogy

Density jump in the jet

Estimate density from Hα brightness

Explanation for the profiles along the slit

Trapping of the ionization front in the jet

Oxygen ionization ratios

Notebook in Google Colab

Extinction correction from C(Hb)

[O III] versus Hα emission

Return to this issue for the HH 204 paper

pyds9 commands to open the files

Plan to get the [O III] rim analysis done

Align the images to Robberto 2005

Correct for continuum using F547M

Make ratio images of [O III] / Ha

First make them 2D images - currently they are 3d

Reproject all to the F502 corrected WCS

Comparison with UVES values

Remove upstream nebula from the 2D spectra

Do plots of profiles through the shock

Look at the 2015 WFC3 images

Do proper motions for HH 204

Reproject HH 204 images onto the 2005 coordinate grid

Extract ACS image around HH 204

Put them all on the same grid

Make ratio images

Do flct to get proper motions

Combine the two PM intervals to get variance

Make a figure of HH 204 proper motions

Measure per-feature motions for HH 204

Align in frame of reference of the bow shock nose

Make high-pass filtered versions

Make GIFs

Save images of the FITS files

GIFs in bow shock nose frame

Redo proper motions in bow shock frame

Redo the combining of the two intervals

Make a figure of the proper motions in the bow shock frame

Conclusions from proper motions and new directions

Document with Will’s additions to paper on HH 204

Looking for evidence of the jet in HH 204

Results from Adal’s previous UVES observations of HH202-S

Dust destruction

ADFs

First look at data

HH 529 dataset

Eduardo talk [2021-06-17 Thu]

Will presentation of Eduardo

Eduardo’s text

Bob comments [2021-07-13 Tue]

HH 204 jets

Original message

Proper motions

HH 529

Relevant programs

Plan of action

WFPC2 1996 → ACS 2005 → WFC3 2015 F656N

Extract just a window of the ACS 2005 Ha+[N II] image