In [1]:
# Import the usual libraries
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
#import matplotlib.patches as mpatches
# Enable inline plotting
%matplotlib inline
#from IPython.display import display, Latex, clear_output
In [2]:
# WebbPSF
import webbpsf
# Interpolation and extrapolation
from scipy.interpolate import griddata, RegularGridInterpolator
from scipy.ndimage import rotate
In [3]:
# pySIAF stuff for plotting
import pysiaf
from pysiaf.siaf import Siaf
from pysiaf.siaf import plot_main_apertures
siaf = Siaf('NIRCam')
siaf.generate_toc()
In [4]:
from astropy.table import Table
In [5]:
import pynrc
from pynrc import nrc_utils
from pynrc.nrc_utils import read_filter, gen_psf_coeff, gen_image_coeff
from pynrc.nrc_utils import gen_webbpsf_siwfe, field_coeff_resid, field_coeff_func
from pynrc.nrc_utils import Tel2Sci_info, NIRCam_V2V3_limits, radial_std
from pynrc.nrc_utils import jl_poly_fit, jl_poly
In [6]:
from pynrc.nrc_utils import frebin, dist_image, hist_indices, binned_statistic
from astropy.convolution import convolve, convolve_fft, Gaussian1DKernel, Gaussian2DKernel
def do_contrast(im1, im2, header, nsig=5, supersample=True, func=np.std):
data = im1 - im2
pixscale = header['PIXELSCL']
oversample = header['OVERSAMP']
rr, stds = radial_std(data, pixscale=pixscale, oversample=oversample, supersample=supersample, func=func)
imavg = (im1+im2)/2
if supersample or oversample==1:
psf_max = imavg.max()
else:
psf_max = frebin(imavg, scale=1/oversample).max()
return rr, nsig*stds/psf_max
import matplotlib.ticker as mticker
class MathTextSciFormatter(mticker.Formatter):
def __init__(self, fmt="%.2e"):
self.fmt = fmt
def __call__(self, x, pos=None):
s = self.fmt % x
decimal_point = '.'
positive_sign = '+'
tup = s.split('e')
significand = tup[0].rstrip(decimal_point)
sign = tup[1][0].replace(positive_sign, '')
exponent = tup[1][1:].lstrip('0')
if exponent:
exponent = '10^{%s%s}' % (sign, exponent)
if significand and exponent:
s = r'%s{\times}%s' % (significand, exponent)
else:
s = r'%s%s' % (significand, exponent)
return "${}$".format(s)
In [7]:
# PSF setup
filter = 'F430M'
kwargs = {}
kwargs['pupil'] = None #'CIRCLYOT'
kwargs['force'] = True
kwargs['save'] = False
kwargs['save_name'] = None
bp = read_filter(filter)
channel = 'SW' if bp.avgwave() < 24000 else 'LW'
module = kwargs.get('module', 'A') # If not specified, choose 'A'
kwargs['module'] = module
# Check if coronagraphy
pupil = kwargs.get('pupil', 'CLEAR') # If not specified, choose 'A'
kwargs['pupil'] = pupil
kwargs['detector'] = None
kwargs['detector_position'] = None
kwargs['include_si_wfe'] = True
kwargs['use_legendre'] = True
fov_pix = kwargs['fov_pix'] = 127
oversample = kwargs['oversample'] = 4
kwargs['opd'] = ('OPD_RevW_ote_for_NIRCam_requirements.fits.gz', 0)
kwargs['jitter'] = 'gaussian'
kwargs['jitter_sigma'] = 0.007
apname = 'NRCA5_FULL'
In [8]:
kwargs
Out[8]:
In [9]:
# SI Zernike data
if (pupil is not None) and ('LYOT' in pupil):
zfile = 'si_zernikes_coron_wfe.fits'
else:
zfile = 'si_zernikes_isim_cv3.fits'
data_dir = webbpsf.utils.get_webbpsf_data_path() + '/'
zernike_file = data_dir + zfile
ztable_full = Table.read(zernike_file)
mod = channel + module
ind_nrc = ['NIRCam'+mod in row['instrument'] for row in ztable_full]
ind_nrc = np.where(ind_nrc)
v2_all = np.array(ztable_full[ind_nrc]['V2'].tolist())
v3_all = np.array(ztable_full[ind_nrc]['V3'].tolist())
In [10]:
# WFE Drift Values
wfe_list = [0, 2, 5, 10]
# Do a subset, first 3 and last 3
# v2_all = np.append(v2_all[0:3], v2_all[-3:])
# v3_all = np.append(v3_all[0:3], v3_all[-3:])
v2_calc = v2_all[0:3]
v3_calc = v3_all[0:3]
ap = siaf[apname]
xvals = np.random.randint(0,2048,size=4)
yvals = np.random.randint(0,2048,size=4)
#xvals, yvals = np.array([800,1000]), np.array([1200,1000])
v2_calc, v3_calc = np.array(ap.sci_to_tel(xvals, yvals)) / 60
In [11]:
%%time
nrc = pynrc.NIRCam(filter=filter, use_legendre=True, include_si_wfe=True, apname=apname,
fov_pix=fov_pix, oversample=oversample)
In [29]:
%%time
nrc.wfe_drift = True
In [30]:
%%time
nrc.wfe_field = True
# cf_fields_grid = nrc._psf_coeff_mod['si_field']
# v2grid = nrc._psf_coeff_mod['si_field_v2grid']
# v3grid = nrc._psf_coeff_mod['si_field_v3grid']
In [15]:
%%time
nrc = pynrc.NIRCam(filter=filter, use_legendre=True, include_si_wfe=True, apname=apname,
fov_pix=fov_pix, oversample=oversample)
In [16]:
cf = nrc.psf_coeff
cf_hdr = nrc._psf_coeff_hdr
In [17]:
# varr = [(518,510), (500,500), (520,530), (520,505)]
varr = [(nx/2-4,ny/2+4), (nx/2-14,ny/2-14), (nx/2+16,ny/2+6), (nx/2-9,ny/2+6)]
for (xind, yind) in varr:
xind, yind = int(xind), int(yind)
vals = np.abs(cf[:,yind, xind])
label='Pixel ({:.0f}, {:.0f})'.format((xind-nx/2), (yind-ny/2))
plt.semilogy(vals / np.sum(vals), label=label)
plt.legend()
plt.title('Power spectrum')
plt.xlabel('Coefficients')
Out[17]:
In [106]:
from matplotlib.colors import SymLogNorm, LogNorm
fig, axes = plt.subplots(2,3, figsize=(6,4), sharex=True, sharey=True)
axes = axes.flatten()
extent = 0.5 * np.array([-1,1,-1,1]) * fov_pix * nrc.pixelscale
xylim = np.array([-1,1]) * 1.2
for i in range(len(axes)):
ax = axes[i]
ax.imshow(cf[i], extent=extent)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i>2: ax.set_xlabel('Arcsec')
if i % 3 == 0: ax.set_ylabel('Arcsec')
ax.tick_params(axis='both', color='white', which='both')
for k in ax.spines.keys():
ax.spines[k].set_color('white')
ax.xaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
ax.yaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
fig.tight_layout()
fig.subplots_adjust(wspace=0.01, hspace=0.01, left=0, right=1, top=1, bottom=0)
In [84]:
inst = webbpsf.NIRCam()
inst.filter = filter
inst.SHORT_WAVELENGTH_MIN = inst.LONG_WAVELENGTH_MIN = 1e-7
inst.SHORT_WAVELENGTH_MAX = inst.LONG_WAVELENGTH_MAX = 10e-6
w1, w2 = (2.4,5.1)
ndeg = 9
dn = 20
npsf = np.ceil(dn * (w2-w1))
npsf = 5 if npsf<5 else int(npsf)
npsf = ndeg+1 if npsf<=ndeg else int(npsf)
waves = np.linspace(w1, w2, npsf)
In [ ]:
pynrc.setup_logging('WARN')
worker_arguments = [(inst, wlen, fov_pix, oversample) for wlen in waves]
hdu_arr = []
for wa in worker_arguments:
hdu = nrc_utils._wrap_coeff_for_mp(wa)
if hdu is None:
raise RuntimeError('Returned None values. Issue with multiprocess or WebbPSF??')
hdu_arr.append(hdu)
pynrc.setup_logging('INFO')
In [111]:
images = []
for hdu in hdu_arr:
images.append(hdu.data)
images = np.array(images)
# Simultaneous polynomial fits to all pixels using linear least squares
coeff_all = jl_poly_fit(waves, images, deg=ndeg, use_legendre=True, lxmap=[w1,w2])
In [129]:
fig, ax = plt.subplots(1,1, figsize=(8,5))
# varr = [(518,510), (500,500), (520,530), (520,505)]
varr = [(nx/2-4,ny/2+4), (nx/2-14,ny/2-14), (nx/2+16,ny/2+6), (nx/2-9,ny/2+6)]
colarr = ['C0', 'C1', 'C2', 'C3']
for i, (xind, yind) in enumerate(varr):
xind, yind = int(xind), int(yind)
vals = images[:,yind,xind]
label='Pixel ({:.0f}, {:.0f})'.format((xind-nx/2), (yind-ny/2))
ax.plot(waves, vals, label=label, marker='.', ls='none', color=colarr[i])
vals = jl_poly(waves, coeff_all[:,yind,xind], use_legendre=True, lxmap=[w1,w2])
ax.plot(waves, vals, label=label+' Fit', color=colarr[i])
ax.legend()
ax.set_xlabel('Wavelength')
ax.set_ylabel('Flux Value')
ax.set_title('Example Pixel Fits')
Out[129]:
In [36]:
fig, ax = plt.subplots(1,1)
nrc.siaf_ap.plot(fill=False, ax=ax)
ax.plot(v2_all*60, v3_all*60, marker='x', ls='none')
ax.set_title('CV3 WFE Measurements - Module A')
Out[36]:
In [11]:
pynrc.setup_logging('WARN')
In [12]:
%%time
nrc = pynrc.NIRCam(filter=filter, use_legendre=True, include_si_wfe=True, apname=apname,
fov_pix=fov_pix, oversample=oversample)
In [13]:
%%time
nrc.wfe_drift = True
nrc.wfe_field = True
In [14]:
from pynrc.nrc_utils import field_coeff_func
cf_fit = nrc._psf_coeff_mod['si_field']
v2grid = nrc._psf_coeff_mod['si_field_v2grid']
v3grid = nrc._psf_coeff_mod['si_field_v3grid']
res = field_coeff_func(v2grid, v3grid, cf_fit, 2, -8)
In [15]:
cf_fit.shape
Out[15]:
In [19]:
fig, axes = plt.subplots(2,2, figsize=(6.2,6.5))
axes = axes.flatten()
ny, nx = cf_fit.shape[3:]
varr = [(nx/2-4,ny/2+4), (nx/2-14,ny/2-14), (nx/2+16,ny/2+6), (nx/2-9,ny/2+6)]
extent = np.array([v2grid.min(), v2grid.max(), v3grid.min(), v3grid.max()]) * 60
for i, (xind,yind) in enumerate(varr):
ax = axes[i]
xind, yind = int(xind), int(yind)
ax.imshow(cf_fit[:,:,1,yind, xind], extent=extent, interpolation='bicubic')
nrc.siaf_ap.plot(fill=False, ax=ax, color='C1')
ax.plot(v2_all*60, v3_all*60, marker='x', ls='none', color='C3')
ax.set_title('PSF Pixel ({:.0f}, {:.0f})'.format(xind-nx/2,yind-ny/2))
fig.suptitle('Coefficient Residuals for $P_{1}(x)$', fontsize=16)
fig.tight_layout()
fig.subplots_adjust(top=0.9)
In [20]:
v23_0 = np.array(ap.sci_to_tel(1024, 1024)) / 60
%time hdul0 = gen_webbpsf_siwfe(filter, v23_0, wfe_drift=0, **kwargs)
In [177]:
psf0_dict = {}
for wfe in wfe_list:
hdul = gen_webbpsf_siwfe(filter, v23_0, wfe_drift=wfe, **kwargs)
k_wfe = 'wfe_{:.0f}nm'.format(wfe)
psf0_dict[k_wfe] = hdul
In [21]:
pynrc.setup_logging('WARN', verbose=False)
psf_dict = {}
for (v2, v3) in zip(v2_calc, v3_calc):
print(v2,v3)
d = {}
for wfe in wfe_list:
hdul = gen_webbpsf_siwfe(filter, (v2,v3), wfe_drift=wfe, **kwargs)
k_wfe = 'wfe_{:.0f}nm'.format(wfe)
d[k_wfe] = hdul
k_si = 'V2V3_({:.2f},{:.2f})'.format(v2,v3)
psf_dict[k_si] = d
In [13]:
# i1, i2 = (0,1) # PSF postions
# j1, j2 = (0,0) # WFE drifts
# keys1 = list(psf_dict.keys())
# d1 = psf_dict[keys1[i1]]
# d2 = psf_dict[keys1[i2]]
# keys2 = list(d1.keys())
# hdul1 = d1[keys2[j1]]
# hdul2 = d2[keys2[j2]]
# print(keys1[i1], keys2[j1])
# print(keys1[i2], keys2[j2])
# print(hdul1[1].data.sum(), hdul2[1].data.sum())
In [22]:
%%time
psf0, psf0_over = nrc.gen_psf(wfe_drift=0, coord_vals=(1024,1024), coord_frame='sci',
return_oversample=True)
In [205]:
%%time
psf_list = []
psf_over_list = []
for (v2, v3) in zip(v2_calc, v3_calc):
psf, psf_over = nrc.gen_psf(wfe_drift=0, coord_vals=(v2,v3), coord_frame='tel',
return_oversample=True)
psf_list.append(psf)
psf_over_list.append(psf_over)
In [26]:
from matplotlib.colors import SymLogNorm, LogNorm
fig, axes = plt.subplots(2,4,figsize=(14,6))
im_list = [psf_dict[k]['wfe_0nm'][0].data for k in psf_dict.keys()]
extent = 0.5 * np.array([-1,1,-1,1]) * fov_pix * nrc.pixelscale
for i, im in enumerate(im_list):
ax = axes[0,i]
pos = ax.imshow(im, extent=extent, norm=LogNorm(vmin=1e-7))
ax.set_title('WebbPSF: pix=({:.0f},{:.0f})'.format(xvals[i],yvals[i]))
for i, im in enumerate(psf_over_list):
ax = axes[1,i]
pos = ax.imshow(im, extent=extent, norm=LogNorm(vmin=1e-7))
ax.set_title('pyNRC: pix=({:.0f},{:.0f})'.format(xvals[i],yvals[i]))
xylim = np.array([-1,1]) * 1.5
for i, ax in enumerate(axes.flatten()):
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i>3: ax.set_xlabel('Arcsec')
if i % 4 == 0: ax.set_ylabel('Arcsec')
ax.tick_params(axis='both', color='white', which='both')
for k in ax.spines.keys():
ax.spines[k].set_color('white')
ax.xaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
ax.yaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
fig.suptitle('PSF Comparisons ({})'.format(nrc.filter), fontsize=15)
fig.tight_layout()
fig.subplots_adjust(top=0.88)
cb = fig.colorbar(pos, ax=axes.ravel().tolist())
cb.ax.yaxis.set_major_formatter(MathTextSciFormatter("%.1e"))
In [27]:
from matplotlib.colors import SymLogNorm, LogNorm
fig, axes = plt.subplots(2,4,figsize=(14,6))
im_list = [psf_dict[k]['wfe_0nm'][0].data for k in psf_dict.keys()]
extent = 0.5 * np.array([-1,1,-1,1]) * fov_pix * nrc.pixelscale
vmin, vmax = np.array([-1,1])*1e-5
vlevels_list = []
for i, im in enumerate(im_list):
ax = axes[0,i]
diff = im - hdul0[0].data
# vmin, vmax = np.array([-1,1]) * np.mean(np.abs([diff.min(),diff.max()])) / 2
vlevels_list.append((vmin,vmax))
pos = ax.imshow(diff, extent=extent, vmin=vmin, vmax=vmax, cmap='RdBu_r')
ax.set_title('WebbPSF: pix=({:.0f},{:.0f})'.format(xvals[i],yvals[i]))
# fig.colorbar(pos, ax=ax, format='%.1e')
for i, im in enumerate(psf_over_list):
ax = axes[1,i]
diff = im - psf0_over
# vmin, vmax = vlevels_list[i]
pos = ax.imshow(diff, extent=extent, vmin=vmin, vmax=vmax, cmap='RdBu_r')
ax.set_title('pyNRC: pix=({:.0f},{:.0f})'.format(xvals[i],yvals[i]))
# fig.colorbar(pos, ax=ax, format='%.1e')
xylim = np.array([-1,1]) * 1.5
for i, ax in enumerate(axes.flatten()):
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i>3: ax.set_xlabel('Arcsec')
if i % 4 == 0: ax.set_ylabel('Arcsec')
ax.tick_params(axis='both', color='white', which='both')
for k in ax.spines.keys():
ax.spines[k].set_color('white')
ax.xaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
ax.yaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
fig.suptitle('{} Differenced Images (ref pixel = (1024,1024))'.format(nrc.filter), fontsize=15)
fig.tight_layout()
fig.subplots_adjust(top=0.88)
cb = fig.colorbar(pos, ax=axes.ravel().tolist())
cb.ax.yaxis.set_major_formatter(MathTextSciFormatter("%.1e"))
In [207]:
fig, axes = plt.subplots(1,4,figsize=(14,4))
im_list = [psf_dict[k]['wfe_0nm'][0].data for k in psf_dict.keys()]
header = hdul0[0].header
for i, im in enumerate(im_list):
ax = axes[i]
diff = im - hdul0[0].data
rr, cont = do_contrast(im, hdul0[0].data, header, nsig=5, supersample=False, func=robust.std)
ax.semilogy(rr, cont, label='WebbPSF')
ax.set_title('pix=({:.0f},{:.0f})'.format(xvals[i],yvals[i]))
for i, im in enumerate(psf_over_list):
ax = axes[i]
diff = im - hdul0[0].data
rr, cont = do_contrast(im, psf0_over, header, nsig=5, supersample=False, func=robust.std)
ax.semilogy(rr, cont, label='pyNRC')
axes[0].set_ylabel('Contrast')
axes[0].legend()
for ax in axes:
ax.set_ylim([1e-6,2e-2])
ax.set_xlabel('Arcsec')
fig.suptitle('Contrast Curves (ref pixel = (1024,1024))', fontsize=15)
fig.tight_layout()
fig.subplots_adjust(top=0.84)
In [50]:
kwargs
Out[50]:
In [55]:
# Cycle through WFE drifts for fitting
wfe_list = np.array([0,1,2,5,10,20,30,40])
nwfe = len(wfe_list)
cf_wfe = []
for wfe in wfe_list:
print('WFE Drift: {} nm'.format(wfe))
cf, _ = gen_psf_coeff(nrc.bandpass, wfe_drift=wfe, **kwargs)
cf_wfe.append(cf)
cf_wfe = np.array(cf_wfe)
# Get residuals
cf_wfe = cf_wfe - cf_wfe[0]
In [30]:
cf_resid = nrc._psf_coeff_mod['wfe_drift']
lxmap = nrc._psf_coeff_mod['wfe_drift_lxmap']
In [96]:
fig, axes = plt.subplots(3,3, figsize=(6,6))
indy, indx = int(ny/2+2),int(nx/2-1)
for i, ax in enumerate(axes.flatten()):
cf = cf_resid[:,i, indy,indx]
wfe_arr = np.array([0,1,2,5,10,20,40])
cf_vals = jl_poly(wfe_arr, cf, use_legendre=True, lxmap=lxmap)
# ax.plot(wfe_arr, cf_vals, ls='none', marker='o')
ax.plot(wfe_list, cf_wfe[:,i,indy,indx], ls='none', marker='o')
wfe_arr = np.arange(0,41)
cf_vals = jl_poly(wfe_arr, cf, use_legendre=True, lxmap=lxmap)
ax.plot(wfe_arr, cf_vals)
labels = [item.get_text() for item in ax.get_yticklabels()]
empty_string_labels = ['']*len(labels)
ax.set_yticklabels(empty_string_labels)
labels = [item.get_text() for item in ax.get_xticklabels()]
empty_string_labels = ['']*len(labels)
if i<6: ax.set_xticklabels(empty_string_labels)
fig.text(0.5, 0., 'WFE Drift (nm RMS)', ha='center', va='bottom', fontsize=14)
fig.text(0., 0.5, 'Relative Coefficient Modifications', va='center', ha='left', fontsize=14, rotation=90)
txt = 'Pixel ({:.0f}, {:.0f}) for LWA (coeffs 0-8)'.format(indx-nx/2,indy-nx/2)
fig.suptitle(txt, fontsize=16)
fig.tight_layout()
fig.subplots_adjust(bottom=0.08, left=0.05, top=0.92)
In [115]:
from matplotlib.colors import SymLogNorm, LogNorm
fig, axes = plt.subplots(2,3, figsize=(6,4), sharex=True, sharey=True)
axes = axes.flatten()
extent = 0.5 * np.array([-1,1,-1,1]) * fov_pix * nrc.pixelscale
xylim = np.array([-1,1]) * 1.2
for i in range(len(axes)):
ax = axes[i]
ax.imshow(cf_resid[1,i,:,:], extent=extent)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i>2: ax.set_xlabel('Arcsec')
if i % 3 == 0: ax.set_ylabel('Arcsec')
ax.tick_params(axis='both', color='white', which='both')
for k in ax.spines.keys():
ax.spines[k].set_color('white')
ax.xaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
ax.yaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
fig.tight_layout()
fig.subplots_adjust(wspace=0.01, hspace=0.01, left=0, right=1, top=1, bottom=0)
In [142]:
keys_v2v3 = list(psf_dict.keys())
keys_wfed = list(psf_dict[keys_v2v3[0]].keys())
print(keys_v2v3)
print(keys_wfed)
In [139]:
wfe_list = [0, 2, 5, 10]
In [210]:
%%time
psf_list = []
psf_over_list = []
for wfed in wfe_list:
psf, psf_over = nrc.gen_psf(wfe_drift=wfed, coord_vals=(1024,1024), coord_frame='sci',
return_oversample=True)
psf_list.append(psf)
psf_over_list.append(psf_over)
In [200]:
from matplotlib.colors import SymLogNorm, LogNorm
fig, axes = plt.subplots(2,3,figsize=(11,6))
im_list = [psf0_dict[k][0].data for k in keys_wfed]
extent = 0.5 * np.array([-1,1,-1,1]) * fov_pix * nrc.pixelscale
vmin, vmax = np.array([-1,1])*2e-6
vlevels_list = []
for i, im in enumerate(im_list):
if i==0:
pass
ax = axes[0,i-1]
diff = im - im_list[0]
vlevels_list.append((vmin,vmax))
pos = ax.imshow(diff, extent=extent, vmin=vmin, vmax=vmax, cmap='RdBu_r')
ax.set_title('WebbPSF: $\Delta$WFE={}nm'.format(wfe_list[i]))
for i, im in enumerate(psf_over_list):
if i==0:
pass
ax = axes[1,i-1]
diff = im - psf_over_list[0]
pos = ax.imshow(diff, extent=extent, vmin=vmin, vmax=vmax, cmap='RdBu_r')
ax.set_title('pyNRC: $\Delta$WFE={}nm'.format(wfe_list[i]))
xylim = np.array([-1,1]) * 1.5
for i, ax in enumerate(axes.flatten()):
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i>2: ax.set_xlabel('Arcsec')
if i % 3 == 0: ax.set_ylabel('Arcsec')
ax.tick_params(axis='both', color='white', which='both')
for k in ax.spines.keys():
ax.spines[k].set_color('white')
ax.xaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
ax.yaxis.get_major_locator().set_params(nbins=9, steps=[1, 2, 5, 10])
fig.suptitle('{} Differenced Images -- Det Pixel Loc (1024,1024)'.format(nrc.filter), fontsize=15)
fig.tight_layout()
fig.subplots_adjust(top=0.88)
cb = fig.colorbar(pos, ax=axes.ravel().tolist())
cb.ax.yaxis.set_major_formatter(MathTextSciFormatter("%.1e"))
In [217]:
from pynrc import robust
fig, ax = plt.subplots(1,1, figsize=(8,5))
lin_vals = np.linspace(0.2,0.8,len(wfe_list)-1)
c1 = plt.cm.Blues_r(lin_vals)[::-1]
c2 = plt.cm.Reds_r(lin_vals)[::-1]
im_list = [psf0_dict[k][0].data for k in keys_wfed]
wfe_rev = wfe_list[1:][::-1]
im_list_rev = im_list[1:][::-1]
for i, im in enumerate(im_list_rev):
rr, cont = do_contrast(im, im_list[0], header, nsig=5, supersample=False, func=robust.std)
ax.semilogy(rr, cont, color=c1[i], label='WebbPSF: {} nm'.format(wfe_rev[i]))
im_list_rev = psf_over_list[1:][::-1]
for i, im in enumerate(im_list_rev):
rr, cont = do_contrast(im, psf_over_list[0], header, nsig=5, supersample=False, func=robust.std)
ax.semilogy(rr, cont, color=c2[i], label='pyNRC: {} nm'.format(wfe_rev[i]))
ax.legend(ncol=2)
ax.set_xlabel('Arcsec')
ax.set_ylabel('5-$\sigma$ Contrast')
fig.tight_layout()
In [ ]:
# Plot field contrast curves as a function of distance.
# Plot contrast at 1"
In [225]:
sep_asec = np.array([1,2,5,10,20,30,60])
v2_arr = np.ones([4,len(sep_asec)]) * v23_0[0]
v3_arr = np.ones([4,len(sep_asec)]) * v23_0[1]
v2_arr[0] -= sep_asec/60
v2_arr[1] += sep_asec/60
v3_arr[2] -= sep_asec/60
v3_arr[3] += sep_asec/60
In [229]:
psf0, psf0_over = nrc.gen_psf(wfe_drift=0, coord_vals=(v23_0[0],v23_0[1]), coord_frame='tel',
return_oversample=True)
In [245]:
# Field Dependent PSFs
psf_over_list = []
for (v2, v3) in zip(v2_arr.flatten(), v3_arr.flatten()):
_, psf_over = nrc.gen_psf(wfe_drift=0, coord_vals=(v2,v3), coord_frame='tel',
return_oversample=True)
psf_over_list.append(psf_over)
In [258]:
# WFE Drift PSFs
wfe_list = [1,2,5,10]
psf_over_list2 = []
for wfed in wfe_list:
_, psf_over = nrc.gen_psf(wfe_drift=wfed, coord_vals=(1024,1024), coord_frame='sci',
return_oversample=True)
psf_over_list2.append(psf_over)
In [282]:
# Field Contrasts
cont_arr = []
for im in psf_over_list:
rr, cont = do_contrast(im, psf0_over, header, nsig=5, supersample=False, func=np.std)
cont_arr.append(cont)
cont_arr = np.array(cont_arr).reshape([4,len(sep_asec),-1])
cont_med = np.mean(cont_arr, axis=0)
# WFE Drift Contrasts
cont_arr2 = []
for im in psf_over_list2:
rr, cont = do_contrast(im, psf0_over, header, nsig=5, supersample=False, func=np.std)
cont_arr2.append(cont)
cont_arr2 = np.array(cont_arr2)
In [283]:
print(cont_arr.shape, cont_med.shape, cont_arr2.shape)
In [307]:
fig, axes = plt.subplots(2,1, figsize=(7.5,7))
lin_vals = np.linspace(0.2,0.8,cont_med.shape[0])
c1 = plt.cm.Blues_r(lin_vals)
lin_vals = np.linspace(0.1,0.6,len(wfe_list))
c2 = plt.cm.Reds_r(lin_vals)
ax = axes[0]
for i, cont in enumerate(cont_med):
ax.semilogy(rr, cont, color=c1[i], label='Dist = {}"'.format(sep_asec[i]))
for i, cont in enumerate(cont_arr2):
ax.semilogy(rr, cont, color=c2[i], label='$\Delta$WFE = {} nm'.format(wfe_list[i]), ls='--')
ax.legend(ncol=2)
ax.set_xlabel('Radial Distance (Arcsec)')
ax.set_ylabel('5-$\sigma$ Contrast')
ax.set_title('Contrast Curves')
ax = axes[1]
vals = np.array([np.interp(1, rr, cont) for cont in cont_med])
ax.semilogy(sep_asec, vals, lw=5)
xlim = ax.get_xlim()
for i, cont in enumerate(cont_arr2[::-1]):
xvals = xlim
yvals = np.interp(1,rr,cont) * np.array([1,1])
label='$\Delta$WFE = {} nm'.format(wfe_list[::-1][i])
ax.plot(xvals, yvals, color=c2[::-1][i], label=label, ls='--')
ax.set_xlim(xlim)
ax.set_xlabel('PSF Separation (Arcsec)')
ax.set_ylabel('5-$\sigma$ Contrast at 1"')
ax.set_title('Contrast at 1"')
ax.legend(loc='lower right')
fig.tight_layout()
fig.savefig('contrast_compare.pdf')
In [16]:
xvals = np.random.randint(0,2048,size=1)
yvals = np.random.randint(0,2048,size=1)
v2, v3 = np.array(ap.sci_to_tel(xvals, yvals)) / 60
pts = np.array([v2,v3]).transpose()
cf_fit = nrc._psf_coeff_mod['si_field']
v2grid = nrc._psf_coeff_mod['si_field_v2grid']
v3grid = nrc._psf_coeff_mod['si_field_v3grid']
%time cf_resid = field_coeff_func(v2grid, v3grid, cf_fit, v2, v3)
In [17]:
%time cf_new = nrc._psf_coeff + cf_resid
psf_coeff_hdr = nrc._psf_coeff_hdr
psf_info = nrc._psf_info
In [41]:
bp = nrc.bandpass
waveset = np.copy(bp.wave)
binsize = 3
excess = waveset.size % binsize
waveset = waveset[:waveset.size-excess]
waveset = waveset.reshape(-1,binsize) # Reshape
waveset = waveset[:,binsize//2] # Use the middle values
waveset = np.concatenate(([bp.wave[0]],waveset,[bp.wave[-1]]))
wgood = waveset / 1e4
w1 = wgood.min()
w2 = wgood.max()
wrange = w2 - w1
use_legendre = True if psf_coeff_hdr['LEGNDR'] else False
lxmap = [psf_coeff_hdr['WAVE1'], psf_coeff_hdr['WAVE2']]
%time psf_fit = nrc_utils.jl_poly(wgood, cf_new, dim_reorder=True, use_legendre=use_legendre, lxmap=lxmap)
#%time psf_fit[psf_fit<=0] = np.min(psf_fit[psf_fit>0]) / 10
In [28]:
psf_fit.shape
Out[28]:
In [22]:
from pynrc.nrc_utils import S
In [33]:
sp_flat = S.ArraySpectrum(waveset, 0*waveset + 10.)
sp_norm = sp_flat.renorm(bp.unit_response(), 'flam', bp)
sp_norm = [sp_norm]
obs_list = [S.Observation(sp, bp, binset=waveset) for sp in sp_norm]
In [34]:
obs_list[0].binflux.shape
Out[34]:
In [31]:
psf_temp = psf_fit.reshape([356,-1])
In [36]:
%time test = psf_temp * obs_list[0].binflux.reshape([-1,1])
In [40]:
psf_list = []
psf_temp = psf_fit.reshape([356,-1])
for obs in obs_list:
psf_list.appen(psf_temp)
In [42]:
%time psf_list = [psf_fit*obs.binflux for obs in obs_list]
In [80]:
test = psf_fit[:,:,3]
len(test[test<0])
Out[80]:
In [83]:
%time fin = nrc_utils.krebin(data_over, (257,257))
In [81]:
data_over = psf_fit.sum(axis=2)
data_over[data_over<=0]
Out[81]:
In [69]:
data_over.shape
Out[69]:
In [65]:
len(ind)
Out[65]:
In [45]:
%%time
psf = gen_image_coeff(nrc.bandpass,
coeff=cf_new, coeff_hdr=psf_coeff_hdr,
fov_pix=psf_info['fov_pix'], oversample=psf_info['oversample'],
return_oversample=False)
In [194]:
coeff, head = nrc_utils.gen_psf_coeff('F444W', fov_pix=fov_pix, oversample=oversample, quick=True, save=False, force=True,
use_legendre=True, ndeg=11)
In [199]:
print(head['NPSF'], head['NDEG'])
lxmap = [head['WAVE1'], head['WAVE2']]
test = nrc_utils.jl_poly(np.array([4.4]), coeff, dim_reorder=True, use_legendre=True, lxmap=lxmap)
test[test<=0] = np.min(test[test>0])
In [200]:
print(test.min(), test.sum())
In [201]:
plt.imshow(test[:,:,0], norm=LogNorm(vmin=1e-8))
Out[201]:
In [47]:
hdr = hdul1[0].header
rr, cont = do_contrast(hdul1[0].data, hdul2[0].data, hdr, supersample=False)
plt.semilogy(rr, cont)
rr, cont = do_contrast(psf1_over, psf2_over, hdr, supersample=False)
plt.semilogy(rr, cont)
#plt.ylim([1e-4,1e-2])
Out[47]:
In [49]:
hdr = hdul1[0].header
rr, cont = do_contrast(hdul1[0].data, hdul2[0].data, hdr, supersample=False)
plt.semilogy(rr, cont)
rr, cont = do_contrast(psf1_over, psf2_over, hdr, supersample=False)
plt.semilogy(rr, cont)
#plt.ylim([1e-4,1e-2])
Out[49]:
In [8]:
# SI Zernike data
data_dir = webbpsf.utils.get_webbpsf_data_path() + '/'
zernike_file = data_dir + 'si_zernikes_isim_cv3.fits'
# Read in table data
ztable_full = Table.read(zernike_file)
In [6]:
# PSF setup
filter = 'F405N'
kwargs = {}
kwargs['pupil'] = None #'CIRCLYOT'
kwargs['force'] = False
kwargs['save'] = True
kwargs['save_name'] = None
bp = read_filter(filter)
channel = 'SW' if bp.avgwave() < 24000 else 'LW'
module = kwargs.get('module', 'A') # If not specified, choose 'A'
kwargs['module'] = module
# Check if coronagraphy
pupil = kwargs.get('pupil', 'CLEAR') # If not specified, choose 'A'
kwargs['pupil'] = pupil
kwargs['detector'] = None
kwargs['detector_position'] = None
kwargs['include_si_wfe'] = True
kwargs['use_legendre'] = False
fov_pix = kwargs['fov_pix'] = 127
oversample = kwargs['oversample'] = 4
kwargs['opd'] = ('OPD_RevW_ote_for_NIRCam_requirements.fits.gz', 0)
kwargs['jitter'] = 'gaussian'
kwargs['jitter_sigma'] = 0.007
In [7]:
# SI Zernike data
if (pupil is not None) and ('LYOT' in pupil):
zfile = 'si_zernikes_coron_wfe.fits'
else:
zfile = 'si_zernikes_isim_cv3.fits'
data_dir = webbpsf.utils.get_webbpsf_data_path() + '/'
zernike_file = data_dir + zfile
ztable_full = Table.read(zernike_file)
mod = channel + module
ind_nrc = ['NIRCam'+mod in row['instrument'] for row in ztable_full]
ind_nrc = np.where(ind_nrc)
v2_all = np.array(ztable_full[ind_nrc]['V2'].tolist())
v3_all = np.array(ztable_full[ind_nrc]['V3'].tolist())
# Add detector corners
v2_min, v2_max, v3_min, v3_max = NIRCam_V2V3_limits(module, channel=channel, pupil=pupil, rederive=True, border=1)
igood = v3_all > v3_min
v2_all = np.append(v2_all[igood], [v2_min, v2_max, v2_min, v2_max])
v3_all = np.append(v3_all[igood], [v3_min, v3_min, v3_max, v3_max])
In [8]:
coeff0, cf_hdr = gen_psf_coeff(filter, **kwargs)
In [9]:
kwargs
Out[9]:
In [10]:
from pynrc.nrc_utils import field_coeff_resid
cf_fields, v2_all, v3_all = field_coeff_resid(filter, coeff0, return_raw=True, **kwargs)
# v2_corn, v3_corn = NIRCam_V2V3_limits(module, channel=channel, pupil=pupil, rederive=True,
# return_corners=True, border=0)
In [9]:
pynrc.setup_logging('INFO')
In [11]:
plt.plot(v2_all, v3_all, marker='o', ls='none')
v2_corn, v3_corn = NIRCam_V2V3_limits(module, channel=channel, pupil=pupil, rederive=True, return_corners=True, border=1)
plt.plot(v2_corn, v3_corn, marker='o', ls='none')
Out[11]:
In [12]:
# Let's work with a single pixel and its coefficients to look at how it changes w.r.t. position
yind = xind = int(fov_pix*oversample / 2)
# Create fine mesh grid
v2_min, v2_max = v2_all.min(), v2_all.max()
v3_min, v3_max = v3_all.min(), v3_all.max()
dstep = 1. / 60. # 1" steps
xgrid = np.arange(v2_min, v2_max+dstep, dstep)
ygrid = np.arange(v3_min, v3_max+dstep, dstep)
V2_interp, V3_interp = np.meshgrid(xgrid,ygrid)
#V2_interp = V2_interp.ravel()
#V3_interp = V3_interp.ravel()
cf_vals_one = cf_fields[:,:,0:5,0:5]
print(cf_fields.shape, cf_vals_one.shape, v2_all.shape, V2_interp.shape, cf_vals_one.shape)
In [162]:
# ap = siaf['NRCA5_FULL']
# nvals = 1
# xvals = np.random.randint(0,2048,size=nvals)
# yvals = np.random.randint(0,2048,size=nvals)
# v2_new, v3_new = np.array(ap.sci_to_tel(xvals, yvals)) / 60
# isort = np.argsort(v2_new)
# v2_new = v2_new[isort]
# v3_new = v3_new[isort]
# cf_temp = cf_fields[:,0,:,:]#.reshape([17,-1,508,508])
# %time cf_grid = field_coeff_interp(v2_all, v3_all, cf_temp, v2_new, v3_new)
# print(cf_grid.shape)
In [36]:
cfnorm = np.abs(coeff0.reshape([11,-1]))
cfnorm = cfnorm / cfnorm.sum(axis=0)
cfstd = np.std(cfnorm, axis=1)
cfavg = np.mean(cfnorm, axis=1)
plt.semilogy(cfavg)
plt.plot(cfavg+cfstd)
plt.plot(cfavg-cfstd)
In [50]:
# Create an evenly spaced grid of V2/V3 coordinates
nv23 = 8
v2grid = np.linspace(v2_min, v2_max, num=nv23)
v3grid = np.linspace(v3_min, v3_max, num=nv23)
v2_new, v3_new = np.meshgrid(v2grid,v3grid)
sh = cf_fields.shape
cf_fields_grid = np.zeros([nv23,nv23,sh[1],sh[2],sh[3]])
# Cycle through each coefficient to interpolate onto V2/V3 grid
for i in range(sh[1]):
print(i)
cf_fields_grid[:,:,i,:,:] = griddata((v2_all, v3_all), cf_fields[:,i,:,:], (v2_new, v3_new), method='cubic')
In [55]:
print(cf_fields_grid.shape)
func = RegularGridInterpolator((v3grid, v2grid), cf_fields_grid[:,:,0,yind,xind], method='linear', bounds_error=False, fill_value=None)
pts = np.array([V3_interp, V2_interp]).transpose()
im_interp = func(pts).squeeze()
In [52]:
extent = np.array([V2_interp.min(), V2_interp.max(), V3_interp.min(), V3_interp.max()])
plt.imshow(im_interp, extent=extent)
Out[52]:
In [40]:
extent = np.array([V2_interp.min(), V2_interp.max(), V3_interp.min(), V3_interp.max()])
plt.imshow(im_interp, extent=extent)
Out[40]:
In [ ]:
In [53]:
extent = [v2grid.min(),v2grid.max(),v3grid.min(),v3grid.max()]
plt.imshow(cf_fields_grid[:,:,0,yind,xind], extent=extent)
plt.plot(v2_all, v3_all, marker='o', ls='none', color='C1')
Out[53]:
In [54]:
fig, axes = plt.subplots(1,2, figsize=(12,5))
axes[0].tricontourf(v2_all, v3_all, cf_fields[:,0,yind,xind], levels=100)
axes[1].contourf(v2grid, v3grid, cf_fields_grid[:,:,0,yind,xind], levels=100)
fig.tight_layout()
In [31]:
%time func = RegularGridInterpolator((v3grid, v2grid), cf_fields_grid, method='linear', bounds_error=False, fill_value=None)
In [165]:
pts = np.array([V3_interp, V2_interp]).transpose()
pts.shape
Out[165]:
In [184]:
v2 = np.array([0.5,1,1.5])
v3 = np.array([-9, -8, -7])
pts = np.array([v2,v3]).transpose()
In [185]:
%time test = func(pts).squeeze()
In [187]:
test2 = coeff0 + test
In [194]:
np.size([[5,2],[4,2]])
Out[194]:
In [186]:
print(pts.shape,test.shape)
In [122]:
extent = np.array([V2_interp.min(), V2_interp.max(), V3_interp.min(), V3_interp.max()])
plt.imshow(test[:,:,0], extent=extent)
Out[122]:
In [55]:
from scipy.interpolate import RectBivariateSpline, interp2d
In [74]:
# func = RectBivariateSpline(v3grid, v2grid, cf_grid[:,:,0,yind,xind])
func = interp2d(v2grid, v3grid, cf_grid[:,:,0,yind,xind], kind='cubic')
In [75]:
%time test = func(V2_interp[0,:].ravel(),V3_interp[:,0].ravel())
In [50]:
%time test = func(V3_interp[:,0].ravel(), V2_interp[0,:].ravel())
In [64]:
xt = [0,1,2]; yt = [0,3]; zt = np.array([[1,2,3], [4,5,6]])
zt.s
Out[64]:
In [76]:
plt.imshow(test, extent=extent)
Out[76]:
In [64]:
cf_temp = cf_fields.reshape([17,11,-1])
%time cf_grid = griddata((v2_all, v3_all), cf_temp, (V2_interp[0:100],V3_interp[0:100]), method='cubic')
print(cf_grid.shape)
In [13]:
extent = np.array([xgrid.min(),xgrid.max(),ygrid.min(),ygrid.max()])
plt.imshow(cf_grid[:,:,0], extent=extent)
Out[13]:
In [ ]:
# Work
In [14]:
# least squares estimation
# X*A = Z
X = v2_all
Y = v3_all
x_fl = X.flatten().reshape([X.size, 1])
y_fl = Y.flatten().reshape([Y.size, 1])
z_ones = np.ones([X.size,1])
XX = np.hstack((x_fl**2, y_fl**2, x_fl, y_fl, z_ones))
In [15]:
Z = cf_fields.ravel()
Z = Z.reshape([X.size,-1])
#A_lsq = np.linalg.lstsq(XX,Z)[0]
q, r = np.linalg.qr(XX, 'reduced')
qTb = np.dot(q.T, Z)
A_lsq = np.linalg.lstsq(r, qTb)[0]
A_lsq = A_lsq.reshape([-1, coeff0.shape[0], coeff0.shape[1], coeff0.shape[2]])
In [110]:
print(A_lsq.shape, cf_fields.shape)
In [16]:
from pynrc import nrc_utils
A_reshape = A_lsq[:,:,yind,xind].reshape([5,11,1,1])
test1 = field_model(V2_interp, V3_interp, A_reshape, calc_legpoly=False)
#diff1 = cf_fields - test1
#print(diff1.min(), diff1.max())
In [17]:
cf_grid2 = test1.reshape([132,132,-1])
In [27]:
ind_nan = np.isnan(cf_grid)
cf_grid[ind_nan] = cf_grid2[ind_nan]
diff = cf_grid - cf_grid2
In [32]:
plt.plot(diff[10:100,10:100,0].flatten())
plt.plot(cf_grid[10:100,10:100,0].flatten())
plt.plot(cf_grid2[10:100,10:100,0].flatten())
Out[32]:
In [26]:
extent = np.array([xgrid.min(),xgrid.max(),ygrid.min(),ygrid.max()])
plt.imshow(cf_grid[:,:,0], extent=extent)
Out[26]:
In [ ]: