Spatial smoothing: beta=1 error

Issue #2 wontfix
Maria Werhahn created an issue

Doing spatial smoothing with beta=1, wavelength smoothing with F 250 and generating s2n-cube:

lsd_cc_spatial.py -i \
median_filtered_W201_DATACUBE-HDFS-v1.24_incl_effnoise_blowup.fits -S 1 -N 2 \
--beta 1 -p0 1.2 -p1=-1.24e-4 -p2 6.3e-9 -o \
spatial_smoothed_beta1_median_filtered_W201_DATACUBE-HDFS-v1.24_incl_effnoise_blowup.fits

lsd_cc_spectral.py -i \
spatial_smoothed_beta1_median_filtered_W201_DATACUBE-HDFS-v1.24_incl_effnoise_blowup.fits \
-F 250 -o \
wavelength_smooth_F250_spatial_smoothed_beta1_median_filtered_W201_DATACUBE-HDFS-v1.24_incl_effnoise_blowup.fits

s2n-cube.py -i \
wavelength_smooth_F250_spatial_smoothed_beta1_median_filtered_W201_DATACUBE-HDFS-v1.24_incl_effnoise_blowup.fits \
-o s2n_v250_beta1_skpsf_blowup.fits

Saving the maximum S/N values of the Lyman Alpha Emitters (Jarle Table) in a npz-file:

import pylab as p
import numpy as np
from astropy.io import fits

# input datacube
#sn_cube_hdu = fits.open('s2n_v300_beta1_skpsf_blowup.fits')
sn_cube_hdu = fits.open('./s2n_v250_beta1_skpsf_blowup.fits')
sn_cube = sn_cube_hdu[0].data
sn_cube_header = sn_cube_hdu[0].header
cdelt = sn_cube_header['CD3_3']
crval = sn_cube_header['CRVAL3']
crpix = sn_cube_header['CRPIX3']

# input from table
jarle_table_hdu = fits.open('./pfit-results-v03p1.fits')
jarle_table = jarle_table_hdu[1].data

table_x = jarle_table['X']
table_y = jarle_table['Y']
redshift = jarle_table['Z']
lya_flux_all = jarle_table['LY_ALPHA_FLUX']
lya_flux_err_all = jarle_table['LY_ALPHA_FLUX_ERR']
redshift_confidence = jarle_table['Z_CONFIDENCE']
ids = jarle_table['ID']

# Selektor
lya_em_sel = (lya_flux_all > 0) & (redshift_confidence > 5)

lya_x = table_x[lya_em_sel]
lya_y = table_y[lya_em_sel]
lya_redshift = redshift[lya_em_sel]
lya_flux = lya_flux_all[lya_em_sel]
lya_flux_err = lya_flux_err_all[lya_em_sel]
lya_ids = ids[lya_em_sel]

lya_0 = 1215.67
lya_obs = (1 + lya_redshift) * lya_0
lya_z = (lya_obs - crval) / cdelt + crpix

# now searchin for maximum in 10x10x10 cubes centered on
# lya_x,lya_y,ly_z
lya_xmin = lya_x - 5
lya_xmax = lya_x + 5
lya_ymin = lya_y - 5
lya_ymax = lya_y + 5
lya_zmin = lya_z - 5
lya_zmax = lya_z + 5

lya_sn_max = np.zeros_like(lya_xmin)

for i in xrange(len(lya_xmin)):
    lya_sn_max[i] = sn_cube[lya_zmin[i]:lya_zmax[i],
                            lya_ymin[i]:lya_ymax[i],
                            lya_xmin[i]:lya_xmax[i]].max()


#save as npz
np.savez("beta10.npz",lya_ids=lya_ids,ids=ids,lya_sn_max=lya_sn_max, lya_flux=lya_flux,lya_flux_err=lya_flux_err)

loading the file in ipython (import numpy as np):

In [296]:beta10=np.load('beta10.npz')

In [297]: beta10['lya_sn_max'] Out[297]: array([ nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan], dtype=float32)

Comments (3)

  1. Christian Herenz repo owner

    For beta 1 everything in the moffat kernel will be zero, since

    norm = (beta - 1) / (m.pi * m.pow(R,2))
    

    The area under the Moffat is 1 / norm.... I.e. infinity for beta == 1. Thus beta == 1 is an impossible value for the Moffat.

  2. Log in to comment