Basic Usage¶
This tutorial walks through the basic usage of the pynrc
package to calculate sensitivities and saturation limits for NIRCam in a variety of modes.
[1]:
# Import the usual libraries
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
# Enable inline plotting at lower left
%matplotlib inline
Getting Started¶
We assume you have already installed pynrc
as outlined in the documentation.
[2]:
# import main module
import pynrc
from pynrc import nrc_utils
Log messages for pynrc
follow the same the logging functionality included in webbpsf
. Logging levels include DEBUG
, INFO
, WARN
, and ERROR
.
[3]:
pynrc.setup_logging()
pyNRC log messages of level INFO and above will be shown.
pyNRC log outputs will be directed to the screen.
If you get tired of the INFO
level messages, simply type:
pynrc.setup_logging('WARN', verbose=False)
First NIRCam Observation¶
The basic NIRCam object consists of all the instrument settings one would specify for a JWST observation, including filter, pupil, and coronagraphic mask selections along with detector subarray settings and ramp sampling cadence (i.e., MULTIACCUM
).
The NIRCam class makes use of high order polynomial coefficient maps to quickly generate large numbers of monochromatic PSFs that can be convolved with arbitrary spectra and collapsed into a final broadband PSF (or dispersed with NIRCam’s slitless grisms). The PSF coefficients are calculated from a series of WebbPSF monochromatic PSFs and saved to disk. These polynomial coefficients are further modifed based on focal plane position and drift in the wavefront error relative to nominal OPD mao.
There are a multitude of posssible keywords one can pass upon initialization, including certain detector settings and PSF generation parameters. If not passed initially, then defaults are assumed. The user can update these parameters at any time by either setting attributes directly (e.g., filter
, mask
, pupil
, etc.) along with using the update_detectors()
and update_psf_coeff()
methods.
For instance,
nrc = pynrc.NIRCam('F210M')
nrc.module = 'B'
nrc.update_detectors(read_mode='DEEP8', nint=10, ngroup=5)
is the same as:
nrc = pynrc.NIRCam('F210M', module='B', read_mode='DEEP8', nint=10, ngroup=5)
To start, we’ll set up a simple observation using the F430M
filter. Defaults will be populated for unspecified attributes such as module
, pupil
, mask
, etc.
Check the function docstrings for more detailed information
[4]:
nrc = pynrc.NIRCam(filter='F430M', fov_pix=33, oversample=4)
print('\nFilter: {}; Pupil Mask: {}; Image Mask: {}; Module: {}'\
.format(nrc.filter, nrc.pupil_mask, nrc.image_mask, nrc.module))
[ webbpsf:INFO] NIRCam aperture name updated to NRCA1_FULL
[ webbpsf:INFO] NIRCam pixel scale switched to 0.063000 arcsec/pixel for the long wave channel.
[ webbpsf:INFO] NIRCam aperture name updated to NRCA5_FULL
[ pynrc:INFO] RAPID readout mode selected.
[ pynrc:INFO] Setting ngroup=1, nf=1, nd1=0, nd2=0, nd3=0.
[ pynrc:INFO] Initializing SCA 485/A5
[ pynrc:INFO] Suggested SIAF aperture name: NRCA5_FULL
[ pynrc:INFO] RAPID readout mode selected.
[ pynrc:INFO] Setting ngroup=1, nf=1, nd1=0, nd2=0, nd3=0.
[ pynrc:INFO] Initializing SCA 485/A5
[webbpsf_ext:INFO] Loading /Users/jarron/NIRCam/webbpsf_ext_data/psf_coeffs/NIRCam/LWA_F430M_CLEAR_NONE_pix33_os4_jsig5_r0.00_th+0.0_RevAAslice0_siwfe_distort_legendre.fits
[ pynrc:INFO] Suggested SIAF aperture name: NRCA5_FULL
Filter: F430M; Pupil Mask: None; Image Mask: None; Module: A
Keyword information for detector settings are stored in the det_info
dictionary. These cannot be modified directly, but instead are updated via the update_detectors()
methods.
[5]:
print('Detector Info Keywords:')
print(nrc.det_info)
Detector Info Keywords:
{'wind_mode': 'FULL', 'nout': 4, 'xpix': 2048, 'ypix': 2048, 'x0': 0, 'y0': 0, 'read_mode': 'RAPID', 'nint': 1, 'ngroup': 1, 'nf': 1, 'nd1': 0, 'nd2': 0, 'nd3': 0}
PSF settings are stored and modified in the same manner as WebbPSF (https://webbpsf.readthedocs.io/en/stable/usage.html) with a couple minor modifications: 1. The calculated field of view is specified as nrc.fov_pix
along with the nrc.oversample
attributes. 2. In addition, you can turn on/off distortions via the nrc.include_distortions
attribute. 3. The primary method for calculating PSFs has changed to nrc.calc_psf_from_coeff
.
You can quickly obtain a brief overview of important settings through the psf_info
dictionary property. These properties can also be updated via the nrc.update_psf_coeff()
function, which immediately generates a new set of PSF coefficients.
[6]:
print('PSF Info:')
print(nrc.psf_info)
PSF Info:
{'fov_pix': 33, 'oversample': 4, 'npsf': 7, 'ndeg': 6, 'include_si_wfe': True, 'include_distortions': True, 'jitter': 'gaussian', 'jitter_sigma': 0.005, 'offset_r': 0, 'offset_theta': 0, 'bar_offset': None, 'pupil': '/Users/jarron/NIRCam/webbpsf-data/jwst_pupil_RevW_npix1024.fits.gz', 'pupilopd': ('JWST_OTE_OPD_RevAA_prelaunch_predicted.fits.gz', 0)}
PSF coefficient information is stored in the psf_coeff
attribute. An associated header file exists in the psf_coeff_header
, showing all of the parameters used to generate that data. This data is accessed by many of the NIRCam class functions to generate PSFs with arbitrary wavelength weights, such as the calc_psf_from_coeff()
function.
[7]:
# Demonstrate the color difference of the PSF for different spectral types, same magnitude
sp_M0V = pynrc.stellar_spectrum('M0V', 10, 'vegamag', nrc.bandpass)
sp_A0V = pynrc.stellar_spectrum('A0V', 10, 'vegamag', nrc.bandpass)
# Generate oversampled PSFs (counts/sec)
hdul_M0V = nrc.calc_psf_from_coeff(sp=sp_M0V, return_oversample=True)
hdul_A0V = nrc.calc_psf_from_coeff(sp=sp_A0V, return_oversample=True)
psf_M0V = hdul_M0V[0].data
psf_A0V = hdul_A0V[0].data
fig, axes = plt.subplots(1,3, figsize=(12,4))
axes[0].imshow(psf_M0V**0.5)
axes[0].set_title('M0V PSF ({})'.format(nrc.filter))
axes[1].imshow(psf_A0V**0.5)
axes[1].set_title('A0V PSF ({})'.format(nrc.filter))
diff = psf_M0V - psf_A0V
minmax = np.abs(diff).max() / 2
axes[2].imshow(diff, cmap='RdBu', vmin=-minmax, vmax=minmax)
axes[2].set_title('Difference')
fig.tight_layout()
Bandpass information is stored in the bandpass
attribute and can be plotted with the convenience function plot_bandpass()
.
[8]:
fig, ax = plt.subplots(1,1)
nrc.plot_bandpass(ax=ax, color='C3')
fig.tight_layout()
1. Saturation Limits¶
One of the most basic functions is to determine the saturation limit of a CDS observation, so let’s try this for the current filter selection. Generally, saturation is considered to be 80% of the full well, but can go as high as 95%.
[9]:
# Turn off those pesky informational texts
pynrc.setup_logging('WARN', verbose=False)
# Configure the observation for CDS frames (ngroup=2)
# Print out frame and ramp information using verbose=True
nrc.update_detectors(ngroup=2, verbose=True)
New Ramp Settings
read_mode : RAPID
nf : 1
nd2 : 0
ngroup : 2
nint : 1
New Detector Settings
wind_mode : FULL
xpix : 2048
ypix : 2048
x0 : 0
y0 : 0
New Ramp Times
t_group : 10.737
t_frame : 10.737
t_int : 21.474
t_int_tot1 : 21.474
t_int_tot2 : 0.000
t_exp : 21.474
t_acq : 21.479
The sat_limits()
function returns a dictionary of results. There’s the option in include a Pysynphot spectrum, but if None is specificed then it defaults to a G2V star.
[10]:
# Set verbose=True to print results in a user-friendly manner
sat_lims = nrc.sat_limits(verbose=True)
# Dictionary information
print("\nDictionary Info:", sat_lims)
F430M Saturation Limit assuming G2V source (point source): 12.18 vegamag
F430M Saturation Limit assuming G2V source (extended): 8.48 vegamag/arcsec^2
Dictionary Info: ({'satlim': 12.179917081677289, 'units': 'vegamag', 'bp_lim': 'F430M', 'Spectrum': 'G2V'}, {'satlim': 8.478633456035855, 'units': 'vegamag/arcsec^2', 'bp_lim': 'F430M', 'Spectrum': 'G2V'})
By default, the function sat_limits()
uses a G2V stellar spectrum, but any arbritrary spectrum can be passed via the sp
keyword. In addition, using the bp_lim
keyword, you can use spectral information to determine the brightness in some other bandpass that saturates the source within the NIRCam filter.
[11]:
# Spectrum of an M0V star (not normalized)
sp_M0V = pynrc.stellar_spectrum('M0V')
# 2MASS Ks Bandpass
bp_k = pynrc.bp_2mass('K')
sat_lims = nrc.sat_limits(sp=sp_M0V, bp_lim=bp_k, verbose=True)
Ks-Band Saturation Limit for F430M assuming M0V source (point source): 12.31 vegamag
Ks-Band Saturation Limit for F430M assuming M0V source (extended): 8.61 vegamag/arcsec^2
Now, let’s get the same saturation limit assuming a 128x128 detector subarray (faster frame rate).
[12]:
nrc.update_detectors(wind_mode='WINDOW', xpix=128, ypix=128)
sat_lims = nrc.sat_limits(sp=sp_M0V, bp_lim=bp_k, verbose=True)
Ks-Band Saturation Limit for F430M assuming M0V source (point source): 7.89 vegamag
Ks-Band Saturation Limit for F430M assuming M0V source (extended): 4.18 vegamag/arcsec^2
You can also use the saturation_levels()
function to generate an image of a point source indicating the fractional well fill level.
[13]:
# Spectum of A0V star with Ks = 8 mag
sp = pynrc.stellar_spectrum('M0V', 8, 'vegamag', bp_k)
sat_levels = nrc.saturation_levels(sp, full_size=False, ngroup=nrc.det_info['ngroup'])
print('Max Well Fraction: {:.2f}'.format(sat_levels.max()))
Max Well Fraction: 0.72
[14]:
# Plot the well fill levels for each pixel
fig, axes = plt.subplots(1,2, figsize=(12,5))
for i,ax in enumerate(axes):
extent = 0.5 * nrc.psf_info['fov_pix'] * np.array([-1,1,-1,1])
if i==0:
cax = ax.imshow(sat_levels, extent=extent, vmin=0, vmax=1)
else:
norm = matplotlib.colors.LogNorm(vmin=0.001, vmax=1)
cax = ax.imshow(sat_levels, extent=extent, norm=norm)
ax.set_xlabel('Pixels')
ax.set_ylabel('Pixels')
ax.set_title('Well Fraction in {} of $K_s = 5$ M0V star'.format(nrc.filter))
cbar = fig.colorbar(cax, ax=ax)
cbar.set_label('Well Fill Fraction')
ax.tick_params(axis='both', color='white', which='both')
for k in ax.spines.keys():
ax.spines[k].set_color('white')
fig.tight_layout()
Information for slitless grism observations show wavelength-dependent results.
[15]:
nrc = pynrc.NIRCam(filter='F444W', pupil_mask='GRISM0', ngroup=2, wind_mode='STRIPE', ypix=128)
sat_lims = nrc.sat_limits(sp=sp_M0V, bp_lim=bp_k, verbose=True)
Ks-Band Saturation Limit for F444W assuming M0V source:
Wave Sat Limit (vegamag)
--------- -------------------
3.90 4.44
4.00 4.44
4.10 4.33
4.20 4.21
4.30 4.06
4.40 3.82
4.50 3.64
4.60 3.47
4.70 3.30
4.80 3.13
4.90 2.92
5.00 2.35
2. Sensitivity Limits¶
Similarly, we can determine sensitivity limits of point sources (and extended sources) for the defined instrument configuration. By default, the sensitivity()
function uses a flat spectrum. In this case, let’s find the sensitivities NIRCam can reach in a single ~1000sec integration with the F430M filter. Noise values will depend on the exact MULTIACCUM
settings.
[16]:
nrc = pynrc.NIRCam(filter='F430M')
nrc.update_detectors(read_mode='MEDIUM8', ngroup=10)
# The multiaccum_times attribute describes the various timing information
print(nrc.multiaccum_times)
{'t_frame': 10.73677, 't_group': 107.3677, 't_int': 1052.20346, 't_exp': 1052.20346, 't_acq': 1052.2087, 't_int_tot1': 1052.20346, 't_int_tot2': 0.0}
[17]:
sens = nrc.sensitivity(nsig=5, units='vegamag', verbose=True)
Point Source Sensitivity (5-sigma): 23.35 vegamag
Surface Brightness Sensitivity (5-sigma): 21.23 vegamag/arcsec^2
The sensitivity function also includes a keyword forwardSNR
, which allows the user to pass a normalized spectrum and estimate the SNR For some extraction aperture.
[18]:
sp = pynrc.stellar_spectrum('M0V', 20, 'vegamag', nrc.bandpass)
snr = nrc.sensitivity(sp=sp, forwardSNR=True, units='vegamag', verbose=True)
Point Source SNR (20.00 vegamag): 77.08 sigma
Surface Brightness SNR (20.00 vegamag/arcsec^2): 14.59 sigma
3. Ramp Optimization¶
Armed with these two basic functions, we can attempt to determine the best instrument settings to optimize for SNR and efficiency. In these types of optimizations, we must consider observational constraints such as saturation levels, SNR requirements, and limits on acquisition time.
Note: The reported acquisition times do 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 and Fast Row Resets).
For instance, we want to observe an M-Dwarf (K=18 mag) in the F430M filter. What is the most efficient configuration to obtain an SNR of 100?
[19]:
# Setup observation
nrc = pynrc.NIRCam(filter='F430M', wind_mode='WINDOW', xpix=160, ypix=160)
# Spectrum of an M2V star
bp_k = pynrc.bp_2mass('K')
sp_M0V = pynrc.stellar_spectrum('M0V', 18, 'vegamag', bp_k)
[20]:
# Run optimizer. Result is a sorted by efficiency.
tbl = nrc.ramp_optimize(sp_M0V, snr_goal=100, ng_min=5, nint_min=10, verbose=True)
BRIGHT1
BRIGHT2
DEEP2
DEEP8
MEDIUM2
MEDIUM8
RAPID
SHALLOW2
SHALLOW4
Pattern NGRP NINT t_int t_exp t_acq SNR Well eff
---------- ---- ---- --------- --------- --------- -------- -------- --------
DEEP8 5 10 24.52 245.20 248.04 115.7 0.005 7.344
MEDIUM8 7 11 18.95 208.42 211.54 101.1 0.004 6.949
MEDIUM8 7 12 18.95 227.37 230.78 105.6 0.004 6.949
MEDIUM2 9 10 22.85 228.48 231.32 101.7 0.005 6.687
MEDIUM2 9 11 22.85 251.33 254.46 106.7 0.005 6.687
SHALLOW4 10 18 13.65 245.76 250.87 100.3 0.003 6.334
SHALLOW4 10 19 13.65 259.41 264.81 103.1 0.003 6.334
DEEP2 5 11 22.85 251.33 254.46 99.2 0.005 6.218
... ... ... ... ... ... ... ... ...
RAPID 10 634 2.79 1766.58 1946.51 101.5 0.001 2.300
RAPID 10 635 2.79 1769.36 1949.58 101.6 0.001 2.300
RAPID 10 636 2.79 1772.15 1952.65 101.7 0.001 2.300
RAPID 10 637 2.79 1774.94 1955.72 101.7 0.001 2.300
RAPID 10 638 2.79 1777.72 1958.79 101.8 0.001 2.300
RAPID 10 639 2.79 1780.51 1961.86 101.9 0.001 2.300
RAPID 10 640 2.79 1783.30 1964.93 102.0 0.001 2.300
RAPID 10 641 2.79 1786.08 1968.00 102.1 0.001 2.300
Length = 55 rows
For a slightly more complicated scenario, consider an additional foreground source. In this scenario, the F0V star will saturate much more quickly compared to the fainter M2V, so it limits which ramp settings we may want to use (assuming we want unsaturated frames, which isn’t always necessarily true).
[21]:
sp_F0V = pynrc.stellar_spectrum('F0V', 10, 'vegamag', bp_k)
tbl = nrc.ramp_optimize(sp_M0V, sp_bright=sp_F0V, snr_goal=100, ng_min=5, nint_min=10,
well_frac_max=0.9, verbose=True)
BRIGHT1
BRIGHT2
DEEP2
DEEP8
MEDIUM2
MEDIUM8
RAPID
SHALLOW2
SHALLOW4
Pattern NGRP NINT t_int t_exp t_acq SNR Well eff
---------- ---- ---- --------- --------- --------- -------- -------- --------
RAPID 10 615 2.79 1713.64 1888.17 100.0 0.778 2.300
RAPID 10 616 2.79 1716.42 1891.24 100.0 0.778 2.300
RAPID 10 617 2.79 1719.21 1894.31 100.1 0.778 2.300
RAPID 10 618 2.79 1722.00 1897.38 100.2 0.778 2.300
RAPID 10 619 2.79 1724.78 1900.45 100.3 0.778 2.300
RAPID 10 620 2.79 1727.57 1903.52 100.4 0.778 2.300
RAPID 10 621 2.79 1730.35 1906.59 100.5 0.778 2.300
RAPID 10 622 2.79 1733.14 1909.66 100.5 0.778 2.300
... ... ... ... ... ... ... ... ...
BRIGHT1 6 699 3.07 2142.46 2340.84 101.6 0.856 2.099
BRIGHT1 6 700 3.07 2145.53 2344.19 101.6 0.856 2.099
BRIGHT1 6 701 3.07 2148.59 2347.54 101.7 0.856 2.099
BRIGHT1 6 702 3.07 2151.66 2350.89 101.8 0.856 2.099
BRIGHT1 6 703 3.07 2154.72 2354.23 101.8 0.856 2.099
BRIGHT1 6 704 3.07 2157.79 2357.58 101.9 0.856 2.099
BRIGHT1 6 705 3.07 2160.85 2360.93 102.0 0.856 2.099
BRIGHT1 6 706 3.07 2163.92 2364.28 102.1 0.856 2.099
Length = 85 rows
If there are no objections to saturating the bright source, then we can set the well_frac_max
parameter to something like 5 times the hard saturation limit. This allows for much more efficient exposure settings.
[22]:
tbl = nrc.ramp_optimize(sp_M0V, sp_bright=sp_F0V, snr_goal=100, ng_min=5, nint_min=10,
well_frac_max=5, verbose=True)
BRIGHT1
BRIGHT2
DEEP2
DEEP8
MEDIUM2
MEDIUM8
RAPID
SHALLOW2
SHALLOW4
Pattern NGRP NINT t_int t_exp t_acq SNR Well eff
---------- ---- ---- --------- --------- --------- -------- -------- --------
MEDIUM8 6 14 16.16 226.26 230.23 100.1 4.511 6.597
MEDIUM8 6 15 16.16 242.42 246.67 103.6 4.511 6.597
SHALLOW4 10 18 13.65 245.76 250.87 100.3 3.811 6.334
SHALLOW4 10 19 13.65 259.41 264.81 103.1 3.811 6.334
MEDIUM8 5 21 13.37 280.87 286.83 103.9 3.733 6.135
MEDIUM2 7 16 17.28 276.41 280.95 100.1 4.822 5.971
MEDIUM2 7 17 17.28 293.69 298.51 103.2 4.822 5.971
SHALLOW2 10 22 13.10 288.11 294.36 98.8 3.656 5.759
... ... ... ... ... ... ... ... ...
RAPID 10 634 2.79 1766.58 1946.51 101.5 0.778 2.300
RAPID 10 635 2.79 1769.36 1949.58 101.6 0.778 2.300
RAPID 10 636 2.79 1772.15 1952.65 101.7 0.778 2.300
RAPID 10 637 2.79 1774.94 1955.72 101.7 0.778 2.300
RAPID 10 638 2.79 1777.72 1958.79 101.8 0.778 2.300
RAPID 10 639 2.79 1780.51 1961.86 101.9 0.778 2.300
RAPID 10 640 2.79 1783.30 1964.93 102.0 0.778 2.300
RAPID 10 641 2.79 1786.08 1968.00 102.1 0.778 2.300
Length = 54 rows
[ ]: