In [1]:
import time
import os
import sys
import numpy as np
import uncertainties
#import ipympl
import matplotlib
matplotlib.use('nbagg')
#from matplotlib import style
#style.use('ggplot')
import matplotlib.pyplot as plt
import astropy.units as u
from astropy import stats
from astropy.io import fits
from mmtwfs.wfs import *
from mmtwfs.zernike import *
from mmtwfs.telescope import MMT
from pathlib import Path
from astropy.modeling import models, fitting
%load_ext autoreload
%autoreload 2
import lmfit
In [ ]:
plt.close('all')
home = Path(os.environ['HOME'])
mmirs = WFSFactory(wfs="mmirs")
f9wfs = WFSFactory(wfs="newf9")
f5wfs = WFSFactory(wfs="f5")
In [ ]:
%%prun
mmirs_file = home / "MMT/mmtwfs/mmtwfs/data/test_data/mmirs_wfs_0150.fits"
mmirs_results = mmirs.measure_slopes(mmirs_file, plot=True)
mmirs_results['figures']['slopes'].show()
In [ ]:
%%prun
zresults = mmirs.fit_wavefront(mmirs_results)
zvec = zresults['zernike']
print(zresults['residual_rms'])
print(zvec.pretty_print())
zresults['resid_plot'].show()
In [ ]:
zvec.fringe_bar_chart().show()
In [ ]:
f9_file = home / "MMT/mmtwfs/mmtwfs/data/test_data/test_newf9.fits"
f9_results = f9wfs.measure_slopes(f9_file, 'blue', plot=True)
f9_results['figures']['slopes'].show()
In [ ]:
%%prun
zresults = f9wfs.fit_wavefront(f9_results)
zvec = zresults['zernike']
print(zresults['residual_rms'])
print(zvec.pretty_print())
zresults['resid_plot'].show()
In [ ]:
psf, psf_fig = f5wfs.telescope.psf(zv=zvec)
psf_fig.show()
In [ ]:
f5_file = home / "MMT/mmtwfs/mmtwfs/data/test_data/auto_wfs_0037_ave.fits"
f5_results = f5wfs.measure_slopes(f5_file, 'hecto', plot=True)
f5_results['figures']['slopes'].show()
In [ ]:
zresults = f5wfs.fit_wavefront(f5_results)
zvec = zresults['zernike']
print(zresults['residual_rms'])
print(zvec.pretty_print())
zresults['resid_plot'].show()
In [ ]:
zvec.fringe_bar_chart().show()
In [ ]:
def make_init_pars(nmodes=21, modestart=2, init_zv=None):
pars = []
for i in range(modestart, modestart+nmodes, 1):
key = "Z{:02d}".format(i)
if init_zv is not None:
val = init_zv[key].value
if val < 1.e-4:
val = 0.0
else:
val = 0.0
zpar = (key, val)
pars.append(zpar)
params = lmfit.Parameters()
params.add_many(*pars)
return params
In [ ]:
def zernike_slopes(parsdict, rho, phi, norm=False):
xslope = 0.
yslope = 0.
for k, v in parsdict.items():
l = int(k.replace("Z", ""))
dwx, dwy = zernike_slope_noll(l, rho, phi, norm=norm)
xslope += v * dwx
yslope += v * dwy
return xslope, yslope
In [ ]:
def slope_chisq(pars, coords, slopes, norm=False):
parsdict = pars.valuesdict()
rho, phi = cart2pol(coords)
xslope = slopes[0]
yslope = slopes[1]
pred_xslope, pred_yslope = zernike_slopes(parsdict, rho, phi, norm=norm)
dist = np.sqrt((xslope - pred_xslope)**2 + (yslope - pred_yslope)**2)
return dist
In [ ]:
#slopes = -f9wfs.tiltfactor * f9_results['slopes']
#pup_coords = f9wfs.modes['blue']['reference'].pup_coords(f9wfs.pup_size/2.)
In [ ]:
slopes = -f5wfs.tiltfactor * f5_results['slopes']
pup_coords = f5wfs.modes['hecto']['reference'].pup_coords(f5wfs.pup_size/2.)
In [ ]:
params = make_init_pars(nmodes=21)
out = lmfit.minimize(slope_chisq, params, args=(pup_coords, slopes))
print(lmfit.fit_report(out))
In [ ]:
zvec = ZernikeVector(coeffs=out) - f5wfs.reference_aberrations('hecto')
#zvec.units = u.m
zvec['Z23'] = 10.
print(zvec.pretty_print(last=37))
In [ ]:
zvec.fringe_bar_chart().show()
In [ ]:
rho, phi = cart2pol(pup_coords)
pred_slopes = np.array(zernike_slopes(ZernikeVector(coeffs=out), rho, phi))
diff = slopes - pred_slopes
rms = (1. / f5wfs.telescope.nmperasec) * np.sqrt((diff[0]**2 + diff[1]**2).mean())
rms
In [ ]:
for p, v in out.params.items():
print(v.stderr)
In [ ]:
(1 * u.nm)**1.3
In [ ]:
zvec.coeffs
In [ ]:
zresults = f9wfs.fit_wavefront(f9_results, plot=True)
print(zresults['residual_rms'])
print(zresults['raw_zernike'].pretty_print())
zresults['resid_plot'].show()
In [ ]:
"{0:>0.04g} +/- {1:0.04g}".format(zvec['Z04'].value, zvec['Z09'])
In [ ]:
testdata, hdr = check_wfsdata("/Users/tim/mmirs_wfs_0120.fits", header=True)
In [ ]:
mmirs.reference_aberrations('mmirs1', hdr=hdr)
In [ ]:
610/40.8
In [ ]:
np.sqrt((415/3)**2 + (1528/3)**2)
In [ ]:
np.sqrt(67**2 + 7**2)
In [ ]:
mmirs.focal_plane_position(hdr)
In [ ]:
zvec = mmirs.reference_aberrations('mmirs1', hdr=hdr)
zvec
In [ ]:
zvec.rotate(angle=-90*u.deg)
zvec
In [ ]:
from uncertainties import ufloat
In [ ]:
u.Quantity(ufloat(100., 10.), u.nm)
In [ ]:
import unyt
In [ ]:
type(ufloat(100., 10.))
In [ ]:
u.Quantity(ufloat(100., 10.), u.nm)
In [ ]:
def ZernikeGrad(zpars, x, y, atype):
if(len(zpars) > 21):
print('ZernikeGrad() is not implemented with >21 terms')
return
x2 = x * x
y2 = y * y
xy = x * y
r2 = x2 + y2
if (atype == 'dx'):
d = 0. * x # to make d an array with the same size as x
d = d + zpars['Z02'] * 1.
d = d + zpars['Z03'] * 0.
d = d + zpars['Z04'] * 4. * x
d = d + zpars['Z05'] * 2. * y
d = d + zpars['Z06'] * 2. * x
d = d + zpars['Z07'] * 6. * xy
d = d + zpars['Z08'] * (9. * x2 + 3. * y2 - 2.)
d = d + zpars['Z09'] * 6. * xy
d = d + zpars['Z10'] * (3. * x2 - 3. * y2)
d = d + zpars['Z11'] * 12. * x * (2. * (x2 + y2) - 1.)
d = d + zpars['Z12'] * x * (16. * x2 - 6.)
d = d + zpars['Z13'] * y * (24. * x2 + 8. * y2 - 6.)
d = d + zpars['Z14'] * 4. * x * (x2 - 3. * y2)
d = d + zpars['Z15'] * 4. * y * (3. * x2 - y2)
d = d + zpars['Z16'] * (x2 * (50. * x2 + 60. * y2 - 36.) + y2 * (10. * y2 - 12.) + 3.)
d = d + zpars['Z17'] * (xy * (40. * r2 - 24.))
d = d + zpars['Z18'] * (x2 * (25. * x2 - 12. - 30. * y2) + y2 * (12. - 15. * y2))
d = d + zpars['Z19'] * (4. * xy * (-6. + 15. * x2 + 5. * y2))
d = d + zpars['Z20'] * 5. * (x2 * (x2 - 6. * y2) + y2 * y2)
d = d + zpars['Z21'] * 20. * xy * (x2 - y2)
d = d + zpars['Z22'] * 24. * x * (1. + x2 * (10. * y2 - 5. + 5. * x2) + y2 * (5. * y2 - 5.))
elif (atype, 'dy'):
d = 0. * y
d = d + zpars['Z02'] * 0.
d = d + zpars['Z03'] * 1.
d = d + zpars['Z04'] * 4. * y
d = d + zpars['Z05'] * 2. * x
d = d + zpars['Z06'] * (-2.) * y
d = d + zpars['Z07'] * (3. * x2 + 9. * y2 - 2.)
d = d + zpars['Z08'] * 6. * xy
d = d + zpars['Z09'] * (3. * x2 - 3. * y2)
d = d + zpars['Z10'] * (-6.) * xy
d = d + zpars['Z11'] * 12. * y * (2. * (x2 + y2) - 1.)
d = d + zpars['Z12'] * y * (6. - 16. * y2)
d = d + zpars['Z13'] * x * (8. * x2 + 24. * y2 - 6.)
d = d + zpars['Z14'] * 4. * y * (y2 - 3. * x2)
d = d + zpars['Z15'] * 4. * x * (x2 - 3. * y2)
d = d + zpars['Z16'] * (xy * (40. * r2 - 24.))
d = d + zpars['Z17'] * (x2 * (10. * x2 + 60. * y2 - 12.) + y2 * (50. * y2 - 36.) + 3.)
d = d + zpars['Z18'] * (4. * xy * (6. - 5. * x2 - 15. * y2))
d = d + zpars['Z19'] * (y2 * (-25. * y2 + 12. + 30. * x2) + x2 * (-12. + 15. * x2))
d = d + zpars['Z20'] * 20. * xy * (y2 - x2)
d = d + zpars['Z21'] * 5. * (x2 * (x2 - 6. * y2) + y2 * y2)
d = d + zpars['Z22'] * 24. * y * (1. + y2 * (10. * x2 - 5. + 5. * y2) + x2 * (5. * x2 - 5.))
return d
In [ ]:
from math import factorial as fac
In [ ]:
rho = np.arange(10)/10
m, n = 0, 6
In [ ]:
karr = np.arange(int((n-m)/2) + 1)
karr
In [ ]:
wf = 0.0
for k in range(int((n - m)/2) + 1):
wf += rho**(n - 2.0*k) * (-1.0)**k * fac(n-k) / (fac(k) * fac((n + m)/2.0 - k) * fac((n - m)/2.0 - k))
In [ ]:
wf
In [ ]:
f = lambda k: rho**(n - 2.0*k) * (-1.0)**k * fac(n-k) / (fac(k) * fac((n + m)/2.0 - k) * fac((n - m)/2.0 - k))
ws = np.fromiter((f(ki) for ki in karr), karr.dtype, count=len(karr))
In [28]:
f9_fl = 45 * u.mm
f9_pitch = 625 * u.um
oldf9_pix = 20 * u.um
newf9_pix = 5.4 * u.um * 3
f5_pix = 20 * u.um
f5_pitch = 0.6 * u.mm
f5_fl = 40 * u.mm
mmirs_pitch = 0.6 * u.mm
mmirs_fl = 40 * u.mm
mmirs_pix = 13 * u.um * 2
def ref_width(d, fl, pix):
w = 550 * u.nm # default wavelength
theta = 1.028 * w / d
res = np.arctan(theta).value * fl
res_pix = res.to(u.um).value / pix.to(u.um).value
return res_pix
In [29]:
ref_width(f9_pitch, f9_fl, oldf9_pix)
Out[29]:
In [30]:
ref_width(mmirs_pitch, mmirs_fl, mmirs_pix)
Out[30]:
In [31]:
ref_width(f5_pitch, f5_fl, f5_pix)
Out[31]:
In [32]:
ref_width(f9_pitch, f9_fl, newf9_pix)
Out[32]:
In [ ]: