This jupyter notebook contains a simple reduction of the spectrophotometric standard star HD209290, observed with EMIR trhough the J, H and K grisms. Note that some steps have been dealt with quickly and in an approximate way.

The starting point are the images created by PyEmir, as explained in the tutorial:

https://guaix-ucm.github.io/pyemir-tutorials/tutorial_mos/ngc7798.html

The final results are response cuves, which are employed in that tutorial to proceed with the absolute flux calibration of NGC7798.

In [1]:
from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np

from numina.array.display.ximshow import ximshow
from numina.array.display.ximplotxy import ximplotxy

HD209290: spectrum extraction¶

Grism J¶

In [2]:
with fits.open('grismJ_hd209290_A-B.fits', mode='readonly') as hdulist:
    image_headerJ = hdulist[0].header
    image_dataJ = hdulist[0].data
    
# image dimensions
naxis2, naxis1 = image_dataJ.shape

# read wavelength calibration from header
crpix1 = image_headerJ['crpix1']
crval1 = image_headerJ['crval1']
cdelt1 = image_headerJ['cdelt1']
exptime = image_headerJ['exptime']
print('crpix1, crval1, cdelt1:', crpix1, crval1, cdelt1)
print('exptime:', exptime)

crpix1, crval1, cdelt1: 1.0 11200.0 0.77
exptime: 9.999289
In [3]:
ximshow(image_dataJ, image_bbox=(600,2750,1050,1110), 
        cbar_orientation='vertical',
        title='HD209290 (grism J)')
In [4]:
ximshow(image_dataJ, image_bbox=(600,2750,850,930), 
        cbar_orientation='vertical',
        title='HD209290 (grism J)')
In [5]:
# extract and coadd individual spectra by coadding rows
sp1 = np.sum(image_dataJ[1060:1100,], axis=0)
sp2 = -np.sum(image_dataJ[870:910,], axis=0)
# sum of two spectra
spectrumJ = sp1 + sp2
# compute counts/second
spectrumJ = spectrumJ / (2 * exptime)
In [6]:
# plot spectrum
fig, ax = plt.subplots()
waveJ = crval1 + (np.arange(1, naxis1 + 1) - crpix1) * cdelt1
ax.plot(waveJ, spectrumJ)
wvminJ = 11710
wvmaxJ = 13270
ax.set_xlim([wvminJ, wvmaxJ])
ax.set_xlabel('wavelength (Angstrom)')
ax.set_ylabel('count rate (ADU/s)')
ax.set_title('HD209290 (grism J)')
Out[6]:
Text(0.5, 1.0, 'HD209290 (grism J)')
In [7]:
# save spectrum into FITS file
hdu_sp = fits.PrimaryHDU(data=spectrumJ, header=image_headerJ)
hdu_sp.writeto('grismJ_hd209290.fits', overwrite=True)

Grism H¶

In [8]:
with fits.open('grismH_hd209290_A-B.fits', mode='readonly') as hdulist:
    image_headerH = hdulist[0].header
    image_dataH = hdulist[0].data
    
# image dimensions
naxis2, naxis1 = image_dataH.shape

# read wavelength calibration from header
crpix1 = image_headerH['crpix1']
crval1 = image_headerH['crval1']
cdelt1 = image_headerH['cdelt1']
exptime = image_headerH['exptime']
print('crpix1, crval1, cdelt1:', crpix1, crval1, cdelt1)
print('exptime:', exptime)

crpix1, crval1, cdelt1: 1.0 14500.0 1.22
exptime: 9.999289
In [9]:
ximshow(image_dataH, image_bbox=(580,2730,1050,1110), 
        cbar_orientation='vertical',
        title='HD209290 (grism H)')
In [10]:
ximshow(image_dataH, image_bbox=(580,2730,850,930), 
        cbar_orientation='vertical',
        title='HD209290 (grism H)')
In [11]:
# extract and coadd individual spectra by coadding rows
sp1 = np.sum(image_dataH[1060:1100,], axis=0)
sp2 = -np.sum(image_dataH[870:910,], axis=0)
# sum of two spectra
spectrumH = sp1 + sp2
# compute counts/second
spectrumH = spectrumH / (2 * exptime)
In [12]:
# plot spectrum
fig, ax = plt.subplots()
waveH = crval1 + (np.arange(1, naxis1 + 1) - crpix1) * cdelt1
ax.plot(waveH, spectrumH)
wvminH = 15280
wvmaxH = 17750
ax.set_xlim([wvminH, wvmaxH])
ax.set_xlabel('wavelength (Angstrom)')
ax.set_ylabel('count rate (ADU/s)')
ax.set_title('HD209290 (grism H)')
Out[12]:
Text(0.5, 1.0, 'HD209290 (grism H)')
In [13]:
# save spectrum into FITS file
hdu_sp = fits.PrimaryHDU(data=spectrumH, header=image_headerH)
hdu_sp.writeto('grismH_hd209290.fits', overwrite=True)

Grism K¶

In [14]:
with fits.open('grismK_hd209290_A-B.fits', mode='readonly') as hdulist:
    image_headerK = hdulist[0].header
    image_dataK = hdulist[0].data
    
# image dimensions
naxis2, naxis1 = image_dataK.shape

# read wavelength calibration from header
crpix1 = image_headerK['crpix1']
crval1 = image_headerK['crval1']
cdelt1 = image_headerK['cdelt1']
exptime = image_headerK['exptime']
print('crpix1, crval1, cdelt1:', crpix1, crval1, cdelt1)
print('exptime:', exptime)

crpix1, crval1, cdelt1: 1.0 19100.0 1.73
exptime: 9.999289
In [15]:
ximshow(image_dataK, image_bbox=(630,2780,1050,1110), 
        cbar_orientation='vertical',
        title='HD209290 (grism K)')
In [16]:
ximshow(image_dataK, image_bbox=(630,2780,850,930), 
        cbar_orientation='vertical',
        title='HD209290 (grism K)')
In [17]:
# extract and coadd individual spectra by coadding rows
sp1 = np.sum(image_dataK[1060:1100,], axis=0)
sp2 = -np.sum(image_dataK[870:910,], axis=0)
# sum of two spectra
spectrumK = sp1 + sp2
# compute counts/second
spectrumK = spectrumK / (2 * exptime)
In [18]:
# plot spectrum
fig, ax = plt.subplots()
waveK = crval1 + (np.arange(1, naxis1 + 1) - crpix1) * cdelt1
ax.plot(waveK, spectrumK)
wvminK = 20280
wvmaxK = 23810
ax.set_xlim([wvminK, wvmaxK])
ax.set_xlabel('wavelength (Angstrom)')
ax.set_ylabel('count rate (ADU/s)')
ax.set_title('HD209290 grism K')
Out[18]:
Text(0.5, 1.0, 'HD209290 grism K')
In [19]:
# save spectrum into FITS file
hdu_sp = fits.PrimaryHDU(data=spectrumK, header=image_headerK)
hdu_sp.writeto('grismK_hd209290.fits', overwrite=True)

HD209290: flux calibration¶

The flux calibrated spectrum of this star can be found in the web page of the NASA Infrared Telescope Facility (IRTF), in particular see the file: http://irtfweb.ifa.hawaii.edu/~spex/IRTF_Spectral_Library/Data/M0.5V_HD209290.txt

In [20]:
# read absolute flux calibration for the star
table_flux_tabulated = np.genfromtxt(
    "http://irtfweb.ifa.hawaii.edu/~spex/IRTF_Spectral_Library/Data/M0.5V_HD209290.txt"
)
print(table_flux_tabulated)

[[ 8.074530e-01  1.484350e-11  3.487900e-14]
 [ 8.076550e-01  1.492060e-11  3.509390e-14]
 [ 8.078580e-01  1.494630e-11  3.503310e-14]
 ...
 [ 5.422893e+00 -9.990000e+02 -9.990000e+02]
 [ 5.423898e+00 -9.990000e+02 -9.990000e+02]
 [ 5.424903e+00 -9.990000e+02 -9.990000e+02]]
In [21]:
# Wavelength (microns): convert to Angstrom
xtab = table_flux_tabulated[:,0] * 10000
# Flux Density (W m-2 um-1)
ytab = table_flux_tabulated[:,1]
In [22]:
# resample previous spectrum to sampling of observed spectrum
from scipy.interpolate import interp1d
funinterp = interp1d(xtab, ytab, kind="linear")

# convolve to match spectral resolution
from scipy.ndimage import gaussian_filter
yfluxJ = funinterp(waveJ)
sigma_tab_J = 0  # by 'eye'
if sigma_tab_J > 0:
    yfluxJ = gaussian_filter(yfluxJ, sigma_tab_J)

yfluxH = funinterp(waveH)
sigma_tab_H = 3  # by 'eye'
if sigma_tab_H > 0:
    yfluxH = gaussian_filter(yfluxH, sigma_tab_H)

