#!/usr/bin/env Python3 import numpy as np import matplotlib.pyplot as plt from matplotlib.ticker import FormatStrFormatter from matplotlib.colors import LogNorm import pickle import sys sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Documents/Published/FourierOptics GitHub release/') from fourierOptics import * ''' @ Author: Joost van den Born @ Description: Script for analysing WFEs produced using either Zernike polynomials or a PSD description. It also compares the shift in centroid position as the complexity of the WFE increases through spatial frequency filtering. ''' def saveVariable(var, filename): # Saving the objects: with open(filename, 'wb') as f: # Python 3: open(..., 'wb') pickle.dump(var, f) def loadVariable(filename): # Getting back the objects: with open(filename, 'rb') as f: # Python 3: open(..., 'rb') obj = pickle.load(f) return obj def phaseTilt(x,y,tx,ty): phaseTilt = x * np.tan(np.deg2rad(tx)) + y * np.tan(np.deg2rad(ty)) return phaseTilt if __name__ == '__main__': plt.rc('font', family='serif') plt.rc('axes', labelsize='large') # Pick a wavelength l = 1.35e-3 # Initiate Fourier Optics fou = FourierOptics(l=l) # Pick sampling parameters dpix = 160 n, m = 1001, 1001 oC = [-.75, .75, -.75, .75] print ('Your resolution is {}x{} with a pixel size of {}'.format(n,m,dpix)) # Initiate the profiles class and pick a profile (tophat, obscured tophat or spider) prof = Profiles() # Aideal = prof.top_hat(38854/2, n, m, dpix=dpix, oversample=1) A = prof.ELT_pupil(n,m,dpix) # Create an input grid x, y = np.meshgrid(np.linspace(-n/2*dpix, n/2*dpix, n), np.linspace(-m/2*dpix, m/2*dpix, m)) # Add PSD wavefront error to the lens wfe = WavefrontError(n,m) wfe.generate_height_map(5e-5, 2, 1) # 50 nm RMS for standard h_data = wfe.data # Setting up empty lists for fitting process popt_airy_list, popt_cust_list = [], [] pcov_airy_list, custom_var_list = [], [] strehl_list = [] Exyzlist, PSDlist, surfacelist, EEFlist = [], [], [], [] # Create a lens with a focal length ~100 mm, with a clear aperture radius of 4 mm and thickness of 1 mm. M1 = ParabolicMirror(dx=0, dy=0, name='ELT-paraxial', CA=38854, f=6.84022e5) # M1.add_wavefront_error(phaseTilt(x,y,tx=7.5, ty=7.5)) t_l = M1.phase_delay(x,y,l) U_in = A * t_l # Propagating the original, unperturbed, system and calculating encircled energy Exyz_og,xp,yp = fou.fresnel_propagation(U_in, x, y, dz=M1.f, show=False, outputCoordinates=oC) eef_og, r_og = fou.encircled_energy_function(np.abs(Exyz_og)**2, xp, yp) plt.imshow(np.abs(Exyz_og), cmap='hot');plt.colorbar();plt.show() zernikeAnalysis = True PSDAnalysis = True if zernikeAnalysis: # For some reason we need to reset the WFE object. wfe = WavefrontError(n,m) wfe.data = h_data psd1d_og = wfe.calculate_1D_PSD(data=h_data) # For each Zernike mode, perform analysis. # Don't touch the RMS values, it will converge to some value # If it goes as expected you'll see some form of convergence, # possibly a spiral, of the centroid location as you add more # and more modes. # Calculate the maximum mode that we want to fit the zernike modes to. # nmax, mmax = 80, 80 nmax, mmax = 20,20 zernikeLimit = int((nmax*(nmax+2)+mmax)/2) # Fit the zernike coefficients for all modes up to zernikeLimit print('Fitting Zernike up to j={:.0f}'.format(zernikeLimit)) zernikeCoefficients = wfe.fit_zernike_coefficients(h_data, modes=zernikeLimit) # Result lists: x0zern, y0zern = [], [] x0zern_var, y0zern_var = [], [] zernCouplingList = [] eef_zern_list = [] zernImageResponse = [] psd1d_list_zern = [] for nmax_i in range(nmax): print (' Zernike analysis: {0}/{1}'.format(nmax_i, nmax), end="\r") zernikeModeMax = int((nmax_i*(nmax_i+2)+nmax_i)/2+1) # The zernike coefficients for this particular mode zC = zernikeCoefficients[:zernikeModeMax] xz, yz = np.mgrid[-1:1:n*1j, -1:1:m*1j] wfeFilteredZern = wfe.zernike_cartesian(zC, xz, yz) # Adding the wavefront error M1.add_wavefront_error(wfeFilteredZern) wfe_t_l = M1.phase_delay(x,y,l) U_in = A * wfe_t_l # Performing the propagation using Fresnel approximation to the focal plane. Exyz_zern,xz,yz = fou.fresnel_propagation(U_in, x, y, dz=M1.f, show=False, outputCoordinates=oC) eef_zern, _ = fou.encircled_energy_function(np.abs(Exyz_zern)**2, xz, yz, originalPSF=np.abs(Exyz_og)**2) # Fitting the centroid location using the coupling integral # Set up fitting object psfFit_zern = PSFFitter(xz.real,yz.real,np.abs(Exyz_zern)**2) custFitted, popt_cust, cust_var, cust_coupling = psfFit_zern.fit_custom_field(Exyz_og, Exyz_zern, return_coupling=True) # Returning the found PSD curve psd1d_list_zern.append(wfe.calculate_1D_PSD(data=wfeFilteredZern)) x0zern.append(popt_cust[0]) y0zern.append(popt_cust[1]) x0zern_var.append(cust_var[0]) y0zern_var.append(cust_var[1]) zernCouplingList.append(cust_coupling) eef_zern_list.append(eef_zern.tolist()) zernImageResponse.append(Exyz_zern) # saveVariable(x0zern, './data/x0zern.pck') # saveVariable(y0zern, './data/y0zern.pck') # saveVariable(x0zern_var, './data/x0zern_var.pck') # saveVariable(y0zern_var, './data/y0zern_var.pck') # saveVariable(zernCouplingList, './data/zernCouplingList.pck') # saveVariable(eef_zern_list, './data/eef_zern_list.pck') # saveVariable(zernImageResponse, './data/zernImageResponse.pck') # saveVariable(psd1d_list_zern, './data/psd1d_list_zern.pck') # saveVariable(zernikeCoefficients, './data/zernikeCoefficients.pcks') if PSDAnalysis: # For each spatial frequency, perform the analysis # Don't touch the RMS value. It will converge to some value # For some reason we need to reset the WFE object. wfe = None wfe = WavefrontError(n,m) wfe.data = h_data psd1d_og = wfe.calculate_1D_PSD(data=h_data) # Result lists: x0psd, y0psd = [], [] x0psd_var, y0psd_var =[], [] psdCouplingList = [] eef_psd_list = [] psdImageResponse = [] frequencyRange = np.asarray([1,2,3,4,5,6,7,8,9] + np.linspace(10,int(np.sqrt(2)*n/2+1),70).tolist()) psd1d_list_psd = [] for lowpassFrequency in frequencyRange: i1 = frequencyRange.tolist().index(lowpassFrequency) i2 = len(frequencyRange) print (' Power Spectral Density analysis: {0}/{1}'.format(i1, i2), end="\r") # Filter the spatial frequencies wfeFilteredPSD = wfe.set_lowpass_filter(lowpassFrequency) # Adding the wavefront error M1.add_wavefront_error(wfeFilteredPSD) wfe_t_l = M1.phase_delay(x,y,l) U_in = A * wfe_t_l # Performing the propagation using Fresnel approximation to the focal plane. Exyz_psd,xp,yp = fou.fresnel_propagation(U_in, x, y, dz=M1.f, show=False, outputCoordinates=oC) eef_psd, _ = fou.encircled_energy_function(np.abs(Exyz_psd)**2, xp, yp, originalPSF=np.abs(Exyz_og)**2) # Fitting the centroid location using the coupling integral # Set up fitting object psfFit_psd = PSFFitter(xp.real,yp.real,np.abs(Exyz_psd)**2) custFitted, popt_cust, cust_var, cust_coupling = psfFit_psd.fit_custom_field(Exyz_og, Exyz_psd, return_coupling=True) # Returning the found PSD curve psd1d_list_psd.append(wfe.calculate_1D_PSD(data=wfeFilteredPSD)) x0psd.append(popt_cust[0]) y0psd.append(popt_cust[1]) x0psd_var.append(cust_var[0]) y0psd_var.append(cust_var[1]) psdCouplingList.append(cust_coupling) eef_psd_list.append(eef_psd) psdImageResponse.append(Exyz_psd) # saveVariable(x0psd, './data/x0psd.pck') # saveVariable(y0psd, './data/y0psd.pck') # saveVariable(x0psd_var, './data/x0psd_var.pck') # saveVariable(y0psd_var, './data/y0psd_var.pck') # saveVariable(psdCouplingList, './data/psdCouplingList.pck') # saveVariable(eef_psd_list, './data/eef_psd_list.pck') # saveVariable(psdImageResponse, './data/psdImageResponse.pck') # saveVariable(psd1d_list_psd, './data/psd1d_list_psd.pck') # Plot some of the results: if PSDAnalysis: # Power Spectral Density distribution for PSD analysis fig = plt.figure() ax = fig.add_subplot(111) ax.loglog(psd1d_og, 'r-') # The original PSD curve # The PSD curve for each lowpass filtered PSD WFE for i in range(len(frequencyRange)): ax.loglog(psd1d_list_psd[i], 'b-', label='data') ax.grid(True, alpha=0.3) ax.set_xlabel('Spatial frequency') ax.set_ylabel('PSD') plt.show(block=False) if zernikeAnalysis: # Power Spectral Density distribution for Zernike Analysis fig = plt.figure() ax = fig.add_subplot(111) # The original PSD curve ax.loglog(psd1d_og, 'r-') # The PSD curve for each of the zernike fits for i in range(nmax): ax.loglog(psd1d_list_zern[i], 'b-', label='data') ax.grid(True, alpha=0.3) ax.set_xlabel('Spatial frequency') ax.set_ylabel('PSD') plt.show() # Plot of the centroid shift fig = plt.figure() ax = fig.add_subplot(111) if zernikeAnalysis: ax.plot(x0zern, y0zern, 'b-', label='data') ax.plot(x0zern, y0zern, 'bx', label='data') ax.plot(x0zern[0], y0zern[0], 'rx') if PSDAnalysis: ax.plot(x0psd, y0psd, 'r-', label='data') ax.plot(x0psd, y0psd, 'rx', label='data') ax.plot(x0psd[0], y0psd[0], 'bx') ax.grid(True, alpha=0.3) ax.set_xlabel('x [mm]') ax.set_ylabel('y [mm]') plt.show()