Calculating WFC3 zeropoints with STSynphot¶
This notebook shows how to calculate photometric zeropoints using the python package stsynphot for any WFC3 detector, filter, date, or aperture. This is especially useful for calculating Vegamag zeropoints which require an input spectrum. The notebook is also useful for computing time-dependent WFC3/UVIS zeropoints for any observation date, as the values listed in WFC3 ISR 2021-04 are defined for the reference epoch. As of mid-2021, the WFC3/IR zeropoints are not time-dependent).
To install stsynphot, activate your conda environment in a bash shell and enter the command pip install stsynphot. More documentation on stsynphot is available here. Using stsynphot requires downloading the throughput curves for the HST instruments and optical path. One method of doing this is shown below. More information can be found here.
1. Downloading throughput tables and defining variables ¶
This section obtains the WFC3 throughput component tables for use with synphot
import os
#cmd_input = 'curl -O ftp://archive.stsci.edu/pub/hst/pysynphot/synphot1.tar.gz'
#os.system(cmd_input)
Once the files are downloaded, unpack the files and set the environment variable PYSYN_CDBS to the path of the unpacked files.
# os.environ['PYSYN_CDBS'] = '/YOUR/PATH/HERE/'
os.environ['PYSYN_CDBS'] = '/grp/hst/cdbs/' # for STScI staff
import numpy as np
import stsynphot as stsyn
from astropy.table import Table
from astropy.time import Time
from synphot import Observation
#Rather than downloading the entire calspec database (synphot6.tar.gz),
#we can point directly to the latest Vega spectrum which is required for computing VEGAMAG
vega_url = 'https://ssb.stsci.edu/trds/calspec/alpha_lyr_stis_010.fits'
stsyn.Vega = stsyn.spectrum.SourceSpectrum.from_file(vega_url)
2. Setting up the 'obsmode' string¶
Parameters to set in the obsmode string include: 1.) detector, 2.) filter, 3.) observation date (UVIS only), and 4.) aperture size (in arcsec).
Note that a 6.0" aperture is considered to be 'infinite', thus containing all of the flux. The zeropoints posted on the WFC3 website are calculated for an infinite aperture, so when computing photometry for smaller radii, aperture corrections must be applied.
The inputs below can be changed to any desired obsmode, with examples of alternate parameters shown as commented lines.
#Detector:
detectors = ['uvis1']
#detectors = ['uvis1', 'uvis2'] # both UVIS chips
#detectors = ['ir'] # if using IR, must update the filtnames below
#Filters:
filtnames = ['f200lp','f218w','f225w','f275w','f280n','f300x', 'f336w','f343n','f350lp',
'f373n', 'f390m','f390w','f395n','f410m','f438w', 'f467m','f469n','f475w',
'f475x', 'f487n','f502n','f547m','f555w','f600lp','f606w','f621m','f625w',
'f631n', 'f645n','f656n','f657n','f658n','f665n', 'f673n','f680n','f689m',
'f763m', 'f775w','f814w','f845m','f850lp','f953n']
#filtnames = ['f606w']
# For IR filters, must set detectors = ['ir'] above
#filtnames = ['f098m','f105w','f110w','f125w','f126n','f127m','f128n','f130n','f132n','f139m','f140w','f153m','f160w','f164n','f167n']
#Date
mjd = '55008' # WFC3/UVIS reference epoch (26Jun2009)
# mjd = str(Time.now().mjd) # Time right now
#Aperture Radius
aper = '6.0' # 151 pixels (infinity) [default behavior]
#aper = '0.396' # 10 pixels for UVIS
#aper = '0.385' # 3 pixels for IR
3. Basic usage for a single 'obsmode' ¶
The calculation of the zeropoints starts with creating a specific bandpass object. Bandpasses generally consist of at least an instrument name, detector name, and filter name, though other parameters (such as the MJD and aperture radius shown above) are optional. For example:
obsmode = 'wfc3,uvis1,f200lp'
bp = stsyn.band(obsmode) # creates bandpass object
Optional parameters are supplied on the end of the basic bandpass:
obsmode = 'wfc3,uvis1,f200lp,mjd#55008,aper#6.0'
# or to use parameters above:
obsmode = 'wfc3,{},{},mjd#{},aper#{}'.format(detectors[0],filtnames[0],mjd,aper)
bp = stsyn.band(obsmode)
4. Computing zeropoints and other photometric properties ¶
With the bandpass objects we can now calculate zeropoints, pivot wavelengths, and photometric bandwidths. To calculate Vegamag zeropoints, we use the Vega spectrum to calculate the flux in a given bandpass.
def calculate_values(detector, filt, mjd, aper):
# parameters can be removed from obsmode as needed
obsmode = 'wfc3,{},{},mjd#{},aper#{}'.format(detector, filt, mjd, aper)
bp = stsyn.band(obsmode)
# STMag
photflam = bp.unit_response(stsyn.conf.area) # inverse sensitivity in flam
stmag = -21.1 -2.5 * np.log10(photflam.value)
# Pivot Wavelength and bandwidth
photplam = bp.pivot() # pivot wavelength in angstroms
bandwidth = bp.photbw() # bandwidth in angstroms
# ABMag
abmag = stmag - 5 * np.log10(photplam.value) + 18.6921
# Vegamag
obs = Observation(stsyn.Vega, bp, binset=bp.binset) # synthetic observation of vega in bandpass using vega spectrum
vegamag = -obs.effstim(flux_unit='obmag', area=stsyn.conf.area)
return obsmode, photplam.value, bandwidth.value, photflam.value, stmag, abmag, vegamag.value
obsmode, photplam, bandwidth, photflam, stmag, abmag, vegamag = calculate_values(detectors[0], filtnames[0], mjd, aper)
# print values
print('Obsmode PivotWave Photflam STMAG ABMAG VEGAMAG')
print(f'{obsmode}, {photplam:.1f}, {photflam:.4e}, {stmag:.3f}, {abmag:.3f}, {vegamag:.3f}')
5. Iterating over multiple 'obsmodes' ¶
To calculate zeropoints for multiple detectors and/or filters:
oms, pivots, bws, pfs, st, ab, vm = [], [], [], [], [], [], []
print('Obsmode PivotWave Photflam STMAG ABMAG VEGAMAG')
for detector in detectors:
for filt in filtnames:
res = calculate_values(detector, filt, mjd, aper)
obsmode, photplam, bandwidth, photflam, stmag, abmag, vegamag = res # solely for readability
# print values
print(f'{obsmode}, {photplam:.1f}, {photflam:.4e}, {stmag:.3f}, {abmag:.3f}, {vegamag:.3f}')
oms.append(obsmode)
pivots.append(photplam)
bws.append(bandwidth)
pfs.append(photflam)
st.append(stmag)
ab.append(abmag)
vm.append(vegamag)
Values can also be written into an astropy table:
tbl = Table([oms, pivots, bws, pfs, st, ab, vm],
names=['Obsmode', 'Pivot Wave', 'Bandwidth', 'Photflam', 'STMag', 'ABMag', 'VegaMag'])
# Just for rounding columns to smaller number of decimals
for col in tbl.itercols():
if col.name == 'Photflam':
col.info.format = '.4e'
elif col.info.dtype.kind == 'f':
col.info.format = '.3f'
# Show table
tbl
# Write to a file
tbl.write('uvis_zp_tbl.txt', format='ascii.commented_header')
6. Creating and plotting 'total system throughput' tables¶
def calculate_bands(bp, save=False):
# Pass in bandpass object as bp
waves = bp.waveset
throughput = bp(waves)
if save:
tmp = Table([waves, throughput], names=['WAVELENGTH', 'THROUGHPUT'])
tmp.write(','.join(bp.obsmode.modes)+'.txt', format='ascii.commented_header')
return (waves, throughput)
The above function returns a tuple containing two objects, the first being an array of wavelengths, and the second being the throughput at each of those wavelengths. The result can be plotted via:
wl, tp = calculate_bands(bp)
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10,5))
plt.plot(wl, tp)
plt.xlim(1500, 11000)
plt.xlabel('Wavelength [Angstroms]')
plt.ylabel('Throughput')
plt.title('WFC3,UVIS1,F200LP')
To save the curve in an ascii table, simply pass the argument save=True:
calculate_bands(bp, save=True)
ls *txt
To save curves for all obsmodes in the input list:
for det in detectors:
for filt in filtnames:
obsmode = 'wfc3,{},{}'.format(det, filt)
bp = stsyn.band(obsmode)
calculate_bands(bp, save=True)
ls wfc3*txt