yfluxK = funinterp(waveK)
sigma_tab_K = 3  # by 'eye'
if sigma_tab_K > 0:
    yfluxK = gaussian_filter(yfluxK, sigma_tab_K)
In [23]:
# plot observed spectrum and tabulated data
for grism, wvmin, wvmax, wave, spectrum, yflux in zip(['J', 'H', 'K'],
                                                      [wvminJ, wvminH, wvminK],
                                                      [wvmaxJ, wvmaxH, wvmaxK],
                                                      [waveJ, waveH, waveK],
                                                      [spectrumJ, spectrumH, spectrumK],
                                                      [yfluxJ, yfluxH, yfluxK]):
    fig, ax = plt.subplots()
    ax.plot(wave, spectrum, 'C0-', label='observed spectrum')
    ax.set_xlabel('Wavelength (Angstrom)')
    ax.set_ylim([0, spectrum.max()*1.05])
    ax.set_xlim([wvmin, wvmax])
    ax.set_ylabel('count rate (ADU/s)', color='C0')
    ax.legend(loc=4)    
    ax2 = ax.twinx()
    ax2.plot(wave, yflux, 'C1-', label='tabulated spectrum')
    ax2.set_ylim([0, yflux.max()*1.05])
    ax2.set_xlabel('wavelength (Angstrom)')
    ax2.set_ylabel(r'tabulated $F_\lambda$ (W m$^{-2} \mu\rm{m}^{-1}$)', color='C1')
    ax2.set_title('HD209290 (grism ' + grism + ')')
    ax2.legend(loc=3)

Since the spectral resolution of the tabulated spectrum is worse than the one in the observed spectrum, we convolve the observed spectrum (using gaussian kernel).

In [24]:
sigma_obs_J = 3  # by 'eye'
if sigma_obs_J > 0:
    yfluxJ = gaussian_filter(yfluxJ, sigma=sigma_obs_J)
    
sigma_obs_H = 0  # by 'eye'
if sigma_obs_H > 0:
    yfluxH = gaussian_filter(yfluxH, sigma=sigma_obs_H)
    
sigma_obs_K = 3  # by 'eye'
if sigma_obs_K > 0:
    yfluxK = gaussian_filter(yfluxK, sigma=sigma_obs_K)

# repeat previous plot
for grism, wvmin, wvmax, wave, spectrum, yflux in zip(['J', 'H', 'K'],
                                                      [wvminJ, wvminH, wvminK],
                                                      [wvmaxJ, wvmaxH, wvmaxK],
                                                      [waveJ, waveH, waveK],
                                                      [spectrumJ, spectrumH, spectrumK],
                                                      [yfluxJ, yfluxH, yfluxK]):
    fig, ax = plt.subplots()
    ax.plot(wave, spectrum, 'C0-', label='observed spectrum')
    ax.set_xlabel('wavelength (Angstrom)')
    ax.set_ylim([0, spectrum.max()*1.05])
    ax.set_xlim([wvmin, wvmax])
    ax.set_ylabel('count rate (ADU/s)', color='C0')
    ax.legend(loc=4)    
    ax2 = ax.twinx()
    ax2.plot(wave, yflux, 'C1-', label='tabulated spectrum')
    ax2.set_ylim([0, yflux.max()*1.05])
    ax2.set_xlabel('wavelength (Angstrom)')
    ax2.set_ylabel(r'tabulated $F_\lambda$ (W m$^{-2} \mu\rm{m}^{-1}$)', color='C1')
    ax2.set_title('HD209290 (grism ' + grism + ')')
    ax2.legend(loc=3)

Download the following atmospheric transmission file: http://nartex.fis.ucm.es/~ncl/emir/skycalc_transmission_R20000.txt

In [25]:
# preliminary response curve
yratioJ = spectrumJ / yfluxJ
yratioH = spectrumH / yfluxH
yratioK = spectrumK / yfluxK

# read telluric transmission
telluric_tabulated = np.genfromtxt("skycalc_transmission_R20000.txt")
xtab = telluric_tabulated[:,0] * 10  # convert from nm to Angstrom
ytab = telluric_tabulated[:,1]
# resample previous spectrum to sampling of observed spectrum
funinterp = interp1d(xtab, ytab, kind="linear")
ytelluricJ = funinterp(waveJ)
ytelluricH = funinterp(waveH)
ytelluricK = funinterp(waveK)
# broaden resolution
sigma_tel_J = 3  # by 'eye'
if sigma_tel_J > 0:
    ytelluricJ = gaussian_filter(ytelluricJ, sigma=sigma_tel_J)
