This notebook outlines an example to optimize the ramp settings for a few different types of observations.
In these types of optimizations, we must consider observations constraints such as saturation levels, SNR requirements, and limits on acquisition time.
Note: The reported acquisition time does not include obsevatory and instrument-level overheads, such as slew times, filter changes, script compilations, etc. It only includes detector readout times (including reset frames).
In [1]:
# Makes print and division act like Python 3
from __future__ import print_function, division
# Import the usual libraries
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
# Enable inline plotting at lower left
%matplotlib inline
from IPython.display import display, Latex, clear_output
In [2]:
import pynrc
from pynrc import nrc_utils
from pynrc.nrc_utils import S
from pynrc.pynrc_core import table_filter
pynrc.setup_logging('WARNING', verbose=False)
from astropy.table import Table
We want to observe an M-Dwarf companion (K=18 mag) in the vicinity of a brighter F0V (K=13 mag) in the F430M filter. Assume the M-Dwarf flux is not significantly impacted by the brighter PSF (ie., in the background limited regime). In this scenario, the F0V star will saturate much more quickly compared to the fainter companion, so it limits which ramp settings we can use.
We will test a couple different types of observations (direct imaging vs coronagraphy).
In [3]:
# Get stellar spectra and normalize at K-Band
# The stellar_spectrum convenience function creates a Pysynphot spectrum
bp_k = S.ObsBandpass('k')
sp_M2V = pynrc.stellar_spectrum('M2V', 18, 'vegamag', bp_k)#, catname='ck04models')
sp_F0V = pynrc.stellar_spectrum('F0V', 13, 'vegamag', bp_k)#, catname='ck04models')
In [4]:
# Initiate a NIRCam observation
nrc = pynrc.NIRCam('F430M', wind_mode='WINDOW', xpix=160, ypix=160)
In [5]:
# Set some observing constraints
# Let's assume we want photometry on the primary to calibrate the M-Dwarf for direct imaging
# - Set well_frac_max=0.75
# Want a SNR~100 in the F430M filter
# - Set snr_goal=100
res = nrc.ramp_optimize(sp_M2V, sp_bright=sp_F0V, snr_goal=100, well_frac_max=0.75, verbose=True)
In [6]:
# Print the Top 2 settings for each readout pattern
res2 = table_filter(res, 2)
print(res2)
In [7]:
# Do the same thing, but for coronagraphic mask instead
#pynrc.setup_logging('DEBUG', verbose=False)
nrc = pynrc.NIRCam('F430M', mask='MASK430R', pupil='CIRCLYOT',
wind_mode='WINDOW', xpix=320, ypix=320)
# We assume that longer ramps will give us the best SNR for time
patterns = ['MEDIUM8', 'DEEP8']
res = nrc.ramp_optimize(sp_M2V, sp_bright=sp_F0V, snr_goal=100,
patterns=patterns, even_nints=True)
# Take the Top 2 settings for each readout pattern
res2 = table_filter(res, 2)
print(res2)
In [8]:
# Get stellar spectra and normalize at K-Band
# The stellar_spectrum convenience function creates a Pysynphot spectrum
bp_k = pynrc.bp_2mass('ks')
sp_G0V = pynrc.stellar_spectrum('G0V', 4, 'vegamag', bp_k)
# Choose a representative planet spectrum
planet = pynrc.planets_sb12(atmo='hy3s', mass=8, age=200, entropy=8, distance=17.5)
sp_pl = planet.export_pysynphot()
# Renormalize to F360M = 18.8
bp_l = pynrc.read_filter('F360M') #
sp_pl = sp_pl.renorm(18.8, 'vegamag', bp_l)
In [9]:
# Initiate a NIRCam observation
nrc = pynrc.NIRCam('F444W', pupil='CIRCLYOT', mask='MASK430R', wind_mode='WINDOW', xpix=320, ypix=320)
In [10]:
# Set even_nints=True assume 2 roll angles
res = nrc.ramp_optimize(sp_pl, sp_bright=sp_G0V, tacq_max=3600, tacq_frac=0.05,
even_nints=True, verbose=True)
In [11]:
# Take the Top 2 settings for each readout pattern
res2 = table_filter(res, 2)
print(res2)
In [12]:
# The DEEP and MEDIUMs are very similar.
# Let's go with MEDIUM2 for more GROUPS & INTS
# and slightly better efficiency over MEDIUM8
nrc.update_detectors(read_mode='MEDIUM2', ngroup=10, nint=36)
keys = list(nrc.multiaccum_times.keys())
keys.sort()
for k in keys:
print("{:<10}: {: 12.5f}".format(k, nrc.multiaccum_times[k]))
In [13]:
# Background sensitivity (5 sigma)
sens_dict = nrc.sensitivity(nsig=5, units='vegamag', verbose=True)
In [14]:
# M9V star at K=12 mags
bp_k = S.ObsBandpass('k')
sp_M9V = pynrc.stellar_spectrum('M9V', 12, 'vegamag', bp_k)
In [15]:
nrc = pynrc.NIRCam('F444W', pupil='GRISM0', wind_mode='STRIPE', ypix=64)
In [16]:
# Set a minimum of 10 integrations to be robust against cosmic rays
# Also set a minimum of 10 groups for good ramp sampling
res = nrc.ramp_optimize(sp_M9V, snr_goal=100, nint_min=10, ng_min=10, verbose=True)
In [17]:
# Print the Top 2 settings for each readout pattern
res2 = table_filter(res, 2)
print(res2)
In [18]:
# Let's say we choose SHALLOW4, NGRP=10, NINT=10
# Update detector readout
nrc.update_detectors(read_mode='SHALLOW4', ngroup=10, nint=10)
keys = list(nrc.multiaccum_times.keys())
keys.sort()
for k in keys:
print("{:<10}: {: 12.5f}".format(k, nrc.multiaccum_times[k]))
In [19]:
# Print final wavelength-dependent SNR
# For spectroscopy, the snr_goal is the median over the bandpass
snr_dict = nrc.sensitivity(sp=sp_M9V, forwardSNR=True, units='mJy', verbose=True)
Mock observed spectrum
Create a series of ramp integrations based on the current NIRCam settings. The gen_exposures() function creates a series of mock observations in raw DMS format by default. By default, it's point source objects centered in the observing window.
In [20]:
# Ideal spectrum and wavelength solution
wspec, imspec = nrc.gen_psf(sp=sp_M9V)
# Resize to detector window
nx = nrc.det_info['xpix']
ny = nrc.det_info['ypix']
In [21]:
# Create a series of ramp integrations based on the current NIRCam settings
# Output is 10 HDULists
im_slope = imspec + nrc.bg_zodi()
res = nrc.gen_exposures(im_slope=im_slope, return_results=True, targ_name='sp_M9V')
In [22]:
header = res[0]['PRIMARY'].header
tvals = (np.arange(header['NGROUPS'])+1) * header['TGROUP']
slope_list = []
for hdul in res:
header = hdul['PRIMARY'].header
data = hdul['SCI'].data
ref = pynrc.ref_pixels.NRC_refs(data, header, DMS=True, do_all=False)
ref.calc_avg_amps()
ref.correct_amp_refs()
# Linear fit to determine slope image
cf = nrc_utils.jl_poly_fit(tvals, ref.data, deg=1)
slope_list.append(cf[1])
# Create a master averaged slope image
slopes_all = np.array(slope_list)
slope_sim = slopes_all.mean(axis=0) * nrc.Detectors[0].gain
In [23]:
# Expand wspec to nx (fill value of 0)
# Then shrink to a size excluding wspec=0
# This assumes simulated spectrum is centered
wspec = nrc_utils.pad_or_cut_to_size(wspec, nx)
ind = wspec>0
# Estimate background emission and subtract from slope_sim
bg = np.median(slope_sim[:,~ind])
slope_sim = slope_sim[:,ind] - bg
wspec = wspec[ind]
In [24]:
# Extract 2 spectral x 5 spatial pixels
# First, cut out the central 5 pixels
slope_sub = nrc_utils.pad_or_cut_to_size(slope_sim, (5,slope_sim.shape[1]))
slope_sub_ideal = nrc_utils.pad_or_cut_to_size(imspec, (5,imspec.shape[1]))
# Sum along the spatial axis
spec = slope_sub.sum(axis=0)
spec_ideal = slope_sub_ideal.sum(axis=0)
spec_ideal_rebin = nrc_utils.frebin(spec_ideal, scale=0.5, total=False)
# Build a quick RSRF from extracted ideal spectral slope
sp_M9V.convert('mjy')
rsrf = spec_ideal / sp_M9V.sample(wspec*1e4)
# Rebin along spectral direction
wspec_rebin = nrc_utils.frebin(wspec, scale=0.5, total=False)
spec_rebin_cal = nrc_utils.frebin(spec/rsrf, scale=0.5, total=False)
In [25]:
# Expected noise per extraction element
snr_interp = np.interp(wspec_rebin, snr_dict['wave'], snr_dict['snr'])
_spec_rebin = spec_ideal_rebin / snr_interp
_spec_rebin_cal = _spec_rebin / nrc_utils.frebin(rsrf, scale=0.5, total=False)
In [26]:
fig, ax = plt.subplots(1,1, figsize=(12,8))
ax.plot(sp_M9V.wave/1e4, sp_M9V.flux, label='Input Spectrum')
ax.plot(wspec_rebin, spec_rebin_cal, alpha=0.7, label='Extracted Observation')
ax.errorbar(wspec_rebin, spec_rebin_cal, yerr=_spec_rebin_cal, zorder=3,
fmt='none', label='Expected Error Bars', alpha=0.7)
ax.set_ylim([0,10])
ax.set_xlim([3.7,5.1])
ax.set_xlabel('Wavelength ($\mu m$)')
ax.set_ylabel('Flux (mJy)')
ax.set_title('Simulated Spectrum')
ax.legend(loc='upper right');
In [27]:
nrc = pynrc.NIRCam('F322W2', pupil='GRISM0', wind_mode='STRIPE', ypix=64)
In [28]:
# K6V star at K=8.4 mags
bp_k = S.ObsBandpass('k')
sp_K6V = pynrc.stellar_spectrum('K6V', 8.4, 'vegamag', bp_k)
In [29]:
# Constraints
well = 0.5 # Keep well below 50% full
tacq = 2.1*3600. # 2.1 hour transit duration
ng_max = 30 # Transit spectroscopy allows for up to 30 groups per integrations
nint_max = int(1e6) # Effectively no limit on number of integrations
# Let's bin the spectrum to R~100
# dw_bin is a passable parameter for specifiying spectral bin sizes
R = 100
dw_bin = (nrc.bandpass.avgwave() / 10000) / R
In [35]:
res = nrc.ramp_optimize(sp_K6V, tacq_max=tacq, nint_max=nint_max,
ng_min=10, ng_max=ng_max, well_frac_max=well,
dw_bin=dw_bin, verbose=True)
In [31]:
# Print the Top 2 settings for each readout pattern
res2 = table_filter(res, 2)
print(res2)
In [32]:
# Even though BRIGHT1 has a slight efficiency preference over RAPID
# and BRIGHT2, we decide to choose RAPID, because we are convinced
# that saving all data (and no coadding) is a better option.
# If APT informs you that the data rates or total data shorage is
# an issue, you can select one of the other options.
# Update to RAPID, ngroup=30, nint=715 and plot PPM
nrc.update_detectors(read_mode='RAPID', ngroup=30, nint=715)
snr_dict = nrc.sensitivity(sp=sp_K6V, dw_bin=dw_bin, forwardSNR=True, units='Jy')
wave = np.array(snr_dict['wave'])
snr = np.array(snr_dict['snr'])
# Let assume bg subtraction of something with similar noise
snr /= np.sqrt(2.)
ppm = 1e6 / snr
# NOTE: We have up until now neglected to include a "noise floor"
# which represents the expected minimum achievable ppm from
# unknown systematics. To first order, this can be added in
# quadrature to the calculated PPM.
noise_floor = 30 # in ppm
ppm_floor = np.sqrt(ppm**2 + noise_floor**2)
In [33]:
plt.plot(wave, ppm, marker='o', label='Calculated PPM')
plt.plot(wave, ppm_floor, marker='o', label='PPM + Noise Floor')
plt.xlabel('Wavelength ($\mu m$)')
plt.ylabel('Noise Limit (PPM)')
plt.xlim([2.4,4.1])
plt.ylim([20,100])
plt.legend()
Out[33]:
In [34]:
# Detection bandpass is F200W
nrc = pynrc.NIRCam('F200W')
# Flat spectrum (in photlam) with ABMag = 25 in the NIRCam bandpass
sp = pynrc.stellar_spectrum('flat', 25, 'abmag', nrc.bandpass)
In [35]:
res = nrc.ramp_optimize(sp, is_extended=True, tacq_max=10000, tacq_frac=0.05, verbose=True)
In [36]:
# Print the Top 2 settings for each readout pattern
res2 = table_filter(res, 2)
print(res2)
In [37]:
# MEDIUM8 10 10 looks like a good option
nrc.update_detectors(read_mode='MEDIUM8', ngroup=10, nint=10, verbose=True)
In [38]:
# Calculate flux/mag for various nsigma detection limits
tbl = Table(names=('Sigma', 'Point (nJy)', 'Extended (nJy/asec^2)',
'Point (AB Mag)', 'Extended (AB Mag/asec^2)'))
tbl['Sigma'].format = '.0f'
for k in tbl.keys()[1:]:
tbl[k].format = '.2f'
for sig in [1,3,5,10]:
snr_dict1 = nrc.sensitivity(nsig=sig, units='nJy', verbose=False)
snr_dict2 = nrc.sensitivity(nsig=sig, units='abmag', verbose=False)
tbl.add_row([sig, snr_dict1[0]['sensitivity'], snr_dict1[1]['sensitivity'],
snr_dict2[0]['sensitivity'], snr_dict2[1]['sensitivity']])
In [39]:
tbl
Out[39]:
In [ ]: