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
from matplotlib.backends.backend_pdf import PdfPages
In [2]:
import pynrc
from pynrc import nrc_utils
from pynrc.nrc_utils import S
from pynrc.obs_nircam import model_to_hdulist
pynrc.setup_logging('WARNING', verbose=False)
In [77]:
from pynrc.nrc_utils import pad_or_cut_to_size
In [3]:
fov_pix = 65
oversample = 2
psf_info={'fov_pix':fov_pix, 'oversample':oversample,
'save':False, 'force':True, 'jitter':None}
In [4]:
cf_ote = nrc_utils.psf_coeff('F444W', **psf_info)
psf_ote, psf_ote_over = nrc_utils.gen_image_coeff('F444W', coeff=cf_ote,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [5]:
cf_off = nrc_utils.psf_coeff('F444W', 'CIRCLYOT', **psf_info)
psf_off, psf_off_over = nrc_utils.gen_image_coeff('F444W', 'CIRCLYOT', coeff=cf_off,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [6]:
cf_cen = nrc_utils.psf_coeff('F444W', pupil='CIRCLYOT', mask='MASK335R', **psf_info)
psf_cen, psf_cen_over = nrc_utils.gen_image_coeff('F444W', 'CIRCLYOT', 'MASK335R', coeff=cf_cen,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [7]:
fov_asec = fov_pix*nrc_utils.pixscale_LW
pixscale = nrc_utils.pixscale_LW / oversample
In [8]:
shift_array = np.arange(0,63,3)
shift_asec = shift_array * pixscale
In [9]:
def mshift(mask, shift_array):
masks = []
for shift in shift_array:
mask_shift = nrc_utils.fshift(mask, delx=-shift)
mask_cut = nrc_utils.pad_or_cut_to_size(mask_shift, psf_ote_over.shape)
masks.append(mask_cut)
return masks
In [10]:
mask = nrc_utils.coron_trans('MASK335R', fov=2*fov_asec, pixscale=pixscale)
masks = mshift(mask, shift_array)
In [11]:
shift_array
Out[11]:
In [12]:
plt.imshow(masks[-2])
Out[12]:
In [13]:
plt.imshow((psf_ote_over*masks[3])**0.2)
Out[13]:
In [14]:
psf_list = []
psf_over_list = []
for i, offset in enumerate(shift_asec):
print(shift_array[i])
nrc = pynrc.NIRCam('F444W', pupil='CIRCLYOT', mask='MASK335R',
offset_r=offset, offset_theta=-90,
fov_pix=2*fov_pix+1, oversample=2,
force=True, save=False)
psf0, psf1 = nrc.gen_psf(return_oversample=True)
psf_list.append(psf0)
psf_over_list.append(psf1)
In [15]:
plt.imshow(psf_over_list[-1])
Out[15]:
In [16]:
for i, psf in enumerate(psf_over_list):
psf_list[i] = nrc_utils.fshift(psf_list[i], delx=-shift_array[i]/2)
psf_over_list[i] = nrc_utils.fshift(psf_over_list[i], delx=-shift_array[i]*oversample/2)
In [17]:
masks[0].shape
Out[17]:
In [18]:
fig, axes = plt.subplots(3, 7, figsize=(13,5.5))
extent = np.array([-1,1,-1,1])*pixscale*masks[0].shape[0]/2
xylim = [-1.5,1.5]
for i, ax in enumerate(axes[0]):
ax.imshow(masks[i], extent=extent)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel('MASK335R')
ax.set_title('{:.2f} arcsec'.format(shift_asec[i]))
extent = np.array([-1,1,-1,1])*pixscale*psf_ote_over.shape[0]/2
for i, ax in enumerate(axes[1]):
ax.imshow((masks[i]*psf_ote_over)**0.2, extent=extent, vmin=0, vmax=psf_ote_over.max()**0.2)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel(r'OTE PSF $\times$ Mask ')
extent = np.array([-1,1,-1,1])*pixscale*psf_over_list[0].shape[0]/2
for i, ax in enumerate(axes[2]):
ax.imshow((psf_over_list[i]**0.2), extent=extent, vmin=0, vmax=psf_off_over.max()**0.2)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel('PSF @ Detector')
fig.tight_layout()
outdir = '/Users/Jarron/Desktop/JWST-PP/'
#fig.savefig(outdir+'M335R_PSF_offset.pdf')
In [96]:
from pynrc.maths.image_manip import align_LSQ, fourier_imshift, fshift
In [ ]:
from scipy.optimize import leastsq, least_squares
def f(x, a, b):
return a*x[0] + b*x[1]
def residual(p, x, y):
return (y - f(x, *p)).ravel()
In [19]:
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)
xdata = np.array([psf_off_over, psf_cen_over]) * 1e5
avals = []
bvals = []
for i, psf in enumerate(psf_over_list):
#print(shift_array[i])
#ydata = fshift(psf, -0.00957687, -0.01525233)
ydata = nrc_utils.pad_or_cut_to_size(psf, psf_off_over.shape) * 1e5
res_best = [0,0]
chi_fin = 1e10
a0 = np.interp(shift_array[i], xv, mask[nx//2,:]**2)
b0 = 1-a0
for a0 in np.linspace(0,1,10):
b0 = 1-a0
res = least_squares(residual, [a0,b0], args=(xdata, ydata), diff_step=0.1,
loss='soft_l1', f_scale=1.0, bounds=([0,0],[1,1]))
chi_new = np.sum((ydata - f(xdata, *res.x))**2)
if chi_new < chi_fin:
chi_fin = chi_new
res_best = res.x
# ymod = f(xdata, *res_best)
# diff = ydata - ymod
# diff_rebin = nrc_utils.frebin(ydata, scale=1/oversample) - nrc_utils.frebin(ymod, scale=1/oversample)
# print(np.abs(diff_rebin**2).sum(), np.abs(diff_rebin).max())
ymod = f(xdata, *res_best) / 1e5
ydata /= 1e5
diff = ydata - ymod
ydata_rebin = nrc_utils.frebin(ydata, scale=1/oversample)
ymod_rebin = nrc_utils.frebin(ymod, scale=1/oversample)
pnoise = np.sqrt(ymod_rebin)
ind = ymod_rebin>0
diff_rebin = ymod_rebin - ydata_rebin
chi1 = np.sum((diff_rebin[ind])**2 / ymod_rebin[ind])
chi2 = np.sum(diff_rebin[ind]**2 / pnoise[ind]**2)
print(shift_array[i], chi1, chi2/ymod_rebin.sum(), np.abs(diff_rebin).max(), res_best)
avals.append(res_best[0])
bvals.append(res_best[1])
In [29]:
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)
xdata = np.array([psf_off_over, psf_cen_over])
for i, psf in enumerate(psf_over_list):
#print(shift_array[i])
#ydata = fshift(psf, -0.00957687, -0.01525233)
ydata = nrc_utils.pad_or_cut_to_size(psf, psf_off_over.shape)
res_best = [0,0]
chi_fin = 1e10
a0 = np.interp(shift_array[i], xv, mask[nx//2,:]**2)
b0 = 1-np.interp(shift_array[i], xv, mask[nx//2,:]**4)
res_best = (a0, b0)
ymod = f(xdata, *res_best)
diff = ydata - ymod
ydata_rebin = nrc_utils.frebin(ydata, scale=1/oversample)
ymod_rebin = nrc_utils.frebin(ymod, scale=1/oversample)
pnoise = np.sqrt(ymod_rebin)
ind = ymod_rebin>0
diff_rebin = ymod_rebin - ydata_rebin
chi1 = np.sum((diff_rebin[ind])**2 / ymod_rebin[ind])
chi2 = np.sum(diff_rebin[ind]**2 / pnoise[ind]**2)
print(shift_asec[i], chi1, chi2/ymod_rebin.sum(), np.abs(diff_rebin).max(), res_best)
In [30]:
fig,ax = plt.subplots(1,1, figsize=(8,5))
ax.plot(shift_asec, avals, label=r'Off-axis PSF scale (a vals)')
ax.plot(shift_asec, bvals, label=r'Centered PSF scale (b vals)')
#ax.plot(rad_bins[1:], rad_mn**2, label='Test')
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)*pixscale
ax.plot(xv, mask[nx//2,:]**2, label=r'Squared Throughput ($T^2$)')
ax.plot(xv, 1-mask[nx//2,:]**2, label=r'$1-T^2$')
ax.set_xlim([-0.1,2])
ax.set_xlabel('Radial Offset (arcsec)')
ax.set_ylabel('Scale Factor')
ax.set_title('PSF Coefficient Relationship')
ax.legend()
fig.tight_layout()
#fig.savefig(outdir+'psf_coeff_relation.pdf')
In [195]:
fig, ax = plt.subplots(1,3, figsize=(12,5))
i = 1
psf = psf_over_list[i]
ydata = nrc_utils.pad_or_cut_to_size(psf, psf_off_over.shape)
ymod = f(xdata, avals[i], bvals[i])
diff = ydata - ymod
print(np.abs(diff).max())
extent = np.array([-1,1,-1,1])*pixscale*ydata.shape[0]/2
ax[0].imshow(ydata**0.2, extent=extent)
ax[0].set_title(r'WebbPSF r_offset = {:.1f}$^{{\prime\prime}}$'.format(shift_asec[i]))
ax[1].imshow(ymod**0.2, extent=extent)
ax[1].set_title('Linear Combination')
ax[2].imshow(diff, extent=extent)
ax[2].set_title('Residuals')
mm_str = r'min, max = ({:.1e}, {:.1e})'.format(diff.min(),diff.max())
ax[2].text(-1.9, 1.6, mm_str, fontsize=12)
fig.tight_layout()
fig.savefig(outdir + 'psf_model.pdf')
In [ ]:
In [3]:
fov_pix = 129
oversample = 2
psf_info={'fov_pix':fov_pix, 'oversample':oversample,
'save':False, 'force':True, 'jitter':None}
In [4]:
filt, pupil, pmask = ('F430M', 'WEDGELYOT', 'MASKLWB')
#filt, pupil, pmask = ('F430M', 'CIRCLYOT', 'MASK335R')
In [5]:
cf_ote = nrc_utils.psf_coeff(filt, **psf_info)
psf_ote, psf_ote_over = nrc_utils.gen_image_coeff(filt, coeff=cf_ote,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [6]:
cf_off = nrc_utils.psf_coeff(filt, pupil=pupil, **psf_info)
psf_off, psf_off_over = nrc_utils.gen_image_coeff(filt, pupil=pupil, coeff=cf_off,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [7]:
cf_cen = nrc_utils.psf_coeff(filt, pupil=pupil, mask=pmask, **psf_info)
psf_cen, psf_cen_over = nrc_utils.gen_image_coeff(filt, pupil=pupil, mask=pmask, coeff=cf_cen,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [8]:
fov_asec = fov_pix*nrc_utils.pixscale_LW
pixscale = nrc_utils.pixscale_LW / oversample
In [9]:
shift_array = np.arange(0,84,3)
shift_asec = shift_array * pixscale
In [10]:
def mshift(mask, shift_array):
masks = []
for shift in shift_array:
mask_shift = nrc_utils.fshift(mask, dely=shift)
mask_cut = nrc_utils.pad_or_cut_to_size(mask_shift, psf_ote_over.shape)
masks.append(mask_cut)
return masks
In [11]:
mask = nrc_utils.coron_trans(pmask, fov=2*fov_asec, pixscale=pixscale)
masks = mshift(mask, shift_array)
In [12]:
plt.imshow(masks[-5])
Out[12]:
In [13]:
psf_list = []
psf_over_list = []
for i, offset in enumerate(shift_asec):
print(shift_array[i])
nrc = pynrc.NIRCam(filt, pupil=pupil, mask=pmask,
offset_r=offset, offset_theta=0,
fov_pix=2*fov_pix+1, oversample=oversample,
force=True, save=False)
psf0, psf1 = nrc.gen_psf(return_oversample=True)
psf_list.append(psf0)
psf_over_list.append(psf1)
In [14]:
plt.imshow(psf_over_list[-5])
Out[14]:
In [15]:
for i, psf in enumerate(psf_over_list):
psf_list[i] = nrc_utils.fshift(psf_list[i], dely=-shift_array[i]/2)
psf_over_list[i] = nrc_utils.fshift(psf_over_list[i], dely=-shift_array[i]*oversample/2)
In [16]:
masks[0].shape
Out[16]:
In [17]:
fig, axes = plt.subplots(3, 7, figsize=(13,5.5))
extent = np.array([-1,1,-1,1])*pixscale*masks[0].shape[0]/2
xylim = [-2,2]
for i, ax in enumerate(axes[0]):
ax.imshow(masks[i], extent=extent)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel(pmask)
ax.set_title('{:.2f} arcsec'.format(shift_asec[i]))
extent = np.array([-1,1,-1,1])*pixscale*psf_ote_over.shape[0]/2
for i, ax in enumerate(axes[1]):
ax.imshow((masks[i]*psf_ote_over)**0.2, extent=extent, vmin=0, vmax=psf_ote_over.max()**0.2)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel(r'OTE PSF $\times$ Mask ')
extent = np.array([-1,1,-1,1])*pixscale*psf_over_list[0].shape[0]/2
for i, ax in enumerate(axes[2]):
ax.imshow((psf_over_list[i]**0.2), extent=extent, vmin=0, vmax=psf_off_over.max()**0.2)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel('PSF @ Detector')
fig.tight_layout()
#outdir = '/Users/Jarron/Desktop/JWST-PP/'
#fig.savefig(outdir+'M335R_PSF_offset.pdf')
In [18]:
from pynrc.maths.image_manip import align_LSQ, fourier_imshift, fshift
In [19]:
from scipy.optimize import leastsq, least_squares
def f(x, a, b):
return a*x[0] + b*x[1]
def residual(p, x, y):
return (y - f(x, *p)).ravel()
In [26]:
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)
xdata = np.array([nrc_utils.pad_or_cut_to_size(psf_off_over, (64,64)),
nrc_utils.pad_or_cut_to_size(psf_cen_over, (64,64))])
avals = []
bvals = []
for i, psf in enumerate(psf_over_list):
#print(shift_array[i])
#ydata = fshift(psf, -0.00957687, -0.01525233)
ydata = nrc_utils.pad_or_cut_to_size(psf, (64,64))
res_best = [0,0]
chi_fin = 1e10
a0 = np.interp(shift_array[i], xv, mask[:,nx//2]**2)
b0 = 1-a0
for a0 in np.linspace(0,1,10):
b0 = 1-a0
res = least_squares(residual, [a0,b0], args=(xdata, ydata), diff_step=0.1,
bounds=([0,0],[1,1]))
chi_new = np.sum((ydata - f(xdata, *res.x))**2)
if chi_new < chi_fin:
chi_fin = chi_new
res_best = res.x
ymod = f(xdata, *res_best)
diff = ydata - ymod
ydata_rebin = nrc_utils.frebin(ydata, scale=1/oversample)
ymod_rebin = nrc_utils.frebin(ymod, scale=1/oversample)
pnoise = np.sqrt(ymod_rebin)
ind = ymod_rebin>0
diff_rebin = ymod_rebin - ydata_rebin
chi1 = np.sum((diff_rebin[ind])**2 / ymod_rebin[ind])
chi2 = np.sum(diff_rebin[ind]**2 / pnoise[ind]**2)
print(shift_array[i], chi1, chi2/ydata_rebin.sum(), np.abs(diff_rebin).max(), res_best)
avals.append(res_best[0])
bvals.append(res_best[1])
In [27]:
fig,ax = plt.subplots(1,1, figsize=(8,5))
ax.plot(shift_asec, avals, label=r'Off-axis PSF scale (a vals)')
ax.plot(shift_asec, bvals, label=r'Centered PSF scale (b vals)')
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)*pixscale
ax.plot(xv, mask[:,nx//2]**2, label=r'Squared Throughput ($T^2$)')
ax.plot(xv, 1-mask[:,nx//2]**2, label=r'$1-T^2$')
ax.set_xlim([-0.1,2.5])
ax.set_xlabel('Radial Offset (arcsec)')
ax.set_ylabel('Scale Factor')
ax.set_title('PSF Coefficient Relationship')
ax.legend()
fig.tight_layout()
#fig.savefig(outdir+'psf_coeff_relation.pdf')
In [44]:
fig, ax = plt.subplots(1,3, figsize=(12,5))
i = 5
xdata = np.array([psf_off_over, psf_cen_over])
psf = psf_over_list[i]
ydata = nrc_utils.pad_or_cut_to_size(psf, psf_off_over.shape)
ymod = f(xdata, avals[i], bvals[i])
diff = ydata - ymod
extent = np.array([-1,1,-1,1])*pixscale*ydata.shape[0]/2
ax[0].imshow(ydata**0.2, extent=extent, vmin=0, vmax=ydata.max()**0.2)
ax[0].set_title(r'WebbPSF r_offset = {:.1f}$^{{\prime\prime}}$'.format(shift_asec[i]))
ax[1].imshow(ymod**0.2, extent=extent, vmin=0, vmax=ydata.max()**0.2)
ax[1].set_title('Linear Combination')
ax[2].imshow(diff, extent=extent)
ax[2].set_title('Residuals')
print(np.abs(diff).max())
mm_str = r'min, max = ({:.1e}, {:.1e})'.format(diff.min()/ydata.max(), diff.max()/ydata.max())
ax[2].text(-2, 1.6, mm_str, fontsize=12)
fig.tight_layout()
#fig.savefig(outdir + 'psf_model.pdf')
In [45]:
fov_pix = 321
oversample = 2
psf_info={'fov_pix':fov_pix, 'oversample':oversample,
'save':False, 'force':True, 'jitter':None}
In [46]:
filt, pupil, pmask = ('F430M', 'WEDGELYOT', 'MASKLWB')
#filt, pupil, pmask = ('F430M', 'CIRCLYOT', 'MASK335R')
In [47]:
cf_ote = nrc_utils.psf_coeff(filt, **psf_info)
psf_ote, psf_ote_over = nrc_utils.gen_image_coeff(filt, coeff=cf_ote,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [48]:
cf_off = nrc_utils.psf_coeff(filt, pupil=pupil, **psf_info)
psf_off, psf_off_over = nrc_utils.gen_image_coeff(filt, pupil=pupil, coeff=cf_off,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [49]:
cf_cen = nrc_utils.psf_coeff(filt, pupil=pupil, mask=pmask, **psf_info)
psf_cen, psf_cen_over = nrc_utils.gen_image_coeff(filt, pupil=pupil, mask=pmask, coeff=cf_cen,
fov_pix=fov_pix, oversample=oversample,
return_oversample=True)
In [50]:
fov_asec = fov_pix*nrc_utils.pixscale_LW
pixscale = nrc_utils.pixscale_LW / oversample
In [51]:
shift_array = np.arange(-288,300,12)
shift_asec = shift_array * pixscale
print(shift_asec)
In [52]:
def mshift(mask, shift_array):
masks = []
for shift in shift_array:
mask_shift = nrc_utils.fshift(mask, delx=-shift)
mask_cut = nrc_utils.pad_or_cut_to_size(mask_shift, psf_ote_over.shape)
masks.append(mask_cut)
return masks
In [53]:
mask = nrc_utils.coron_trans(pmask, fov=2*fov_asec, pixscale=pixscale)
masks = mshift(mask, shift_array)
In [54]:
plt.imshow(masks[0])
Out[54]:
In [55]:
psf_list = []
psf_over_list = []
for i, offset in enumerate(shift_asec):
print(shift_array[i])
nrc = pynrc.NIRCam(filt, pupil=pupil, mask=pmask,
offset_r=offset, offset_theta=-90,
fov_pix=401, oversample=oversample,
force=True, save=False)
psf0, psf1 = nrc.gen_psf(return_oversample=True)
psf_list.append(psf0)
psf_over_list.append(psf1)
In [67]:
for psf in psf_over_list:
print(psf.sum())
In [56]:
for i, psf in enumerate(psf_over_list):
psf_list[i] = nrc_utils.fshift(psf_list[i], delx=-shift_array[i]/2)
psf_over_list[i] = nrc_utils.fshift(psf_over_list[i], delx=-shift_array[i]*oversample/2)
In [57]:
fig, axes = plt.subplots(3, 7, figsize=(13,5.5))
extent = np.array([-1,1,-1,1])*pixscale*masks[0].shape[0]/2
xylim = [-2,2]
for i, ax in enumerate(axes[0]):
ii = 5*i
ax.imshow(masks[ii], extent=extent)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel(pmask)
ax.set_title('{:.2f} arcsec'.format(shift_asec[ii]))
extent = np.array([-1,1,-1,1])*pixscale*psf_ote_over.shape[0]/2
for i, ax in enumerate(axes[1]):
ii = 5*i
ax.imshow((masks[ii]*psf_ote_over)**0.2, extent=extent, vmin=0, vmax=psf_ote_over.max()**0.2)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel(r'OTE PSF $\times$ Mask ')
extent = np.array([-1,1,-1,1])*pixscale*psf_over_list[0].shape[0]/2
for i, ax in enumerate(axes[2]):
ii = 5*i
ax.imshow((psf_over_list[ii]**0.2), extent=extent, vmin=0, vmax=psf_cen_over.max()**0.2)
ax.set_xlim(xylim)
ax.set_ylim(xylim)
if i==0: ax.set_ylabel('PSF @ Detector')
fig.tight_layout()
#outdir = '/Users/Jarron/Desktop/JWST-PP/'
#fig.savefig(outdir+'M335R_PSF_offset.pdf')
In [58]:
from pynrc.maths.image_manip import align_LSQ, fourier_imshift, fshift
In [66]:
from astropy.modeling import models, fitting
ny,nx = masks[0].shape
xv = (np.arange(nx)-nx/2)*pixscale
yv = (np.arange(ny)-ny/2)*pixscale
sig_list = []
for m in masks:
mv = m[:,nx//2]
ind_temp = mv>0.98
ind = np.abs(yv) < np.min(np.abs(yv[ind_temp]))
g_init = models.Gaussian1D(amplitude=1., mean=0, stddev=1.)
fit_g = fitting.LevMarLSQFitter()
g = fit_g(g_init, yv[ind], 1-mv[ind])
sig_list.append(np.abs(g.stddev.value))
sig_vals = np.array(sig_list)
In [67]:
sig_vals
Out[67]:
In [71]:
from scipy.optimize import leastsq, least_squares
def f(x, a, b, c):
return a*x[0] + b*x[1] + c*x[2]
def residual(p, x, y):
return ((y - f(x, *p))**1).ravel()
In [72]:
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)
xdata = np.array([nrc_utils.pad_or_cut_to_size(psf_over_list[5], (64,64)),
nrc_utils.pad_or_cut_to_size(psf_cen_over, (64,64)),
nrc_utils.pad_or_cut_to_size(psf_over_list[-6], (64,64))]) * 1e5
avals = []
bvals = []
cvals = []
res0_1 = []
res0_2 = []
for i, psf in enumerate(psf_over_list):
#print(shift_array[i])
#ydata = fshift(psf, -0.00957687, -0.01525233)
#ydata = nrc_utils.pad_or_cut_to_size(psf, psf_off_over.shape)
ydata = nrc_utils.pad_or_cut_to_size(psf, (64,64))*1e5
res_best = [0,0]
chi_fin = 1e10
for b0 in np.linspace(0,2,10):
a0= (2-b0)/2
c0= 2 - a0 - b0
res = least_squares(residual, [a0, b0, c0], args=(xdata, ydata), diff_step=0.01,
bounds=([0,0,0],[2,2,2]), method='trf', verbose=0)
chi_new = np.sum((ydata - f(xdata, *res.x))**2)
if chi_new < chi_fin:
chi_fin = chi_new
res_best = res.x
ymod = f(xdata, *res_best) / 1e5
ydata /= 1e5
diff = ydata - ymod
ydata_rebin = nrc_utils.frebin(ydata, scale=1/oversample)
ymod_rebin = nrc_utils.frebin(ymod, scale=1/oversample)
pnoise = np.sqrt(ymod_rebin)
ind = ymod_rebin>0
diff_rebin = ymod_rebin - ydata_rebin
chi1 = np.sum((diff_rebin[ind])**2 / ymod_rebin[ind])
chi2 = np.sum(diff_rebin[ind]**2 / pnoise[ind]**2)
print(shift_array[i], chi1, chi1/ydata_rebin.sum(), np.abs(diff_rebin).max(), res_best)
#print(shift_array[i], np.abs(diff_rebin**2).sum(), np.abs(diff_rebin).max(), res_best)
avals.append(res_best[0])
bvals.append(res_best[1])
cvals.append(res_best[2])
res0_1.append(np.abs(diff_rebin**2).sum())
res0_2.append(np.abs(diff_rebin).max())
In [73]:
fig,ax = plt.subplots(1,1, figsize=(8,5))
ax.plot(shift_asec, avals, label=r'Off-axis PSF scale (a vals)')
ax.plot(shift_asec, bvals, label=r'Centered PSF scale (b vals)')
ax.plot(shift_asec, cvals, label=r'Centered PSF scale (c vals)')
ax.plot(shift_asec, sig_vals, label=r'Sigma')
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)*pixscale
#ax.plot(xv, mask[:,nx//2]**2, label=r'Squared Throughput ($T^2$)')
#ax.plot(xv, 1-mask[:,nx//2]**2, label=r'$1-T^2$')
ax.set_xlim([-10,10])
ax.set_xlabel('Radial Offset (arcsec)')
ax.set_ylabel('Scale Factor')
ax.set_title('PSF Coefficient Relationship')
ax.legend()
fig.tight_layout()
#fig.savefig(outdir+'psf_coeff_relation.pdf')
In [75]:
psf_over_list[5].shape
Out[75]:
In [76]:
psf_cen_over.shape
Out[76]:
In [81]:
plt.imshow(psf_over_list[5]**0.2)
Out[81]:
In [84]:
fig, ax = plt.subplots(1,3, figsize=(12,5))
i = 40
sh = (256,256)
xdata = np.array([pad_or_cut_to_size(psf_over_list[5], sh),
pad_or_cut_to_size(psf_cen_over, sh),
pad_or_cut_to_size(psf_over_list[-6], sh)])
psf = psf_over_list[i]
ydata = pad_or_cut_to_size(psf, sh)
ymod = f(xdata, avals[i], bvals[i], cvals[i])
diff = ydata - ymod
extent = np.array([-1,1,-1,1])*pixscale*ydata.shape[0]/2
ax[0].imshow(ydata**0.2, extent=extent, vmin=0, vmax=ydata.max()**0.2)
ax[0].set_title(r'WebbPSF r_offset = {:.1f}$^{{\prime\prime}}$'.format(shift_asec[i]))
ax[1].imshow(ymod**0.2, extent=extent, vmin=0, vmax=ydata.max()**0.2)
ax[1].set_title('Linear Combination')
ax[2].imshow(diff, extent=extent)
ax[2].set_title('Residuals')
print(np.abs(diff).max())
mm_str = r'min, max = ({:.1e}, {:.1e})'.format(diff.min()/ydata.max(), diff.max()/ydata.max())
ax[2].text(-2, 1.6, mm_str, fontsize=12)
fig.tight_layout()
#fig.savefig(outdir + 'psf_model.pdf')
In [64]:
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)
xdata = np.array([nrc_utils.pad_or_cut_to_size(psf_over_list[5], (64,64)),
nrc_utils.pad_or_cut_to_size(psf_cen_over, (64,64)),
nrc_utils.pad_or_cut_to_size(psf_over_list[-6], (64,64))]) * 1e5
avals = []
bvals = []
cvals = []
res0_1 = []
res0_2 = []
for i, psf in enumerate(psf_over_list):
#print(shift_array[i])
#ydata = fshift(psf, -0.00957687, -0.01525233)
#ydata = nrc_utils.pad_or_cut_to_size(psf, psf_off_over.shape)
ydata = nrc_utils.pad_or_cut_to_size(psf, (64,64))*1e5
res_best = [0,0]
chi_fin = 1e10
for a0 in np.linspace(0,2,10):
b0= (2 - a0) / 2
c0= 2 - a0 - b0
res = least_squares(residual, [a0, 0, c0], args=(xdata, ydata), diff_step=0.01,
bounds=([0,0,0],[2,1e-10,2]), method='trf', verbose=0)
chi_new = np.sum((ydata - f(xdata, *res.x))**2)
if chi_new < chi_fin:
chi_fin = chi_new
res_best = res.x
ymod = f(xdata, *res_best) / 1e5
ydata /= 1e5
diff = ydata - ymod
ydata_rebin = nrc_utils.frebin(ydata, scale=1/oversample)
ymod_rebin = nrc_utils.frebin(ymod, scale=1/oversample)
pnoise = np.sqrt(ymod_rebin)
ind = ymod_rebin>0
diff_rebin = ymod_rebin - ydata_rebin
chi1 = np.sum((diff_rebin[ind])**2 / ymod_rebin[ind])
chi2 = np.sum(diff_rebin[ind]**2 / pnoise[ind]**2)
print(shift_array[i], chi1, chi1/ydata_rebin.sum(), np.abs(diff_rebin).max(), res_best)
avals.append(res_best[0])
bvals.append(res_best[1])
cvals.append(res_best[2])
res0_1.append(np.abs(diff_rebin**2).sum())
res0_2.append(np.abs(diff_rebin).max())
In [69]:
from scipy.optimize import leastsq, least_squares
def f(x, a, b):#, c):
return a*x[0] + b*x[1] #+ c*x[2]
def residual(p, x, y):
return ((y - f(x, *p))**1).ravel()
In [70]:
nx = mask.shape[0]
xv = (np.arange(nx)-nx/2)
xdata = np.array([nrc_utils.pad_or_cut_to_size(psf_over_list[5], (64,64)),
nrc_utils.pad_or_cut_to_size(psf_over_list[-6], (64,64))])
sig_data = sig_vals**0.5
sig_data = sig_data - sig_data.min()
sig_data = sig_data / sig_data.max()
res2_1 = []
res2_2 = []
for i, psf in enumerate(psf_over_list):
ydata = nrc_utils.pad_or_cut_to_size(psf, (64,64))
a0 = sig_data[i]
b0 = 1-a0
res_best = (a0, b0)
ymod = f(xdata, *res_best)
diff = ydata - ymod
ydata_rebin = 1e10 * nrc_utils.frebin(ydata, scale=1/oversample)
ymod_rebin = 1e10 * nrc_utils.frebin(ymod, scale=1/oversample)
pnoise = np.sqrt(ymod_rebin)
ind = ymod_rebin>0
diff_rebin = ymod_rebin - ydata_rebin
chi1 = np.sum((diff_rebin[ind])**2 / ymod_rebin[ind])
chi2 = np.sum(diff_rebin[ind]**2 / pnoise[ind]**2)
print(shift_array[i], chi1, chi1/ydata_rebin.sum(), np.abs(diff_rebin).max(), res_best)
res2_1.append(np.abs(diff_rebin**2).sum())
res2_2.append(np.abs(diff_rebin).max())
In [227]:
psf_cen.sum()
Out[227]:
In [224]:
plt.semilogy(shift_asec, np.array(res0_2) / np.abs(psf_cen).max())
plt.semilogy(shift_asec, np.array(res1_2) / np.abs(psf_cen).max())
plt.semilogy(shift_asec, np.array(res2_2) / np.abs(psf_cen).max())
#plt.ylim([1e-9,1e-5])
#plt.semilogy(shift_asec, res1_2)
#plt.semilogy(shift_asec, res2_2)
Out[224]:
In [216]:
plt.semilogy(shift_asec, res0_2)
plt.semilogy(shift_asec, res1_2)
plt.semilogy(shift_asec, res2_2)
plt.ylim([1e-9,1e-5])
#plt.semilogy(shift_asec, res1_2)
#plt.semilogy(shift_asec, res2_2)
Out[216]:
In [184]:
plt.plot(yv, m[:,nx//2])
#plt.plot(xv, 1-g(xv))
Out[184]:
In [ ]:
In [ ]:
def func_opt(params, xdata, ydata):
model = xdata[0] * params[0] + ydata[0] * params[1]
return ydata - model
res = least_squares(shift_subtract, init_pars, diff_step=0.1,
loss='soft_l1', f_scale=1.0, args=(reference,target),
kwargs={'mask':mask,'pad':pad,'shift_function':shift_function})
In [86]:
res.x
Out[86]:
In [ ]:
def func_opt(params, xdata, ydata):
model = xdata[0] * params[0] + ydata[0] * params[1]
return ydata - model
#return (ydata - numpy.dot(xdata, params))
leastsq(func, x0, args=(xdata, ydata))