Skip to content

Commit

Permalink
Add redshift correction to surface brightness calculation. Close #73
Browse files Browse the repository at this point in the history
  • Loading branch information
@kmhess
kmhess committed May 25, 2023
1 parent 864f25b commit 7dd13b9
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 15 deletions.
10 changes: 5 additions & 5 deletions src/make_images.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ def make_overlay_usr(source, src_basename, cube_params, patch, opt, base_contour
return

nhi, nhi_label, nhi_labels = sbr2nhi(base_contour, hdulist_hi[0].header['bunit'], cube_params['bmaj'].value,
cube_params['bmin'].value)
cube_params['bmin'].value, source)

try:
hiwcs, cubew = get_wcs_info(src_basename + '_{}_cube.fits'.format(source['id']))
Expand Down Expand Up @@ -167,7 +167,7 @@ def make_overlay(source, src_basename, cube_params, patch, opt, base_contour, su
return

nhi, nhi_label, nhi_labels = sbr2nhi(base_contour, hdulist_hi[0].header['bunit'], cube_params['bmaj'].value,
cube_params['bmin'].value)
cube_params['bmin'].value, source)
try:
hiwcs, cubew = get_wcs_info(src_basename + '_{}_cube.fits'.format(source['id']))
except FileNotFoundError:
Expand Down Expand Up @@ -267,7 +267,7 @@ def make_mom0(source, src_basename, cube_params, patch, opt_head, base_contour,
owcs = WCS(opt_head)

nhi, nhi_label, nhi_labels = sbr2nhi(base_contour, hdulist_hi[0].header['bunit'], cube_params['bmaj'].value,
cube_params['bmin'].value)
cube_params['bmin'].value, source)
fig = plt.figure(figsize=(8, 8))
ax1 = fig.add_subplot(111, projection=owcs)
plot_labels(source, ax1, cube_params['default_beam'], x_color='white')
Expand Down Expand Up @@ -353,7 +353,7 @@ def make_snr(source, src_basename, cube_params, patch, opt_head, base_contour, s
owcs = WCS(opt_head)

nhi, nhi_label, nhi_labels = sbr2nhi(base_contour, hdulist_hi[0].header['bunit'], cube_params['bmaj'].value,
cube_params['bmin'].value)
cube_params['bmin'].value, source)
wa_cmap = colors.ListedColormap(['w', 'royalblue', 'limegreen', 'yellow', 'orange', 'r'])
boundaries = [0, 1, 2, 3, 4, 5, 6]
norm = colors.BoundaryNorm(boundaries, wa_cmap.N, clip=True)
Expand Down Expand Up @@ -621,7 +621,7 @@ def make_color_im(source, src_basename, cube_params, patch, color_im, opt_head,
mom0 = hdulist_hi[0].data

nhi, nhi_label, nhi_labels = sbr2nhi(base_contour, hdulist_hi[0].header['bunit'], cube_params['bmaj'].value,
cube_params['bmin'].value)
cube_params['bmin'].value, source)

owcs = WCS(opt_head)
fig = plt.figure(figsize=(8, 8))
Expand Down
50 changes: 40 additions & 10 deletions src/modules/functions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,11 @@
import numpy as np
from pvextractor import extract_pv_slice, PathFromCenter

HI_restfreq = 1420405751.77 * u.Hz


###################################################################


def chan2freq(channels, fits_name):
"""Convert channels to frequencies.
Expand Down Expand Up @@ -67,34 +72,59 @@ def felo2vel(channels, fits_name):
return velocities


def sbr2nhi(sbr, bunit, bmaj, bmin):
"""Get the HI column density from sbr.
def sbr2nhi(sbr, bunit, bmaj, bmin, source):
"""Get the HI column density from sbr. See Section 15 of Meyer et al (2017) for equations:
https://ui.adsabs.harvard.edu/abs/2017PASA...34...52M/abstract
:param sbr: SBR
:type sbr: float
:param bunit: unit in which sbr is measured
:type bunit: str
:param bmaj: major axis of the beam
:type bmaj: float
:param bmin: minor axis of the bea,
:param bmin: minor axis of the beam
:type bmin: float
:param source: source object
:type source: Astropy table
:return: column density
:rtype: float
"""
# NEED TO ADD UNITS THAT ANNOYINGLY COME OUT OF SPECTRAL CUBE! DONE?

c = const.c.to(u.m/u.s).value

if (bunit == 'Jy/beam*m/s') or (bunit == 'Jy/beam*M/S'):
nhi = 1.104e+21 * sbr / bmaj / bmin
print("\tWARNING: Assumes velocity axis of cube is in the *observed* frame. If cube is in source rest frame, "
"the column density is (1+z) times greater than shown.")
if ('v_rad' in source.colnames): # or (cube_params['spec_axis'] == 'VRAD'): # Taken from make_images.py.
# First convert to observed frequency, then to redshift following these equations:
# https://web-archives.iram.fr/ARN/may95/node4.html
print("\tWARNING: HI column density calculation assumes optical velocity convention. Data is in radio convention!")
vel_sys = source['v_rad']
freq_sys = HI_restfreq * (1 - vel_sys/c)
z = (HI_restfreq - freq_sys) / freq_sys
elif 'v_opt' in source.colnames:
vel_sys = source['v_opt']
z = vel_sys / c
elif 'v_app' in source.colnames:
vel_sys = source['v_app']
z = vel_sys / c
nhi = 1.104e+21 * (1 + z) ** 2 * sbr / bmaj / bmin
elif (bunit == 'Jy/beam*Hz') or (bunit == 'beam-1 Jy*Hz'):
nhi = 2.330e+20 * sbr / bmaj / bmin
freq_sys = source['freq']
z = (HI_restfreq.value - freq_sys) / freq_sys
nhi = 2.330e+20 * (1 + z) ** 4 * sbr / bmaj / bmin
else:
print("\tWARNING: Mom0 imag units are not Jy/beam*m/s or Jy/beam*Hz. Cannot convert to HI column density.")
nhi = sbr
print("\tWARNING: Mom0 image units are not Jy/beam*m/s or Jy/beam*Hz. Cannot convert to HI column density.")
nhi = sbr

if np.isfinite(nhi):
nhi_ofm = np.int(np.floor(np.log10(np.abs(nhi))))
nhi_ofm = np.int(np.floor(np.log10(np.abs(nhi))))
else:
nhi_ofm = 0
nhi_ofm = 0

nhi_label = '$N_\mathrm{{HI}}$ = {0:.1f} x $10^{{ {1:d} }}$ cm$^{{-2}}$'.format(nhi/10**nhi_ofm, nhi_ofm)
nhi_labels = '$N_\mathrm{{HI}}$ = $2^n$ x {0:.1f} x $10^{{ {1:d} }}$ cm$^{{-2}}$ ($n$=0,1,...)'.format(nhi/10**nhi_ofm, nhi_ofm)

return nhi, nhi_label, nhi_labels


Expand Down

0 comments on commit 7dd13b9

Please sign in to comment.