sigma_tel_H = 0  # by 'eye'
if sigma_tel_H > 0:
    ytelluricH = gaussian_filter(ytelluricH, sigma=sigma_tel_H)
sigma_tel_K = 3  # by 'eye'
if sigma_tel_K > 0:
    ytelluricK = gaussian_filter(ytelluricK, sigma=sigma_tel_K)

# plot
for grism, wvmin, wvmax, wave, yratio, ytelluric in zip(['J', 'H', 'K'],
                                                      [wvminJ, wvminH, wvminK],
                                                      [wvmaxJ, wvmaxH, wvmaxK],
                                                      [waveJ, waveH, waveK],
                                                      [yratioJ, yratioH, yratioK],
                                                      [ytelluricJ, ytelluricH, ytelluricK]):
    fig, ax = plt.subplots()
    ax.plot(wave, yratio, 'C2', label='preliminary response curve')
    ax.set_xlim([wvmin, wvmax])
    ax.set_xlabel('wavelength (Angstrom)')
    ax.set_ylabel('preliminary response curve\n' + r'(ADU s$^{-1}$) / (W m$^{-2} \mu\rm{m}^{-1})$', color='g')
    ax.set_title('HD209290 (grism ' + grism +')')
    ax.legend(loc=3)
    ax2 = ax.twinx()
    ax2.plot(wave, ytelluric, 'C3', label='skycal_transmission')
    ax2.set_xlim([wvmin, wvmax])
    ax2.set_xlabel('wavelength (Angstrom)')
    ax2.set_ylabel('telluric transmission', color='C3')
    ax2.legend(loc=4)
In [26]:
# save preliminary response curve
hdu_sp = fits.PrimaryHDU(data=yratioJ, header=image_headerJ)
hdu_sp.writeto('grismJ_preliminary_response_curve_hd209290.fits', overwrite=True)

hdu_sp = fits.PrimaryHDU(data=yratioH, header=image_headerH)
hdu_sp.writeto('grismH_preliminary_response_curve_hd209290.fits', overwrite=True)

hdu_sp = fits.PrimaryHDU(data=yratioK, header=image_headerK)
hdu_sp.writeto('grismK_preliminary_response_curve_hd209290.fits', overwrite=True)
In [27]:
# ratio between response curve and atmosphere transmission
yratioplusJ = yratioJ / ytelluricJ
yratioplusH = yratioH / ytelluricH
yratioplusK = yratioK / ytelluricK

# smoothed version of previous curve
from scipy.signal import savgol_filter
yratiosmoothJ = savgol_filter(yratioplusJ, window_length=101, polyorder=2)
yratiosmoothH = savgol_filter(yratioplusH, window_length=101, polyorder=2)
yratiosmoothK = savgol_filter(yratioplusK, window_length=101, polyorder=2)

# plot
for grism, wvmin, wvmax, wave, yratioplus, yratiosmooth in zip(['J', 'H', 'K'],
                                                               [wvminJ, wvminH, wvminK],
                                                               [wvmaxJ, wvmaxH, wvmaxK],
                                                               [waveJ, waveH, waveK],
                                                               [yratioplusJ, yratioplusH, yratioplusK],
                                                               [yratiosmoothJ, yratiosmoothH, yratiosmoothK]):
    fig, ax = plt.subplots()
    ax.plot(wave, yratioplus, 'C4', label='smoothed curve')
    ax.set_xlim([wvmin, wvmax])
    ax.set_xlabel('wavelength (Angstrom)')
    ax.set_ylabel('response curve\n' + r'(ADU s$^{-1}$) / (W m$^{-2} \mu\rm{m}^{-1})$', color='k')
    ax.set_title('HD209290 (grism ' + grism + ')')
    ax.plot(wave, yratiosmooth, 'k', label='skycal_transmission')
    ax.legend(loc=3)
In [28]:
# save smoothed response curve
hdu_sp = fits.PrimaryHDU(data=yratiosmoothJ, header=image_headerJ)
hdu_sp.writeto('grismJ_response_curve_hd209290.fits', overwrite=True)

hdu_sp = fits.PrimaryHDU(data=yratiosmoothH, header=image_headerH)
hdu_sp.writeto('grismH_response_curve_hd209290.fits', overwrite=True)

hdu_sp = fits.PrimaryHDU(data=yratiosmoothK, header=image_headerK)
hdu_sp.writeto('grismK_response_curve_hd209290.fits', overwrite=True)
In [ ]: