diff --git a/fourierOptics/__pycache__/micado.cpython-38.pyc b/fourierOptics/__pycache__/micado.cpython-38.pyc index b9879b0..168b186 100644 Binary files a/fourierOptics/__pycache__/micado.cpython-38.pyc and b/fourierOptics/__pycache__/micado.cpython-38.pyc differ diff --git a/fourierOptics/__pycache__/micado.cpython-39.pyc b/fourierOptics/__pycache__/micado.cpython-39.pyc index 45c3a4c..360a98e 100644 Binary files a/fourierOptics/__pycache__/micado.cpython-39.pyc and b/fourierOptics/__pycache__/micado.cpython-39.pyc differ diff --git a/fourierOptics/__pycache__/misc.cpython-38.pyc b/fourierOptics/__pycache__/misc.cpython-38.pyc index d80488e..77eb07c 100644 Binary files a/fourierOptics/__pycache__/misc.cpython-38.pyc and b/fourierOptics/__pycache__/misc.cpython-38.pyc differ diff --git a/fourierOptics/__pycache__/misc.cpython-39.pyc b/fourierOptics/__pycache__/misc.cpython-39.pyc index efbc57b..dc429bf 100644 Binary files a/fourierOptics/__pycache__/misc.cpython-39.pyc and b/fourierOptics/__pycache__/misc.cpython-39.pyc differ diff --git a/fourierOptics/__pycache__/optics.cpython-38.pyc b/fourierOptics/__pycache__/optics.cpython-38.pyc index 15681f7..8cbf2f8 100644 Binary files a/fourierOptics/__pycache__/optics.cpython-38.pyc and b/fourierOptics/__pycache__/optics.cpython-38.pyc differ diff --git a/fourierOptics/__pycache__/optics.cpython-39.pyc b/fourierOptics/__pycache__/optics.cpython-39.pyc index 9299800..83ef3f5 100644 Binary files a/fourierOptics/__pycache__/optics.cpython-39.pyc and b/fourierOptics/__pycache__/optics.cpython-39.pyc differ diff --git a/fourierOptics/__pycache__/wavefrontErrors.cpython-38.pyc b/fourierOptics/__pycache__/wavefrontErrors.cpython-38.pyc index 93a2685..af65e31 100644 Binary files a/fourierOptics/__pycache__/wavefrontErrors.cpython-38.pyc and b/fourierOptics/__pycache__/wavefrontErrors.cpython-38.pyc differ diff --git a/fourierOptics/__pycache__/wavefrontErrors.cpython-39.pyc b/fourierOptics/__pycache__/wavefrontErrors.cpython-39.pyc index 5ee68f9..206504b 100644 Binary files a/fourierOptics/__pycache__/wavefrontErrors.cpython-39.pyc and b/fourierOptics/__pycache__/wavefrontErrors.cpython-39.pyc differ diff --git a/fourierOptics/micado.py b/fourierOptics/micado.py index 002af97..6cfff7a 100644 --- a/fourierOptics/micado.py +++ b/fourierOptics/micado.py @@ -244,7 +244,7 @@ class ParaxialMICADO: if print_progress: print('Propagating from ELT-IFP to MCD-COL') Ein = Eout if self.ELT_FP.wavefrontError: - Ein *= self.ELT_FP.phase_delay(x, y, l, ip_type='surface') + Ein *= self.ELT_FP.phase_delay(x, y, l, ip_type='phase') Eout, x, y = fou.fresnel_propagation(Ein, x, y, dz=self.ELT_COL.f, oversample=None, show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s], apply_interpolation=apply_ip[s], tx=tilt_x[s], ty=tilt_y[s]) diff --git a/fourierOptics/misc.py b/fourierOptics/misc.py index 8c211e8..47f5bbc 100644 --- a/fourierOptics/misc.py +++ b/fourierOptics/misc.py @@ -1,7 +1,10 @@ +''' +Summary +''' import numpy as np import pickle import matplotlib.pyplot as plt - +from numba import njit ''' Miscellaneous functions that do not belong specifically to any of the other files. @@ -13,15 +16,31 @@ import matplotlib.pyplot as plt def perform_zero_padding(A, xs, ys): ''' - Add zero padding around the array - - Input: - A - 2D-Input array - xs,ys - The new shape after padding. - - Output: - - Array with the new padded zeros - + Add zero padding around the array + + Input + ----- + A - 2D-Input array + xs,ys - The new shape after padding. + + Output + ------ + - Array with the new padded zeros + + Parameters + ---------- + A : TYPE + Description + xs : TYPE + Description + ys : TYPE + Description + + Returns + ------- + TYPE + Description + ''' assert A.shape[0] <= xs and A.shape[1] <= ys, 'New size is smaller than array size.' x_margin = int((xs - A.shape[0])/2) @@ -31,9 +50,21 @@ def perform_zero_padding(A, xs, ys): def find_nearest(array, value): ''' - For a value it finds the nearest element in an array - Returns the index of that element in the original - array. + For a value it finds the nearest element in an array + Returns the index of that element in the original + array. + + Parameters + ---------- + array : TYPE + Description + value : TYPE + Description + + Returns + ------- + TYPE + Description ''' array = np.asarray(array) diff = (np.abs(array - value)).argmin() @@ -43,14 +74,34 @@ def find_nearest(array, value): def find_maximum(array): ''' - Finds the index of the maximum of the given array + Finds the index of the maximum of the given array + + Parameters + ---------- + array : TYPE + Description + + Returns + ------- + TYPE + Description ''' idx = np.unravel_index(array.argmax(), array.shape) return idx def find_pixel_scale(x): - ''' - Calculates the pixel size from the supplied array. + ''' + Calculates the pixel size from the supplied array. + + Parameters + ---------- + x : TYPE + Description + + Returns + ------- + TYPE + Description ''' _dpix = np.mean(np.diff(x.real)) if _dpix < 1e-99: @@ -60,7 +111,21 @@ def find_pixel_scale(x): def select_subsection(x, y , outputCoordinates): ''' - Returns the indices of a subsection of an x,y meshgrid + Returns the indices of a subsection of an x,y meshgrid + + Parameters + ---------- + x : TYPE + Description + y : TYPE + Description + outputCoordinates : TYPE + Description + + Returns + ------- + TYPE + Description ''' select_1 = find_nearest(x, outputCoordinates[0])[1] select_2 = find_nearest(x, outputCoordinates[1])[1] @@ -71,10 +136,26 @@ def select_subsection(x, y , outputCoordinates): def perform_extrapolation(x, y, n, m): ''' - Performs linear extrapolation of a meshgrid to a new - shape n,m - - Returns the new extrapolated x,y meshgrid + Performs linear extrapolation of a meshgrid to a new + shape n,m + + Returns the new extrapolated x,y meshgrid + + Parameters + ---------- + x : TYPE + Description + y : TYPE + Description + n : TYPE + Description + m : TYPE + Description + + Returns + ------- + TYPE + Description ''' # Determine properties of x,y meshgrid @@ -92,36 +173,52 @@ def perform_extrapolation(x, y, n, m): def plot_results(E, x, y, logscale=False, fig_size=(8,3), grid=False, show=True, colormap='hot'): ''' - Quick function to plot complex valued image data, for example the - amplitude field of a PSF. The absolute value of the complex values - will be used. - If the amplitude values can be smaller than zero, use - plot_real_valued_results(). - - Parameters - ---------- - E - 2D Array - Amplitude or intensity of the field to be shown. - x - Array - X-coordinates of the field to be plotted - y - Array - Y-coordinates of the field to be plotted - logscale - Bool (default=False) - If True, imshow() will use a logarithmic color scale - fig_size - Tuple (default=(8,8)) - Size of the figure - grid - Bool (default=False) - If True, the figure will show gridlines. - show - Bool (default=True) - If True, the figure will be shown before continuing. - colormap - String (default='hot') - String, defines the colormap imshow() will use. - labeltext - String (default=None) - Defines the label to be shown along the colorbar - - Returns - ------- - _fig, _ax - Matplotlib figure and axis object. + Quick function to plot complex valued image data, for example the + amplitude field of a PSF. The absolute value of the complex values + will be used. + If the amplitude values can be smaller than zero, use + plot_real_valued_results(). + + Parameters + ---------- + E : TYPE + Description + x : TYPE + Description + y : TYPE + Description + logscale : bool, optional + Description + fig_size : tuple, optional + Description + grid : bool, optional + Description + show : bool, optional + Description + colormap : str, optional + Description + E - 2D Array + Amplitude or intensity of the field to be shown. + x - Array + X-coordinates of the field to be plotted + y - Array + Y-coordinates of the field to be plotted + logscale - Bool (default=False) + If True, imshow() will use a logarithmic color scale + fig_size - Tuple (default=(8,8)) + Size of the figure + grid - Bool (default=False) + If True, the figure will show gridlines. + show - Bool (default=True) + If True, the figure will be shown before continuing. + colormap - String (default='hot') + String, defines the colormap imshow() will use. + labeltext - String (default=None) + Defines the label to be shown along the colorbar + + Returns + ------- + _fig, _ax - Matplotlib figure and axis object. ''' x = x.real y = y.real @@ -155,33 +252,51 @@ def plot_results(E, x, y, logscale=False, fig_size=(8,3), grid=False, show=True, def plot_real_valued_results(E, x, y, logscale=False, fig_size=(8,8), grid=False, show=True, colormap='hot', labeltext=None): ''' - Quick function to plot real valued image data, for example the - intensity field of a PSF. - - Parameters - ---------- - E - 2D Array - Amplitude or intensity of the field to be shown. - x - Array - X-coordinates of the field to be plotted - y - Array - Y-coordinates of the field to be plotted - logscale - Bool (default=False) - If True, imshow() will use a logarithmic color scale - fig_size - Tuple (default=(8,8)) - Size of the figure - grid - Bool (default=False) - If True, the figure will show gridlines. - show - Bool (default=True) - If True, the figure will be shown before continuing. - colormap - String (default='hot') - String, defines the colormap imshow() will use. - labeltext - String (default=None) - Defines the label to be shown along the colorbar - - Returns - ------- - _fig, _ax - Matplotlib figure and axis object. + Quick function to plot real valued image data, for example the + intensity field of a PSF. + + Parameters + ---------- + E : TYPE + Description + x : TYPE + Description + y : TYPE + Description + logscale : bool, optional + Description + fig_size : tuple, optional + Description + grid : bool, optional + Description + show : bool, optional + Description + colormap : str, optional + Description + labeltext : None, optional + Description + E - 2D Array + Amplitude or intensity of the field to be shown. + x - Array + X-coordinates of the field to be plotted + y - Array + Y-coordinates of the field to be plotted + logscale - Bool (default=False) + If True, imshow() will use a logarithmic color scale + fig_size - Tuple (default=(8,8)) + Size of the figure + grid - Bool (default=False) + If True, the figure will show gridlines. + show - Bool (default=True) + If True, the figure will be shown before continuing. + colormap - String (default='hot') + String, defines the colormap imshow() will use. + labeltext - String (default=None) + Defines the label to be shown along the colorbar + + Returns + ------- + _fig, _ax - Matplotlib figure and axis object. ''' x = x.real @@ -208,8 +323,15 @@ def plot_real_valued_results(E, x, y, logscale=False, fig_size=(8,8), grid=False def save_variable(filename, variable): ''' - Dumps the given variable to a binary pickle file - for later use. + Dumps the given variable to a binary pickle file + for later use. + + Parameters + ---------- + filename : TYPE + Description + variable : TYPE + Description ''' with open(filename, 'wb') as f: pickle.dump(variable, f, protocol=-1 ) @@ -217,7 +339,17 @@ def save_variable(filename, variable): def load_variable(filename): ''' - Loads a binary pickle file for use. + Loads a binary pickle file for use. + + Parameters + ---------- + filename : TYPE + Description + + Returns + ------- + TYPE + Description ''' with open(filename, 'rb') as f: variable = pickle.load(f) @@ -226,26 +358,26 @@ def load_variable(filename): def geometric_distortions(xp, yp, x, y, magnification=1): ''' - Plots a grid showing the geometric distortions. Effects can be magnified. - - Parameters - ---------- - xp : array - Meshgrid array denoting the perturbed x-coordinates - yp : array - Meshgrid array denoting the perturbed y-coordinates - x : array - Meshgrid array denoting the unperturbed x-coordinates - y : array - Meshgrid array denoting the unperturbed y-coordinates - magnification : float (optional) - If given, amplifies the perturbation, so that it is more visible. - - Returns - ------- - fig, ax : matplotlib.pyplot object - Figure that can be shown using plt.show(), or modified by the user. - + Plots a grid showing the geometric distortions. Effects can be magnified. + + Parameters + ---------- + xp : array + Meshgrid array denoting the perturbed x-coordinates + yp : array + Meshgrid array denoting the perturbed y-coordinates + x : array + Meshgrid array denoting the unperturbed x-coordinates + y : array + Meshgrid array denoting the unperturbed y-coordinates + magnification : float (optional) + If given, amplifies the perturbation, so that it is more visible. + + Returns + ------- + fig, ax : matplotlib.pyplot object + Figure that can be shown using plt.show(), or modified by the user. + ''' dx = (xp - x) * magnification dy = (yp - y) * magnification @@ -267,53 +399,53 @@ def geometric_distortions(xp, yp, x, y, magnification=1): def get_prism_rtm(x_in, y_in, theta_x_in, theta_y_in, thickness, a1x, a2x, a1y, a2y, n_glass): ''' - Finds the ray transfer matrix for a prism with the given parameters. - It does not assume small angles. It does not take into account any - (internal) reflections. The results were confirmed using CODEV. - It is assumed that the prism is located inside a vacuum (i.e. n=1 - outside the prism). - - Parameters - ---------- - x_in : float - Input x-coordinate at the point in the plane where the optical - axis reaches the front of the prism - y_in : float - Input y-coordinate at the point in the plane where the optical - axis reaches the front of the prism - theta_x_in : float - Input angle about the y-axis (so it affects the x-coordinates) - theta_y_in : float - Input angle about the x-axis (so it affects the y-coordinates) - thickness : float - Thickness of the prism along the optical axis - a1x : float - Rotation angle of the first prism surface about the y-axis, - so that the x-coordinates are affected. - a2x : float - Rotation angle of the second prism surface about the y-axis, - so that the x-coordinates are affected. - a1y : float - Rotation angle of the first prism surface about the x-axis, - so that the y-coordinates are affected. - a2y : float - Rotation angle of the second prism surface about the x-axis, - so that the y-coordinates are affected. - n_glass : float - Refractive index of the prism glass - - Returns - ------- - x_out : float - Output x-coordinate at the point in the plane where the optical - axis reaches the back of the prism - y_out : float - Output y-coordinate at the point in the plane where the optical - axis reaches the back of the prism - theta_x_out : float - Output angle about the y-axis (so it affects the x-coordinates) - theta_y_out : float - Output angle about the x-axis (so it affects the y-coordinates) + Finds the ray transfer matrix for a prism with the given parameters. + It does not assume small angles. It does not take into account any + (internal) reflections. The results were confirmed using CODEV. + It is assumed that the prism is located inside a vacuum (i.e. n=1 + outside the prism). + + Parameters + ---------- + x_in : float + Input x-coordinate at the point in the plane where the optical + axis reaches the front of the prism + y_in : float + Input y-coordinate at the point in the plane where the optical + axis reaches the front of the prism + theta_x_in : float + Input angle about the y-axis (so it affects the x-coordinates) + theta_y_in : float + Input angle about the x-axis (so it affects the y-coordinates) + thickness : float + Thickness of the prism along the optical axis + a1x : float + Rotation angle of the first prism surface about the y-axis, + so that the x-coordinates are affected. + a2x : float + Rotation angle of the second prism surface about the y-axis, + so that the x-coordinates are affected. + a1y : float + Rotation angle of the first prism surface about the x-axis, + so that the y-coordinates are affected. + a2y : float + Rotation angle of the second prism surface about the x-axis, + so that the y-coordinates are affected. + n_glass : float + Refractive index of the prism glass + + Returns + ------- + x_out : float + Output x-coordinate at the point in the plane where the optical + axis reaches the back of the prism + y_out : float + Output y-coordinate at the point in the plane where the optical + axis reaches the back of the prism + theta_x_out : float + Output angle about the y-axis (so it affects the x-coordinates) + theta_y_out : float + Output angle about the x-axis (so it affects the y-coordinates) ''' theta_x_out = -1 * np.arcsin( -n_glass * np.sin( np.arcsin( 1/n_glass * np.sin(theta_x_in - a1x) ) + a1x - a2x ) ) + a2x theta_y_out = -1 * np.arcsin( -n_glass * np.sin( np.arcsin( 1/n_glass * np.sin(theta_y_in - a1y) ) + a1y - a2y ) ) + a2y @@ -335,40 +467,40 @@ def get_prism_rtm(x_in, y_in, theta_x_in, theta_y_in, thickness, a1x, a2x, a1y, def get_paraxial_lens_rtm(x_in, y_in, theta_x_in, theta_y_in, f_lens_x, f_lens_y): ''' - Finds the ray transfer matrix of a paraxial lens, without small angle - approximations. - - Parameters - ---------- - x_in : float - Input x-coordinate at the point in the plane where the optical - axis reaches the front of the prism - y_in : float - Input y-coordinate at the point in the plane where the optical - axis reaches the front of the prism - theta_x_in : float - Input angle about the y-axis (so it affects the x-coordinates) - theta_y_in : float - Input angle about the x-axis (so it affects the y-coordinates) - f_lens_x : float - Focal length of the lens for the points with a nonzero x-coordinate - value. - f_lens_y : float - Focal length of the lens for the points with a nonzero y-coordinate - value. - - Returns - ------- - x_out : float - Output x-coordinate at the point in the plane where the optical - axis reaches the back of the prism - y_out : float - Output y-coordinate at the point in the plane where the optical - axis reaches the back of the prism - theta_x_out : float - Output angle about the y-axis (so it affects the x-coordinates) - theta_y_out : float - Output angle about the x-axis (so it affects the y-coordinates) + Finds the ray transfer matrix of a paraxial lens, without small angle + approximations. + + Parameters + ---------- + x_in : float + Input x-coordinate at the point in the plane where the optical + axis reaches the front of the prism + y_in : float + Input y-coordinate at the point in the plane where the optical + axis reaches the front of the prism + theta_x_in : float + Input angle about the y-axis (so it affects the x-coordinates) + theta_y_in : float + Input angle about the x-axis (so it affects the y-coordinates) + f_lens_x : float + Focal length of the lens for the points with a nonzero x-coordinate + value. + f_lens_y : float + Focal length of the lens for the points with a nonzero y-coordinate + value. + + Returns + ------- + x_out : float + Output x-coordinate at the point in the plane where the optical + axis reaches the back of the prism + y_out : float + Output y-coordinate at the point in the plane where the optical + axis reaches the back of the prism + theta_x_out : float + Output angle about the y-axis (so it affects the x-coordinates) + theta_y_out : float + Output angle about the x-axis (so it affects the y-coordinates) ''' x_out = x_in y_out = y_in @@ -379,39 +511,94 @@ def get_paraxial_lens_rtm(x_in, y_in, theta_x_in, theta_y_in, f_lens_x, f_lens_y def get_free_space_propagation_rtm(x_in, y_in, theta_x_in, theta_y_in, distance): ''' - Finds the ray transfer matrix of a free space propagation, without small - angle approximations. - - Parameters - ---------- - x_in : float - Input x-coordinate at the point in the plane where the optical - axis reaches the front of the prism - y_in : float - Input y-coordinate at the point in the plane where the optical - axis reaches the front of the prism - theta_x_in : float - Input angle about the y-axis (so it affects the x-coordinates) - theta_y_in : float - Input angle about the x-axis (so it affects the y-coordinates) - distance : float - Distance to travel through free space. - - Returns - ------- - x_out : float - Output x-coordinate at the point in the plane where the optical - axis reaches the back of the prism - y_out : float - Output y-coordinate at the point in the plane where the optical - axis reaches the back of the prism - theta_x_out : float - Output angle about the y-axis (so it affects the x-coordinates) - theta_y_out : float - Output angle about the x-axis (so it affects the y-coordinates) + Finds the ray transfer matrix of a free space propagation, without small + angle approximations. + + Parameters + ---------- + x_in : float + Input x-coordinate at the point in the plane where the optical + axis reaches the front of the prism + y_in : float + Input y-coordinate at the point in the plane where the optical + axis reaches the front of the prism + theta_x_in : float + Input angle about the y-axis (so it affects the x-coordinates) + theta_y_in : float + Input angle about the x-axis (so it affects the y-coordinates) + distance : float + Distance to travel through free space. + + Returns + ------- + x_out : float + Output x-coordinate at the point in the plane where the optical + axis reaches the back of the prism + y_out : float + Output y-coordinate at the point in the plane where the optical + axis reaches the back of the prism + theta_x_out : float + Output angle about the y-axis (so it affects the x-coordinates) + theta_y_out : float + Output angle about the x-axis (so it affects the y-coordinates) ''' x_out = x_in + distance * np.tan(theta_x_in) y_out = y_in + distance * np.tan(theta_y_in) theta_x_out = theta_x_in theta_y_out = theta_y_in return x_out, y_out, theta_x_out, theta_y_out + + +def fourier_series(x, y, c, Px=30, Py=30, xs=10, ys=10): + ''' + Calculates the 2D continous Fourier Series for a set of + coefficients. This allows for interpolation on a random + phase field without artefacts (in theory...). + It is however, really quite slow for large arrays. + + Parameters + ---------- + x : meshgrid array + x-coordinates of the 2D array of interest + y : meshgrid array + y-coordinates of the 2D array of interest + c : 2D array + coefficients of the Fourier Series + Px : int, optional + Period length in x + Py : int, optional + Period length in y + xs : int, optional + Shift of the symmetry point x + ys : int, optional + Shift of the symmetry point in y + + Returns + ------- + array + 2D array of the Fourier Series + ''' + # Shift the coordinates to remove the symmetry point + # from the output field. + xpos = x + np.sqrt(2)*xs + ypos = y + np.sqrt(2)*ys + normf = np.sum(c) + + # Preparation for the calculation, to prevent doing + # the same thing over and over again. + xt = 2*np.pi*xpos/Px + yt = 2*np.pi*ypos/Py + # Doing the calculation using numba + fout = numba_fourier_series(c, xt, yt) + return fout/normf + +@njit +def numba_fourier_series(c, xt, yt): + cs0, cs1 = c.shape + fout = np.zeros(xt.shape, dtype=np.float64) + for j in range(cs1): + for i in range(cs0): + v = c[i,j] * np.exp(1j * (i * xt + j * yt)) + fout += v.real + return fout + diff --git a/fourierOptics/optics.py b/fourierOptics/optics.py index a0de6d4..4605209 100644 --- a/fourierOptics/optics.py +++ b/fourierOptics/optics.py @@ -475,16 +475,15 @@ class OpticalSurface: wfeMap = interp2d(ximg[0,:], yimg[:,0], self.sag, kind=kind, **kwargs) wavefrontErrorMap = wfeMap(x[0,:].real, y[:,0].real) _t_l *= np.exp(1j * _k * wavefrontErrorMap ) - elif ip_type == 'phase': + elif ip_type in ['phase', 'fourier_series']: try: xy = deepcopy((x,y)) xy_img = deepcopy(self.sag_coordinates) wavefrontErrorMap = self.wfe.generate_height_map_window( - xy, xy_img, ip_type='phase', kind=kind) + xy, xy_img, ip_type=ip_type, kind=kind) _t_l *= np.exp(1j * _k * wavefrontErrorMap ) except AttributeError: print('No WavefrontError object detected. Continuing without!') - if self.mask is not None: _t_l *= self.mask return _t_l diff --git a/fourierOptics/wavefrontErrors.py b/fourierOptics/wavefrontErrors.py index d1b76b8..bcfffe5 100644 --- a/fourierOptics/wavefrontErrors.py +++ b/fourierOptics/wavefrontErrors.py @@ -4,6 +4,7 @@ import types from math import ceil from scipy.special import factorial, j1 from scipy import interpolate +from .misc import fourier_series class WavefrontError: ''' @@ -277,7 +278,7 @@ class WavefrontError: return _phaseAngle - def generate_height_map(self, s0, p, rho_c, addToData=False): + def generate_height_map(self, s0, p, rho_c, addToData=False,): ''' This function generates a surface height distribution based on a Lorentzian PSD profile. Following the paper by E. Sidick (2009). @@ -289,9 +290,14 @@ class WavefrontError: _dx = 1 / self.m _dy = 1 / self.n _A = (self.m * _dx) * (self.n * _dy) - _phase = np.random.uniform(0,1,(self.n,self.m)) - self.phase = _phase - + try: + _phase = self.phase + print('Phase Attribute present!') + except AttributeError: + + _phase = np.random.uniform(0,1,(self.n,self.m)) + self.phase = _phase + # Calculating the PSD and phase _PSD, _rho_mn_list = self.power_spectral_density(s0, p, rho_c, _dx, _dy, _A) _phiMap = self.phase_map(_phase) @@ -327,12 +333,15 @@ class WavefrontError: xy - tuple containing the x and y coordinates of the small window xy_img - tuple containing the x and y coordinates of the full plane to which the wfe is applied. - ip_type - 'phase' or 'surface' + ip_type - 'phase' or 'surface' or 'fourier_series' Defines what method to use to obtain the interpolated surface height map. 'surface' does a normal interpolation on the height map. This does not conserve the WFE rms value. Then, 'phase' interpolates on the phase information in the PSD algorithm. This does conserve the WFE rms within the window. + The last option, 'fourier_series' is very slow, but does provide + an analytical and continous solution that does not rely on + interpolation. ''' # normalise the coordinates to -1 to 1 ximg, yimg = xy_img @@ -376,14 +385,47 @@ class WavefrontError: ptvw = np.abs(np.max(_hMapw.real)) + np.abs(np.min(_hMapw.real)) elif ip_type == 'surface': - print('Hello there!') assert hasattr(self, 'data'), 'No wfe data available to interpolate on!' try: _hMapw = self.wfe_ip(x[0,:], y[:,0]) except AttributeError: - print('Hello there 2!') self.wfe_ip = interpolate.interp2d(ximg[0,:], yimg[:,0], self.data, kind=kind) _hMapw = self.wfe_ip(x[0,:], y[:,0]) + + elif ip_type == 'fourier_series': + # Check if the coefficients are present + assert hasattr(self, 'coefficients'), 'No coefficients for the Fourier Series present!' + + # Reinterpreting the coordinates, because they do not need to be scaled here + x, y = xy + x, y = x.real, y.real + + # Calculation of the necessary constants + _dx = 1 / self.m + _dy = 1 / self.n + _A = (self.m * _dx) * (self.n * _dy) + + # Generating the array at coordinates x, y + print('Applying the Fourier series. This might take a while...') + _phase_window = fourier_series(x, y, self.coefficients) + print('Fourier series calculated!') + + # Generating the surface + _PSDw, _rho_mn_listw = self.PSD2D, self.rho_mn_list + _phiMap_window = self.phase_map(_phase_window) + _PSDw = np.fft.ifftshift(_PSDw) + + # Doing the fourier transform of the window + _inw = pyfftw.empty_aligned((self.n,self.m), dtype='complex128') + _inw[:] = np.exp( 1j * _phiMap_window) * np.sqrt(_A * np.abs(_PSDw)) + _outw = pyfftw.empty_aligned((self.n,self.m), dtype='complex128') + _FFTW_Objw = pyfftw.FFTW( _inw, _outw, axes=(0,1), direction='FFTW_BACKWARD', threads=6, flags=('FFTW_ESTIMATE',)) + _FFTW_Objw() + + _hMapw = 1 / ( _dx * _dy ) * _outw + rmsw = np.sqrt(np.mean(np.asarray(_hMapw.real)**2)) + ptvw = np.abs(np.max(_hMapw.real)) + np.abs(np.min(_hMapw.real)) + if not addToData: self.data_window = _hMapw.real diff --git a/scripts/SPIE/Thesis/Adding a WFE in the focal plane.ipynb b/scripts/SPIE/Thesis/Adding a WFE in the focal plane.ipynb index a959bee..5ee447f 100644 --- a/scripts/SPIE/Thesis/Adding a WFE in the focal plane.ipynb +++ b/scripts/SPIE/Thesis/Adding a WFE in the focal plane.ipynb @@ -6634,12 +6634,15 @@ "name": "stdout", "output_type": "stream", "text": [ + "-2.8244144799708755\n", "\n", - "[-2.62118705+0.41501946j -0.70850924-0.41518566j 3.24396455+5.77869097j\n", - " -9.24802109+7.67145042j -5.3838523 -5.17726043j 0.58191832+2.7330733j\n", - " -6.01265766+1.18573479j -2.82441448+0.j -6.01265766-1.18573479j\n", - " 0.58191832-2.7330733j -5.3838523 +5.17726043j]\n", + "[ 0.92804617-0.14694j 0.25085172+0.14699884j -1.14854409-2.04597838j\n", + " 3.27431443-2.71612062j 1.90618351+1.83303848j -0.20603149-0.96766014j\n", + " 2.12881562-0.41981614j 1. -0.j 2.12881562+0.41981614j\n", + " -0.20603149+0.96766014j 1.90618351-1.83303848j]\n", + "-2.8244144799708755\n", "\n", + "-2.8244144799708755\n", "\n" ] }, @@ -7629,7 +7632,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -7674,10 +7677,11 @@ " xs - shift of the full field x-coordinates. To remove symmetry around 0.\n", " '''\n", " xpos = x - xs\n", + " norm_factor = np.sum(c)\n", + " print(norm_factor)\n", " f = np.asarray([c[i]*np.exp(1j*i* 2*np.pi/P * xpos) for i in range(len(c))]).sum(axis=0)\n", - "# f = c*np.exp(1j*i*2*np.pi/P)\n", " print(type(f[0]))\n", - " return f\n", + " return f/norm_factor\n", "\n", "# 1D Fourier series test\n", "np.random.seed(1)\n", @@ -7718,23 +7722,24 @@ "'''" ] }, - { - "cell_type": "markdown", - "id": "94b0d227", - "metadata": {}, - "source": [] - }, { "cell_type": "code", "execution_count": 3, - "id": "b7a59755", + "id": "80f810ef", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Elapsed time: 0.0037 s\n" + "-2.8244144799708755\n", + "\n", + "[ 0.92804617-0.14694j 0.73482068-0.1828004j 0.52660171+0.0308276j\n", + " -0.15380869+0.14587545j -1.03250078-0.52772232j -1.14854409-2.04597838j\n", + " 0.15220504-3.4237516j 2.32344441-3.40161665j 3.87639295-1.69179571j\n", + " 3.65598886+0.5920894j 1.90618351+1.83303848j]\n", + "-2.8244144799708755\n", + "\n" ] }, { @@ -8723,7 +8728,1046 @@ { "data": { "text/html": [ - "" + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# 1D Fourier series test\n", + "np.random.seed(1)\n", + "coeffs = np.random.uniform(-1,1,100)\n", + "# coeffs[0] = 1\n", + "# coeffs[1] = 1\n", + "# coeffs[2] = 1\n", + "\n", + "xx0 = np.linspace(-10,10,11)\n", + "yy0 = fourierSeries1D(xx0, coeffs, P=1000)\n", + "print(yy0)\n", + "xx1 = np.linspace(-5, 15, 11)\n", + "yy1 = fourierSeries1D(xx1,coeffs, P=1000)\n", + "# xx2 = np.linspace(-30, 45, 21)\n", + "# yy2 = fourierSeries1D(xx2,coeffs, P=1000)\n", + "plt.figure()\n", + "plt.plot(xx0,yy0.real, 'g', marker='x', linestyle='--', zorder=2)\n", + "plt.plot(xx1,yy1.real, 'b-', zorder=2)\n", + "\n", + "# plt.plot(xx2,yy2.real, 'r-', zorder=3, linewidth=1)\n", + "plt.grid(alpha=0.3, zorder=0)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "94b0d227", + "metadata": {}, + "source": [] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "b7a59755", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Elapsed time: 0.0068 s\n" + ] + }, + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "/* global mpl */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function () {\n", + " if (typeof WebSocket !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof MozWebSocket !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert(\n", + " 'Your browser does not have WebSocket support. ' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.'\n", + " );\n", + " }\n", + "};\n", + "\n", + "mpl.figure = function (figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = this.ws.binaryType !== undefined;\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById('mpl-warnings');\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent =\n", + " 'This browser does not support binary websocket messages. ' +\n", + " 'Performance may be slow.';\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = document.createElement('div');\n", + " this.root.setAttribute('style', 'display: inline-block');\n", + " this._root_extra_style(this.root);\n", + "\n", + " parent_element.appendChild(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message('supports_binary', { value: fig.supports_binary });\n", + " fig.send_message('send_image_mode', {});\n", + " if (fig.ratio !== 1) {\n", + " fig.send_message('set_dpi_ratio', { dpi_ratio: fig.ratio });\n", + " }\n", + " fig.send_message('refresh', {});\n", + " };\n", + "\n", + " this.imageObj.onload = function () {\n", + " if (fig.image_mode === 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function () {\n", + " fig.ws.close();\n", + " };\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "};\n", + "\n", + "mpl.figure.prototype._init_header = function () {\n", + " var titlebar = document.createElement('div');\n", + " titlebar.classList =\n", + " 'ui-dialog-titlebar ui-widget-header ui-corner-all ui-helper-clearfix';\n", + " var titletext = document.createElement('div');\n", + " titletext.classList = 'ui-dialog-title';\n", + " titletext.setAttribute(\n", + " 'style',\n", + " 'width: 100%; text-align: center; padding: 3px;'\n", + " );\n", + " titlebar.appendChild(titletext);\n", + " this.root.appendChild(titlebar);\n", + " this.header = titletext;\n", + "};\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function (_canvas_div) {};\n", + "\n", + "mpl.figure.prototype._root_extra_style = function (_canvas_div) {};\n", + "\n", + "mpl.figure.prototype._init_canvas = function () {\n", + " var fig = this;\n", + "\n", + " var canvas_div = (this.canvas_div = document.createElement('div'));\n", + " canvas_div.setAttribute(\n", + " 'style',\n", + " 'border: 1px solid #ddd;' +\n", + " 'box-sizing: content-box;' +\n", + " 'clear: both;' +\n", + " 'min-height: 1px;' +\n", + " 'min-width: 1px;' +\n", + " 'outline: 0;' +\n", + " 'overflow: hidden;' +\n", + " 'position: relative;' +\n", + " 'resize: both;'\n", + " );\n", + "\n", + " function on_keyboard_event_closure(name) {\n", + " return function (event) {\n", + " return fig.key_event(event, name);\n", + " };\n", + " }\n", + "\n", + " canvas_div.addEventListener(\n", + " 'keydown',\n", + " on_keyboard_event_closure('key_press')\n", + " );\n", + " canvas_div.addEventListener(\n", + " 'keyup',\n", + " on_keyboard_event_closure('key_release')\n", + " );\n", + "\n", + " this._canvas_extra_style(canvas_div);\n", + " this.root.appendChild(canvas_div);\n", + "\n", + " var canvas = (this.canvas = document.createElement('canvas'));\n", + " canvas.classList.add('mpl-canvas');\n", + " canvas.setAttribute('style', 'box-sizing: content-box;');\n", + "\n", + " this.context = canvas.getContext('2d');\n", + "\n", + " var backingStore =\n", + " this.context.backingStorePixelRatio ||\n", + " this.context.webkitBackingStorePixelRatio ||\n", + " this.context.mozBackingStorePixelRatio ||\n", + " this.context.msBackingStorePixelRatio ||\n", + " this.context.oBackingStorePixelRatio ||\n", + " this.context.backingStorePixelRatio ||\n", + " 1;\n", + "\n", + " this.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband_canvas = (this.rubberband_canvas = document.createElement(\n", + " 'canvas'\n", + " ));\n", + " rubberband_canvas.setAttribute(\n", + " 'style',\n", + " 'box-sizing: content-box; position: absolute; left: 0; top: 0; z-index: 1;'\n", + " );\n", + "\n", + " // Apply a ponyfill if ResizeObserver is not implemented by browser.\n", + " if (this.ResizeObserver === undefined) {\n", + " if (window.ResizeObserver !== undefined) {\n", + " this.ResizeObserver = window.ResizeObserver;\n", + " } else {\n", + " var obs = _JSXTOOLS_RESIZE_OBSERVER({});\n", + " this.ResizeObserver = obs.ResizeObserver;\n", + " }\n", + " }\n", + "\n", + " this.resizeObserverInstance = new this.ResizeObserver(function (entries) {\n", + " var nentries = entries.length;\n", + " for (var i = 0; i < nentries; i++) {\n", + " var entry = entries[i];\n", + " var width, height;\n", + " if (entry.contentBoxSize) {\n", + " if (entry.contentBoxSize instanceof Array) {\n", + " // Chrome 84 implements new version of spec.\n", + " width = entry.contentBoxSize[0].inlineSize;\n", + " height = entry.contentBoxSize[0].blockSize;\n", + " } else {\n", + " // Firefox implements old version of spec.\n", + " width = entry.contentBoxSize.inlineSize;\n", + " height = entry.contentBoxSize.blockSize;\n", + " }\n", + " } else {\n", + " // Chrome <84 implements even older version of spec.\n", + " width = entry.contentRect.width;\n", + " height = entry.contentRect.height;\n", + " }\n", + "\n", + " // Keep the size of the canvas and rubber band canvas in sync with\n", + " // the canvas container.\n", + " if (entry.devicePixelContentBoxSize) {\n", + " // Chrome 84 implements new version of spec.\n", + " canvas.setAttribute(\n", + " 'width',\n", + " entry.devicePixelContentBoxSize[0].inlineSize\n", + " );\n", + " canvas.setAttribute(\n", + " 'height',\n", + " entry.devicePixelContentBoxSize[0].blockSize\n", + " );\n", + " } else {\n", + " canvas.setAttribute('width', width * fig.ratio);\n", + " canvas.setAttribute('height', height * fig.ratio);\n", + " }\n", + " canvas.setAttribute(\n", + " 'style',\n", + " 'width: ' + width + 'px; height: ' + height + 'px;'\n", + " );\n", + "\n", + " rubberband_canvas.setAttribute('width', width);\n", + " rubberband_canvas.setAttribute('height', height);\n", + "\n", + " // And update the size in Python. We ignore the initial 0/0 size\n", + " // that occurs as the element is placed into the DOM, which should\n", + " // otherwise not happen due to the minimum size styling.\n", + " if (fig.ws.readyState == 1 && width != 0 && height != 0) {\n", + " fig.request_resize(width, height);\n", + " }\n", + " }\n", + " });\n", + " this.resizeObserverInstance.observe(canvas_div);\n", + "\n", + " function on_mouse_event_closure(name) {\n", + " return function (event) {\n", + " return fig.mouse_event(event, name);\n", + " };\n", + " }\n", + "\n", + " rubberband_canvas.addEventListener(\n", + " 'mousedown',\n", + " on_mouse_event_closure('button_press')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseup',\n", + " on_mouse_event_closure('button_release')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'dblclick',\n", + " on_mouse_event_closure('dblclick')\n", + " );\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband_canvas.addEventListener(\n", + " 'mousemove',\n", + " on_mouse_event_closure('motion_notify')\n", + " );\n", + "\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseenter',\n", + " on_mouse_event_closure('figure_enter')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseleave',\n", + " on_mouse_event_closure('figure_leave')\n", + " );\n", + "\n", + " canvas_div.addEventListener('wheel', function (event) {\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " on_mouse_event_closure('scroll')(event);\n", + " });\n", + "\n", + " canvas_div.appendChild(canvas);\n", + " canvas_div.appendChild(rubberband_canvas);\n", + "\n", + " this.rubberband_context = rubberband_canvas.getContext('2d');\n", + " this.rubberband_context.strokeStyle = '#000000';\n", + "\n", + " this._resize_canvas = function (width, height, forward) {\n", + " if (forward) {\n", + " canvas_div.style.width = width + 'px';\n", + " canvas_div.style.height = height + 'px';\n", + " }\n", + " };\n", + "\n", + " // Disable right mouse context menu.\n", + " this.rubberband_canvas.addEventListener('contextmenu', function (_e) {\n", + " event.preventDefault();\n", + " return false;\n", + " });\n", + "\n", + " function set_focus() {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "};\n", + "\n", + "mpl.figure.prototype._init_toolbar = function () {\n", + " var fig = this;\n", + "\n", + " var toolbar = document.createElement('div');\n", + " toolbar.classList = 'mpl-toolbar';\n", + " this.root.appendChild(toolbar);\n", + "\n", + " function on_click_closure(name) {\n", + " return function (_event) {\n", + " return fig.toolbar_button_onclick(name);\n", + " };\n", + " }\n", + "\n", + " function on_mouseover_closure(tooltip) {\n", + " return function (event) {\n", + " if (!event.currentTarget.disabled) {\n", + " return fig.toolbar_button_onmouseover(tooltip);\n", + " }\n", + " };\n", + " }\n", + "\n", + " fig.buttons = {};\n", + " var buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'mpl-button-group';\n", + " for (var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " /* Instead of a spacer, we start a new button group. */\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + " buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'mpl-button-group';\n", + " continue;\n", + " }\n", + "\n", + " var button = (fig.buttons[name] = document.createElement('button'));\n", + " button.classList = 'mpl-widget';\n", + " button.setAttribute('role', 'button');\n", + " button.setAttribute('aria-disabled', 'false');\n", + " button.addEventListener('click', on_click_closure(method_name));\n", + " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", + "\n", + " var icon_img = document.createElement('img');\n", + " icon_img.src = '_images/' + image + '.png';\n", + " icon_img.srcset = '_images/' + image + '_large.png 2x';\n", + " icon_img.alt = tooltip;\n", + " button.appendChild(icon_img);\n", + "\n", + " buttonGroup.appendChild(button);\n", + " }\n", + "\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + "\n", + " var fmt_picker = document.createElement('select');\n", + " fmt_picker.classList = 'mpl-widget';\n", + " toolbar.appendChild(fmt_picker);\n", + " this.format_dropdown = fmt_picker;\n", + "\n", + " for (var ind in mpl.extensions) {\n", + " var fmt = mpl.extensions[ind];\n", + " var option = document.createElement('option');\n", + " option.selected = fmt === mpl.default_extension;\n", + " option.innerHTML = fmt;\n", + " fmt_picker.appendChild(option);\n", + " }\n", + "\n", + " var status_bar = document.createElement('span');\n", + " status_bar.classList = 'mpl-message';\n", + " toolbar.appendChild(status_bar);\n", + " this.message = status_bar;\n", + "};\n", + "\n", + "mpl.figure.prototype.request_resize = function (x_pixels, y_pixels) {\n", + " // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n", + " // which will in turn request a refresh of the image.\n", + " this.send_message('resize', { width: x_pixels, height: y_pixels });\n", + "};\n", + "\n", + "mpl.figure.prototype.send_message = function (type, properties) {\n", + " properties['type'] = type;\n", + " properties['figure_id'] = this.id;\n", + " this.ws.send(JSON.stringify(properties));\n", + "};\n", + "\n", + "mpl.figure.prototype.send_draw_message = function () {\n", + " if (!this.waiting) {\n", + " this.waiting = true;\n", + " this.ws.send(JSON.stringify({ type: 'draw', figure_id: this.id }));\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", + " var format_dropdown = fig.format_dropdown;\n", + " var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n", + " fig.ondownload(fig, format);\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_resize = function (fig, msg) {\n", + " var size = msg['size'];\n", + " if (size[0] !== fig.canvas.width || size[1] !== fig.canvas.height) {\n", + " fig._resize_canvas(size[0], size[1], msg['forward']);\n", + " fig.send_message('refresh', {});\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_rubberband = function (fig, msg) {\n", + " var x0 = msg['x0'] / fig.ratio;\n", + " var y0 = (fig.canvas.height - msg['y0']) / fig.ratio;\n", + " var x1 = msg['x1'] / fig.ratio;\n", + " var y1 = (fig.canvas.height - msg['y1']) / fig.ratio;\n", + " x0 = Math.floor(x0) + 0.5;\n", + " y0 = Math.floor(y0) + 0.5;\n", + " x1 = Math.floor(x1) + 0.5;\n", + " y1 = Math.floor(y1) + 0.5;\n", + " var min_x = Math.min(x0, x1);\n", + " var min_y = Math.min(y0, y1);\n", + " var width = Math.abs(x1 - x0);\n", + " var height = Math.abs(y1 - y0);\n", + "\n", + " fig.rubberband_context.clearRect(\n", + " 0,\n", + " 0,\n", + " fig.canvas.width / fig.ratio,\n", + " fig.canvas.height / fig.ratio\n", + " );\n", + "\n", + " fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_figure_label = function (fig, msg) {\n", + " // Updates the figure title.\n", + " fig.header.textContent = msg['label'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_cursor = function (fig, msg) {\n", + " var cursor = msg['cursor'];\n", + " switch (cursor) {\n", + " case 0:\n", + " cursor = 'pointer';\n", + " break;\n", + " case 1:\n", + " cursor = 'default';\n", + " break;\n", + " case 2:\n", + " cursor = 'crosshair';\n", + " break;\n", + " case 3:\n", + " cursor = 'move';\n", + " break;\n", + " }\n", + " fig.rubberband_canvas.style.cursor = cursor;\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_message = function (fig, msg) {\n", + " fig.message.textContent = msg['message'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_draw = function (fig, _msg) {\n", + " // Request the server to send over a new figure.\n", + " fig.send_draw_message();\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_image_mode = function (fig, msg) {\n", + " fig.image_mode = msg['mode'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_history_buttons = function (fig, msg) {\n", + " for (var key in msg) {\n", + " if (!(key in fig.buttons)) {\n", + " continue;\n", + " }\n", + " fig.buttons[key].disabled = !msg[key];\n", + " fig.buttons[key].setAttribute('aria-disabled', !msg[key]);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_navigate_mode = function (fig, msg) {\n", + " if (msg['mode'] === 'PAN') {\n", + " fig.buttons['Pan'].classList.add('active');\n", + " fig.buttons['Zoom'].classList.remove('active');\n", + " } else if (msg['mode'] === 'ZOOM') {\n", + " fig.buttons['Pan'].classList.remove('active');\n", + " fig.buttons['Zoom'].classList.add('active');\n", + " } else {\n", + " fig.buttons['Pan'].classList.remove('active');\n", + " fig.buttons['Zoom'].classList.remove('active');\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.updated_canvas_event = function () {\n", + " // Called whenever the canvas gets updated.\n", + " this.send_message('ack', {});\n", + "};\n", + "\n", + "// A function to construct a web socket function for onmessage handling.\n", + "// Called in the figure constructor.\n", + "mpl.figure.prototype._make_on_message_function = function (fig) {\n", + " return function socket_on_message(evt) {\n", + " if (evt.data instanceof Blob) {\n", + " var img = evt.data;\n", + " if (img.type !== 'image/png') {\n", + " /* FIXME: We get \"Resource interpreted as Image but\n", + " * transferred with MIME type text/plain:\" errors on\n", + " * Chrome. But how to set the MIME type? It doesn't seem\n", + " * to be part of the websocket stream */\n", + " img.type = 'image/png';\n", + " }\n", + "\n", + " /* Free the memory for the previous frames */\n", + " if (fig.imageObj.src) {\n", + " (window.URL || window.webkitURL).revokeObjectURL(\n", + " fig.imageObj.src\n", + " );\n", + " }\n", + "\n", + " fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n", + " img\n", + " );\n", + " fig.updated_canvas_event();\n", + " fig.waiting = false;\n", + " return;\n", + " } else if (\n", + " typeof evt.data === 'string' &&\n", + " evt.data.slice(0, 21) === 'data:image/png;base64'\n", + " ) {\n", + " fig.imageObj.src = evt.data;\n", + " fig.updated_canvas_event();\n", + " fig.waiting = false;\n", + " return;\n", + " }\n", + "\n", + " var msg = JSON.parse(evt.data);\n", + " var msg_type = msg['type'];\n", + "\n", + " // Call the \"handle_{type}\" callback, which takes\n", + " // the figure and JSON message as its only arguments.\n", + " try {\n", + " var callback = fig['handle_' + msg_type];\n", + " } catch (e) {\n", + " console.log(\n", + " \"No handler for the '\" + msg_type + \"' message type: \",\n", + " msg\n", + " );\n", + " return;\n", + " }\n", + "\n", + " if (callback) {\n", + " try {\n", + " // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n", + " callback(fig, msg);\n", + " } catch (e) {\n", + " console.log(\n", + " \"Exception inside the 'handler_\" + msg_type + \"' callback:\",\n", + " e,\n", + " e.stack,\n", + " msg\n", + " );\n", + " }\n", + " }\n", + " };\n", + "};\n", + "\n", + "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n", + "mpl.findpos = function (e) {\n", + " //this section is from http://www.quirksmode.org/js/events_properties.html\n", + " var targ;\n", + " if (!e) {\n", + " e = window.event;\n", + " }\n", + " if (e.target) {\n", + " targ = e.target;\n", + " } else if (e.srcElement) {\n", + " targ = e.srcElement;\n", + " }\n", + " if (targ.nodeType === 3) {\n", + " // defeat Safari bug\n", + " targ = targ.parentNode;\n", + " }\n", + "\n", + " // pageX,Y are the mouse positions relative to the document\n", + " var boundingRect = targ.getBoundingClientRect();\n", + " var x = e.pageX - (boundingRect.left + document.body.scrollLeft);\n", + " var y = e.pageY - (boundingRect.top + document.body.scrollTop);\n", + "\n", + " return { x: x, y: y };\n", + "};\n", + "\n", + "/*\n", + " * return a copy of an object with only non-object keys\n", + " * we need this to avoid circular references\n", + " * http://stackoverflow.com/a/24161582/3208463\n", + " */\n", + "function simpleKeys(original) {\n", + " return Object.keys(original).reduce(function (obj, key) {\n", + " if (typeof original[key] !== 'object') {\n", + " obj[key] = original[key];\n", + " }\n", + " return obj;\n", + " }, {});\n", + "}\n", + "\n", + "mpl.figure.prototype.mouse_event = function (event, name) {\n", + " var canvas_pos = mpl.findpos(event);\n", + "\n", + " if (name === 'button_press') {\n", + " this.canvas.focus();\n", + " this.canvas_div.focus();\n", + " }\n", + "\n", + " var x = canvas_pos.x * this.ratio;\n", + " var y = canvas_pos.y * this.ratio;\n", + "\n", + " this.send_message(name, {\n", + " x: x,\n", + " y: y,\n", + " button: event.button,\n", + " step: event.step,\n", + " guiEvent: simpleKeys(event),\n", + " });\n", + "\n", + " /* This prevents the web browser from automatically changing to\n", + " * the text insertion cursor when the button is pressed. We want\n", + " * to control all of the cursor setting manually through the\n", + " * 'cursor' event from matplotlib */\n", + " event.preventDefault();\n", + " return false;\n", + "};\n", + "\n", + "mpl.figure.prototype._key_event_extra = function (_event, _name) {\n", + " // Handle any extra behaviour associated with a key event\n", + "};\n", + "\n", + "mpl.figure.prototype.key_event = function (event, name) {\n", + " // Prevent repeat events\n", + " if (name === 'key_press') {\n", + " if (event.key === this._key) {\n", + " return;\n", + " } else {\n", + " this._key = event.key;\n", + " }\n", + " }\n", + " if (name === 'key_release') {\n", + " this._key = null;\n", + " }\n", + "\n", + " var value = '';\n", + " if (event.ctrlKey && event.key !== 'Control') {\n", + " value += 'ctrl+';\n", + " }\n", + " else if (event.altKey && event.key !== 'Alt') {\n", + " value += 'alt+';\n", + " }\n", + " else if (event.shiftKey && event.key !== 'Shift') {\n", + " value += 'shift+';\n", + " }\n", + "\n", + " value += 'k' + event.key;\n", + "\n", + " this._key_event_extra(event, name);\n", + "\n", + " this.send_message(name, { key: value, guiEvent: simpleKeys(event) });\n", + " return false;\n", + "};\n", + "\n", + "mpl.figure.prototype.toolbar_button_onclick = function (name) {\n", + " if (name === 'download') {\n", + " this.handle_save(this, null);\n", + " } else {\n", + " this.send_message('toolbar_button', { name: name });\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.toolbar_button_onmouseover = function (tooltip) {\n", + " this.message.textContent = tooltip;\n", + "};\n", + "\n", + "///////////////// REMAINING CONTENT GENERATED BY embed_js.py /////////////////\n", + "// prettier-ignore\n", + "var _JSXTOOLS_RESIZE_OBSERVER=function(A){var t,i=new WeakMap,n=new WeakMap,a=new WeakMap,r=new WeakMap,o=new Set;function s(e){if(!(this instanceof s))throw new TypeError(\"Constructor requires 'new' operator\");i.set(this,e)}function h(){throw new TypeError(\"Function is not a constructor\")}function c(e,t,i,n){e=0 in arguments?Number(arguments[0]):0,t=1 in arguments?Number(arguments[1]):0,i=2 in arguments?Number(arguments[2]):0,n=3 in arguments?Number(arguments[3]):0,this.right=(this.x=this.left=e)+(this.width=i),this.bottom=(this.y=this.top=t)+(this.height=n),Object.freeze(this)}function d(){t=requestAnimationFrame(d);var s=new WeakMap,p=new Set;o.forEach((function(t){r.get(t).forEach((function(i){var r=t instanceof window.SVGElement,o=a.get(t),d=r?0:parseFloat(o.paddingTop),f=r?0:parseFloat(o.paddingRight),l=r?0:parseFloat(o.paddingBottom),u=r?0:parseFloat(o.paddingLeft),g=r?0:parseFloat(o.borderTopWidth),m=r?0:parseFloat(o.borderRightWidth),w=r?0:parseFloat(o.borderBottomWidth),b=u+f,F=d+l,v=(r?0:parseFloat(o.borderLeftWidth))+m,W=g+w,y=r?0:t.offsetHeight-W-t.clientHeight,E=r?0:t.offsetWidth-v-t.clientWidth,R=b+v,z=F+W,M=r?t.width:parseFloat(o.width)-R-E,O=r?t.height:parseFloat(o.height)-z-y;if(n.has(t)){var k=n.get(t);if(k[0]===M&&k[1]===O)return}n.set(t,[M,O]);var S=Object.create(h.prototype);S.target=t,S.contentRect=new c(u,d,M,O),s.has(i)||(s.set(i,[]),p.add(i)),s.get(i).push(S)}))})),p.forEach((function(e){i.get(e).call(e,s.get(e),e)}))}return s.prototype.observe=function(i){if(i instanceof window.Element){r.has(i)||(r.set(i,new Set),o.add(i),a.set(i,window.getComputedStyle(i)));var n=r.get(i);n.has(this)||n.add(this),cancelAnimationFrame(t),t=requestAnimationFrame(d)}},s.prototype.unobserve=function(i){if(i instanceof window.Element&&r.has(i)){var n=r.get(i);n.has(this)&&(n.delete(this),n.size||(r.delete(i),o.delete(i))),n.size||r.delete(i),o.size||cancelAnimationFrame(t)}},A.DOMRectReadOnly=c,A.ResizeObserver=s,A.ResizeObserverEntry=h,A}; // eslint-disable-line\n", + "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Left button pans, Right button zooms\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n", + "\n", + "mpl.extensions = [\"eps\", \"jpeg\", \"pgf\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n", + "\n", + "mpl.default_extension = \"png\";/* global mpl */\n", + "\n", + "var comm_websocket_adapter = function (comm) {\n", + " // Create a \"websocket\"-like object which calls the given IPython comm\n", + " // object with the appropriate methods. Currently this is a non binary\n", + " // socket, so there is still some room for performance tuning.\n", + " var ws = {};\n", + "\n", + " ws.binaryType = comm.kernel.ws.binaryType;\n", + " ws.readyState = comm.kernel.ws.readyState;\n", + " function updateReadyState(_event) {\n", + " if (comm.kernel.ws) {\n", + " ws.readyState = comm.kernel.ws.readyState;\n", + " } else {\n", + " ws.readyState = 3; // Closed state.\n", + " }\n", + " }\n", + " comm.kernel.ws.addEventListener('open', updateReadyState);\n", + " comm.kernel.ws.addEventListener('close', updateReadyState);\n", + " comm.kernel.ws.addEventListener('error', updateReadyState);\n", + "\n", + " ws.close = function () {\n", + " comm.close();\n", + " };\n", + " ws.send = function (m) {\n", + " //console.log('sending', m);\n", + " comm.send(m);\n", + " };\n", + " // Register the callback with on_msg.\n", + " comm.on_msg(function (msg) {\n", + " //console.log('receiving', msg['content']['data'], msg);\n", + " var data = msg['content']['data'];\n", + " if (data['blob'] !== undefined) {\n", + " data = {\n", + " data: new Blob(msg['buffers'], { type: data['blob'] }),\n", + " };\n", + " }\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", + " ws.onmessage(data);\n", + " });\n", + " return ws;\n", + "};\n", + "\n", + "mpl.mpl_figure_comm = function (comm, msg) {\n", + " // This is the function which gets called when the mpl process\n", + " // starts-up an IPython Comm through the \"matplotlib\" channel.\n", + "\n", + " var id = msg.content.data.id;\n", + " // Get hold of the div created by the display call when the Comm\n", + " // socket was opened in Python.\n", + " var element = document.getElementById(id);\n", + " var ws_proxy = comm_websocket_adapter(comm);\n", + "\n", + " function ondownload(figure, _format) {\n", + " window.open(figure.canvas.toDataURL());\n", + " }\n", + "\n", + " var fig = new mpl.figure(id, ws_proxy, ondownload, element);\n", + "\n", + " // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n", + " // web socket which is closed, not our websocket->open comm proxy.\n", + " ws_proxy.onopen();\n", + "\n", + " fig.parent_element = element;\n", + " fig.cell_info = mpl.find_output_cell(\"
\");\n", + " if (!fig.cell_info) {\n", + " console.error('Failed to find cell for figure', id, fig);\n", + " return;\n", + " }\n", + " fig.cell_info[0].output_area.element.on(\n", + " 'cleared',\n", + " { fig: fig },\n", + " fig._remove_fig_handler\n", + " );\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_close = function (fig, msg) {\n", + " var width = fig.canvas.width / fig.ratio;\n", + " fig.cell_info[0].output_area.element.off(\n", + " 'cleared',\n", + " fig._remove_fig_handler\n", + " );\n", + " fig.resizeObserverInstance.unobserve(fig.canvas_div);\n", + "\n", + " // Update the output cell to use the data from the current canvas.\n", + " fig.push_to_output();\n", + " var dataURL = fig.canvas.toDataURL();\n", + " // Re-enable the keyboard manager in IPython - without this line, in FF,\n", + " // the notebook keyboard shortcuts fail.\n", + " IPython.keyboard_manager.enable();\n", + " fig.parent_element.innerHTML =\n", + " '';\n", + " fig.close_ws(fig, msg);\n", + "};\n", + "\n", + "mpl.figure.prototype.close_ws = function (fig, msg) {\n", + " fig.send_message('closing', msg);\n", + " // fig.ws.close()\n", + "};\n", + "\n", + "mpl.figure.prototype.push_to_output = function (_remove_interactive) {\n", + " // Turn the data on the canvas into data in the output cell.\n", + " var width = this.canvas.width / this.ratio;\n", + " var dataURL = this.canvas.toDataURL();\n", + " this.cell_info[1]['text/html'] =\n", + " '';\n", + "};\n", + "\n", + "mpl.figure.prototype.updated_canvas_event = function () {\n", + " // Tell IPython that the notebook contents must change.\n", + " IPython.notebook.set_dirty(true);\n", + " this.send_message('ack', {});\n", + " var fig = this;\n", + " // Wait a second, then push the new image to the DOM so\n", + " // that it is saved nicely (might be nice to debounce this).\n", + " setTimeout(function () {\n", + " fig.push_to_output();\n", + " }, 1000);\n", + "};\n", + "\n", + "mpl.figure.prototype._init_toolbar = function () {\n", + " var fig = this;\n", + "\n", + " var toolbar = document.createElement('div');\n", + " toolbar.classList = 'btn-toolbar';\n", + " this.root.appendChild(toolbar);\n", + "\n", + " function on_click_closure(name) {\n", + " return function (_event) {\n", + " return fig.toolbar_button_onclick(name);\n", + " };\n", + " }\n", + "\n", + " function on_mouseover_closure(tooltip) {\n", + " return function (event) {\n", + " if (!event.currentTarget.disabled) {\n", + " return fig.toolbar_button_onmouseover(tooltip);\n", + " }\n", + " };\n", + " }\n", + "\n", + " fig.buttons = {};\n", + " var buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'btn-group';\n", + " var button;\n", + " for (var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " /* Instead of a spacer, we start a new button group. */\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + " buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'btn-group';\n", + " continue;\n", + " }\n", + "\n", + " button = fig.buttons[name] = document.createElement('button');\n", + " button.classList = 'btn btn-default';\n", + " button.href = '#';\n", + " button.title = name;\n", + " button.innerHTML = '';\n", + " button.addEventListener('click', on_click_closure(method_name));\n", + " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", + " buttonGroup.appendChild(button);\n", + " }\n", + "\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = document.createElement('span');\n", + " status_bar.classList = 'mpl-message pull-right';\n", + " toolbar.appendChild(status_bar);\n", + " this.message = status_bar;\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = document.createElement('div');\n", + " buttongrp.classList = 'btn-group inline pull-right';\n", + " button = document.createElement('button');\n", + " button.classList = 'btn btn-mini btn-primary';\n", + " button.href = '#';\n", + " button.title = 'Stop Interaction';\n", + " button.innerHTML = '';\n", + " button.addEventListener('click', function (_evt) {\n", + " fig.handle_close(fig, {});\n", + " });\n", + " button.addEventListener(\n", + " 'mouseover',\n", + " on_mouseover_closure('Stop Interaction')\n", + " );\n", + " buttongrp.appendChild(button);\n", + " var titlebar = this.root.querySelector('.ui-dialog-titlebar');\n", + " titlebar.insertBefore(buttongrp, titlebar.firstChild);\n", + "};\n", + "\n", + "mpl.figure.prototype._remove_fig_handler = function (event) {\n", + " var fig = event.data.fig;\n", + " if (event.target !== this) {\n", + " // Ignore bubbled events from children.\n", + " return;\n", + " }\n", + " fig.close_ws(fig, {});\n", + "};\n", + "\n", + "mpl.figure.prototype._root_extra_style = function (el) {\n", + " el.style.boxSizing = 'content-box'; // override notebook setting of border-box.\n", + "};\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function (el) {\n", + " // this is important to make the div 'focusable\n", + " el.setAttribute('tabindex', 0);\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " } else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype._key_event_extra = function (event, _name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager) {\n", + " manager = IPython.keyboard_manager;\n", + " }\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which === 13) {\n", + " this.canvas_div.blur();\n", + " // select the cell after this one\n", + " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", + " IPython.notebook.select(index + 1);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", + " fig.ondownload(fig, null);\n", + "};\n", + "\n", + "mpl.find_output_cell = function (html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i = 0; i < ncells; i++) {\n", + " var cell = cells[i];\n", + " if (cell.cell_type === 'code') {\n", + " for (var j = 0; j < cell.output_area.outputs.length; j++) {\n", + " var data = cell.output_area.outputs[j];\n", + " if (data.data) {\n", + " // IPython >= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] === html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "};\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel !== null) {\n", + " IPython.notebook.kernel.comm_manager.register_target(\n", + " 'matplotlib',\n", + " mpl.mpl_figure_comm\n", + " );\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" ], "text/plain": [ "" @@ -8774,7 +9818,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 7, "id": "2d7b525b", "metadata": {}, "outputs": [ @@ -8783,7 +9827,7 @@ "output_type": "stream", "text": [ "Old one: 0.0000 s\n", - "New one: 1.8907 s\n" + "New one: 0.4503 s\n" ] }, { @@ -9772,7 +10816,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -9841,6 +10885,7 @@ " xpos = x + np.sqrt(2)*xs\n", " ypos = y + np.sqrt(2)*ys\n", " cs0, cs1 = c.shape\n", + " normf = np.sum(c)\n", " xt = 2*np.pi*xpos/Px\n", " yt = 2*np.pi*ypos/Py\n", " c_ij = c.flatten()\n", @@ -9848,7 +10893,7 @@ " I = i.flatten()\n", " J = j.flatten()\n", " fout = f(c, i, j, xt, yt)\n", - " return fout\n", + " return fout/normf\n", "\n", "@njit\n", "def f(c, i, j, xt, yt):\n", @@ -9896,10 +10941,10 @@ "\n", "# 2D Fourier series test\n", "np.random.seed(2183)\n", - "coeffs = np.random.uniform(-1,1,(101,101))\n", + "coeffs = np.random.uniform(-1,1,(11,11))\n", "# coeffs = np.ones((101,101))\n", - "x0 = np.linspace(-50,50,1001)\n", - "y0 = np.linspace(-50,50,1001)\n", + "x0 = np.linspace(-50,50,101)\n", + "y0 = np.linspace(-50,50,101)\n", "xx, yy = np.meshgrid(x0,y0)\n", "s0 = time.time()\n", "# zz = fourierSeries2Dopt2(xx, yy, coeffs, )#Px=100, Py=100, xs=0, ys=0)\n", @@ -9935,7 +10980,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 8, "id": "2ad234e1", "metadata": {}, "outputs": [ @@ -10925,7 +11970,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -11972,62 +13017,6 @@ "plt.show()" ] }, - { - "cell_type": "code", - "execution_count": 98, - "id": "cd1b454d", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[[ 0 1 2 ... 97 98 99]\n", - " [ 0 1 2 ... 97 98 99]\n", - " [ 0 1 2 ... 97 98 99]\n", - " ...\n", - " [ 0 1 2 ... 97 98 99]\n", - " [ 0 1 2 ... 97 98 99]\n", - " [ 0 1 2 ... 97 98 99]] [[ 0 0 0 ... 0 0 0]\n", - " [ 1 1 1 ... 1 1 1]\n", - " [ 2 2 2 ... 2 2 2]\n", - " ...\n", - " [97 97 97 ... 97 97 97]\n", - " [98 98 98 ... 98 98 98]\n", - " [99 99 99 ... 99 99 99]]\n", - "[[ 0 0 0 ... 0 0 0]\n", - " [ 1 1 1 ... 1 1 1]\n", - " [ 2 2 2 ... 2 2 2]\n", - " ...\n", - " [97 97 97 ... 97 97 97]\n", - " [98 98 98 ... 98 98 98]\n", - " [99 99 99 ... 99 99 99]] [[ 0 1 2 ... 97 98 99]\n", - " [ 0 1 2 ... 97 98 99]\n", - " [ 0 1 2 ... 97 98 99]\n", - " ...\n", - " [ 0 1 2 ... 97 98 99]\n", - " [ 0 1 2 ... 97 98 99]\n", - " [ 0 1 2 ... 97 98 99]]\n" - ] - } - ], - "source": [ - "x, y = np.meshgrid(range(100), range(100))\n", - "print(x,y)\n", - "\n", - "def meshgrid(x, y):\n", - " xx = np.empty(shape=(len(x), len(y)), dtype=int)\n", - " yy = np.empty(shape=(len(x), len(y)), dtype=int)\n", - " for j in y:\n", - " for k in x:\n", - " xx[j,k] = k # change to x[k] if indexing xy\n", - " yy[j,k] = j # change to y[j] if indexing xy\n", - " return yy, xx\n", - "\n", - "x2, y2 = meshgrid(range(100), range(100))\n", - "print(x2, y2)" - ] - }, { "cell_type": "code", "execution_count": null, diff --git a/scripts/SPIE/Thesis/generate_fp_wfe_phase_array.py b/scripts/SPIE/Thesis/generate_fp_wfe_phase_array.py new file mode 100644 index 0000000..c42f628 --- /dev/null +++ b/scripts/SPIE/Thesis/generate_fp_wfe_phase_array.py @@ -0,0 +1,86 @@ +import numpy as np +import matplotlib.pyplot as plt +import sys +# sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +# sys.path.insert(0, '/home/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Documents/Published/FourierOptics GitHub release/') +from fourierOptics import * +from pdb import set_trace + +''' + @Author: Joost van den Born + @Date: 14 July, 2022 + @Description: I use this script to create a large array (1001x1001 and 3001x3001) + that can be used for the phase screen in the WFE production. + Because it takes very long to generate these at this high resolution, + we will need to save them to a pickle file. + + What is saved: + The array of coefficients + The array of the full phase field + + To Do: + - Normalize/Scale the amplitude in an appropriate manner. +''' + + +if __name__ == '__main__': + # plt.rc('text', usetex=True) + # plt.rc('font', family='serif') + import time + # Setup of simulation + l = 1.49e-3 + # mcd = ParaxialMICADO(n=1001, m=1001, dpix=160) + mcd = ParaxialMICADO(n=3001, m=3001, dpix=50) + # mcd.plot_ray_trace() + wfe = WavefrontError(mcd.n, mcd.m) + + + # ------------------------------------------------------- + ####### GENERATING THE PHASE WITH A FOURIER SERIES ###### + # ------------------------------------------------------- + # Get x3 and y3 (the coordinates to apply the ADC WFE to): + # The focal plane area is about 200x200 mm, so we apply the wavefront error to an area + # of 200x200 mm. + xfp, yfp = np.meshgrid(np.linspace(-150, 150, mcd.n), np.linspace(-150, 150, mcd.m)) + # np.random.seed(2183) + # coeff = np.random.uniform(-1,1,(int(mcd.n/3),int(mcd.m/3))) + np.random.seed(2183) + coeff = np.random.uniform(-1,1,(101,101)) + + if False: + print('Going in...') + s1 = time.time() + fs_array = fourier_series(xfp, yfp, coeff, Px=300, Py=300, xs=150, ys=150) + s2 = time.time() + print('Finished!') + print('That took me {:.4f} s!'.format(s2-s1) ) + save_variable('wfe_fp_phase_array_coefficients_{}_{}.pkl'.format(mcd.n, mcd.m), coeff) + save_variable('wfe_fp_phase_array_{}_{}.pkl'.format(mcd.n, mcd.m), fs_array) + + # ------------------------------------------------------- + # ------------------------------------------------------- + # ------------------------------------------------------- + + if True: + # coeff = load_variable('wfe_fp_phase_array_coefficients_{}_{}.pkl'.format(mcd.n, mcd.m)) + fs_array = load_variable('wfe_fp_phase_array_{}_{}.pkl'.format(mcd.n, mcd.m)) + + wfe.phase = fs_array + surf = wfe.generate_height_map(50e-6, 2, 1) + rms = np.sqrt(np.mean(surf**2)) + print('RMS: {}'.format(rms)) + + + plt.figure() + plt.imshow(fs_array, cmap='jet', extent=[xfp.min(), xfp.max(), yfp.min(), yfp.max()], origin='lower') + plt.colorbar() + plt.show() + + + plt.figure() + plt.imshow(surf, cmap='jet', extent=[xfp.min(), xfp.max(), yfp.min(), yfp.max()], origin='lower') + plt.colorbar() + plt.show() + + diff --git a/scripts/SPIE/Thesis/geometric_distortions_fp_wfe.py b/scripts/SPIE/Thesis/geometric_distortions_fp_wfe.py index ca1db18..9beb73f 100644 --- a/scripts/SPIE/Thesis/geometric_distortions_fp_wfe.py +++ b/scripts/SPIE/Thesis/geometric_distortions_fp_wfe.py @@ -16,14 +16,14 @@ from pdb import set_trace Script made for the SPIE manuscript. ''' -def plot_distortions(): +def plot_distortions(xe, ye, xf, yf): # Coordinates without WFEs - xe = load_variable('../data/GD_IFP_obscured_PSF_xcentroid_expected_zern.pkl') - ye = load_variable('../data/GD_IFP_obscured_PSF_ycentroid_expected_zern.pkl') + # xe = load_variable('../data/GD_IFP_zoom93_obscured_PSF_xcentroid_expected_zern_psd.pkl') + # ye = load_variable('../data/GD_IFP_zoom93_obscured_PSF_ycentroid_expected_zern_psd.pkl') - # Coordinates with WFEs - xf = load_variable('../data/GD_IFP_obscured_PSF_xcentroid_found_with_wfe_airy_zern.pkl') - yf = load_variable('../data/GD_IFP_obscured_PSF_ycentroid_found_with_wfe_airy_zern.pkl') + # # Coordinates with WFEs + # xf = load_variable('../data/GD_IFP_zoom93_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + # yf = load_variable('../data/GD_IFP_zoom93_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') mcd = ParaxialMICADO(n=3001, m=3001, dpix=50, l=1.49e-3) @@ -53,8 +53,8 @@ def plot_distortions(): ax.set_xlabel('x (mm)') ax.set_ylabel('y (mm)') fig.colorbar(im, ax=ax, label='PSF centroid shift ($\\mu$as)') - plt.savefig('../figures/UPDATED_GD_focal_plane_zern.pdf', bbox_inches='tight', transparent=True) - plt.savefig('../figures/UPDATED_GD_focal_plane_zern.png', bbox_inches='tight', transparent=True) + # plt.savefig('../figures/UPDATED_GD_zoom93_focal_plane_zern.pdf', bbox_inches='tight', transparent=True) + # plt.savefig('../figures/UPDATED_GD_zoom93_focal_plane_zern.png', bbox_inches='tight', transparent=True) plt.tight_layout() plt.show() @@ -62,11 +62,10 @@ def plot_distortions(): if __name__ == '__main__': # plt.rc('text', usetex=True) # plt.rc('font', family='serif') - # Setup of simulation l = 1.49e-3 mcd = ParaxialMICADO(n=3001, m=3001, dpix=50) - # mcd.plot_ray_trace() + print(mcd.platescale) wfe = WavefrontError(mcd.n, mcd.m) # Get x3 and y3 (the coordinates to apply the ADC WFE to): @@ -80,34 +79,91 @@ if __name__ == '__main__': print('\nSetting the wavefront errors...') wfe_fp = WavefrontError(mcd.n,mcd.m) wfe_rms = 50e-6 - #PSD Wavefront error - # surface_map_fp = wfe_fp.generate_height_map(wfe_rms, 2, 1) - # plot_real_valued_results(surface_map_fp*1e6, xfp, yfp, logscale=False, fig_size=(8,6.4), grid=False, show=True, - # colormap='jet', labeltext='Surface height (nm)') - # save_variable('../data/MCD_IFP_wfe_object_PSD_50nm.pkl', wfe_fp) - # wfe_fp = load_variable('../data/MCD_IFP_wfe_object_PSD_50nm.pkl') - # surface_map_fp = wfe_fp.data - # Zernike wavefront error - # We re0se the ADC_surface classes because it can be applicable to any surface, if used properly - ifp_wfe = ADC_surfaces(xfp, yfp, clearAperture=300) - ifp_wfe.generate_ADC_wfe(zernRMS=wfe_rms, psdRMS=0.) - mask = np.sqrt(ifp_wfe.xz**2 + ifp_wfe.yz**2) > 1 - - # Plot the wfe - # surf_ifp = ifp_wfe.latest_surface.copy() - # surf_ifp[mask]=np.nan - - # surface_map_fp = surf_ifp - # wfe_fp.data = surf_ifp - wfe_fp = load_variable('../data/MCD_IFP_wfe_object_zern_50nm.pkl') - surface_map_fp = wfe_fp.data - # save_variable('../data/MCD_IFP_wfe_object_zern_50nm.pkl', wfe_fp) + psdWFE = True + zernikeWFE = False + + if psdWFE: + #PSD Wavefront error + # surface_map_fp = wfe_fp.generate_height_map(wfe_rms, 2, 1) + # plot_real_valued_results(surface_map_fp*1e6, xfp, yfp, logscale=False, fig_size=(8,6.4), grid=False, show=True, + # colormap='jet', labeltext='Surface height (nm)') + # save_variable('../data/MCD_IFP_wfe_object_PSD_50nm.pkl', wfe_fp) + # wfe_fp = load_variable('../data/MCD_IFP_wfe_object_PSD_50nm.pkl') + wfe_fp = load_variable('../data/GD Focal Plane WFEs/MCD_IFP_wfe_object_PSD_50nm.pkl') + surface_map_fp = wfe_fp.data + + # from copy import deepcopy + # xy_img = deepcopy((xfp, yfp)) + # x, y = np.meshgrid(np.linspace(-10, 10, mcd.n), np.linspace(-10, 10, mcd.m)) + # xy = deepcopy((x,y)) + # zoom_map = wfe_fp.generate_height_map_window( + # xy, xy_img, ip_type='phase') + + + # plt.figure() + # plt.grid() + # plt.plot(1e3*x[1500,:], 1e6*zoom_map[1500,:], '-') + # # plt.imshow(1e6*zoom_map, cmap='jet', extent=[x.min(), x.max(), y.min(), y.max()], + # # origin='lower', vmin=-150, vmax=150) + # plt.xlabel('x (μm)') + # plt.ylabel('surface heigth (nm)') + # plt.show() + + if zernikeWFE: + # Zernike wavefront error + # We reuse the ADC_surface classes because it can be applicable to any surface, if used properly + ifp_wfe = ADC_surfaces(xfp, yfp, clearAperture=300) + ifp_wfe.generate_ADC_wfe(zernRMS=wfe_rms, psdRMS=0.) + mask = np.sqrt(ifp_wfe.xz**2 + ifp_wfe.yz**2) > 1 + + # Plot the wfe + # surf_ifp = ifp_wfe.latest_surface.copy() + # surf_ifp[mask]=np.nan + + # surface_map_fp = surf_ifp + # wfe_fp.data = surf_ifp + # wfe_fp = load_variable('../data/GD Focal Plane WFEs/MCD_IFP_wfe_object_zern_50nm.pkl') + surface_map_fp = wfe_fp.data + # save_variable('../data/MCD_IFP_wfe_object_zern_50nm.pkl', wfe_fp) # figp1, axp1 = plot_real_valued_results(surface_map_fp*1e6, xfp, yfp, - # colormap='jet', labeltext='Surface height (nm)', fig_size=None, show=True) - plt.show() + # colormap='jet', labeltext='Surface height (nm)', fig_size=None, show=True) + # plt.show() + # plt.close() + # List of performed simulations + # Defined by PSD + ## GD_IFP : -1:1 + ## GD_IFP_zoom : -1/11 :1/11 + ## GD_IFP_zoom2 : -.5/11:.5/11 + ## GD_IFP_zoom3 : -.1/11:.1/11 + ## GD_IFP_zoom7 : -.02/11:.02/11, surface interpolation + ## GD_IFP_zoom8 : -.05/11:.05/11, surface interpolation + ## GD_IFP_zoom9 : -.01/11:.01/11, surface interpolation + ## GD_IFP_zoom91 : -.01/11:.01/11 but hx and hy shifted to an area near the corner of the image, where distortion is supposed to be large + ## GD_IFP_zoom92 : -.01/11:.01/11, phase interpolation + ## GD_IFP_zoom93 : -.004/11:.004/11, phase interpolation + # (implemented)IDEA: Automate it, let it run all simulations for the lower left corner. Take multiples of 11 for the field points + # i.e. 1/11, 1/11**2, 1/11**3, 1/11**4 + # GD_IFP_zoom94 : - 1/11:1/11, phase interpolation, with shift to left bottom corner + # GD_IFP_zoom95 : - 1/11**2:1/11**2, phase interpolation, with shift to left bottom corner + # GD_IFP_zoom96 : - 1/11**3:1/11**3, phase interpolation, with shift to left bottom corner + # GD_IFP_zoom97 : - 1/11:1/11, phase interpolation, center of the field + # GD_IFP_zoom98 : - 1/11**2:1/11**2, phase interpolation, center of the field + # GD_IFP_zoom99 : - 1/11**3:1/11**3, phase interpolation, center of the field + # GD_IFP_zoom100 : - 1/11**4:1/11**4, phase interpolation, center of the field + # GD_IFP_zoom101 : - 1/11**4:1/11**4, phase interpolation, center of the field, with shift to left bottom corner + + # Defined by zernikes + ## GD_IFP_obscured_PSF_xcentroid_expected_zern : 50 nm rms 19 zernike modes : -1 : 1 + ## To do; check PSF shape at center field with zern wfe + ## To do; check interpolated WFE at any zoomed in field point. + ## GD_IFP_zoom4 : -.02/11:.02/11 (actually zernike WFE!), surface interpolation --> CONTINUOUS :D + ## GD_IFP_zoom5 : -.05/11:.05/11 (actually zernike WFE!), surface interpolation --> CONTINUOUS :D + ## GD_IFP_zoom6 : -.05/11:.05/11, surface interpolation, but hx and hy shifted to an area near the corner of the image (actually zern!) + + Elist, Eplist = [], [] xlist, ylist = [], [] xplist, yplist = [], [] @@ -117,96 +173,99 @@ if __name__ == '__main__': # Set field points hxList, hyList = np.mgrid[-1:1:11j,-1:1:11j]-(1/300) # hxList, hyList = np.mgrid[-.1/11:.1/11:11j,-.1/11:.1/11:11j]-(1/300) # Zoomed in to central area + # hxList, hyList = np.mgrid[-.004/11:.004/11:11j,-.004/11:.004/11:11j]-(1/300) # Zoomed in to central area hxList = hxList.flatten() hyList = hyList.flatten() - # List of performed simulations - # Defined by PSD - ## GD_IFP : -1:1 - ## GD_IFP_zoom : -1/11 :1/11 - ## GD_IFP_zoom2 : -.5/11:.5/11 - ## GD_IFP_zoom3 : -.1/11:.1/11 + # Zooming in around -84mm by -84mm on the image plane + ihx = hxList[108] + ihy = hyList[108] - # Defined by zernikes - ## GD_IFP_obscured_PSF_xcentroid_expected_zern : 50 nm rms 19 zernike modes : -1 : 1 - ## To do; check PSF shape at center field with zern wfe - ## To do; check interpolated WFE at any zoomed in field point. + for p in [4]: + x_expected_list, y_expected_list = [], [] + x_centroid_list, y_centroid_list = [], [] - # Propagate through the field positions without any WFEs - for i in range(len(hxList)): - hx = hxList[i] - hy = hyList[i] - print('hx: {} - hy: {}'.format(hx,hy) ) + hxList, hyList = np.mgrid[ihx-1/11**p:ihx+1/11**p:11j, + ihy-1/11**p:ihy+1/11**p:11j]-(1/300) + hxList = hxList.flatten() + hyList = hyList.flatten() - # Initializing MICADO paraxial model - mcd = ParaxialMICADO(n=3001, m=3001, dpix=50, entrance_pupil='obscured') + # Propagate through the field positions without any WFEs + for i in range(len(hxList)): + print('Field point #: {}/{}'.format(i+1, len(hxList))) + hx = hxList[i] + hy = hyList[i] + print('hx: {} - hy: {}'.format(hx,hy) ) - # Propagate through the MICADO system - E, x, y = mcd.propagate(l=l,surface=-1, zoom=False, show=False, - print_progress=False, hx=hx, hy=hy) + # Initializing MICADO paraxial model + mcd = ParaxialMICADO(n=3001, m=3001, dpix=50, entrance_pupil='obscured') - # Find approximation for the centroid and amplitude as a first - # guess for the Obscured Airy pattern fit. - amplitude_guess = np.max(np.abs(E))**2 - xbc, ybc, _, _ = mcd.centroid_shift(method='Barycenter') - print('Barycenter guess: ', xbc, ybc) + # Propagate through the MICADO system + E, x, y = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=hx, hy=hy) - # Do the Obscured Airy Pattern Fit for the unperturbed PSF - xe, ye, _, _ = mcd.centroid_shift(method='Obscured Airy', initialGuess=(amplitude_guess, xbc, ybc, 0.0, 0.03)) - print('Airy result: ', xe, ye) - x_expected_list.append(xe) - y_expected_list.append(ye) + # Find approximation for the centroid and amplitude as a first + # guess for the Obscured Airy pattern fit. + amplitude_guess = np.max(np.abs(E))**2 + xbc, ybc, _, _ = mcd.centroid_shift(method='Barycenter') + print('Barycenter guess: ', xbc, ybc) - + # Do the Obscured Airy Pattern Fit for the unperturbed PSF + xe, ye, _, _ = mcd.centroid_shift(method='Obscured Airy', initialGuess=(amplitude_guess, xbc, ybc, 0.0, 0.03)) + print('Airy result: ', xe, ye) + x_expected_list.append(xe) + y_expected_list.append(ye) - # Load in the WFEs and set them to the corresponding ADC surfaces - print('\nSetting the wavefront errors...') - # null_wfe = np.zeros((3001,3001)) - mcd.ELT_FP.add_wavefront_error(surface_map_fp, xfp, yfp, wfe=wfe_fp) - - - # Propagate through the MICADO system - print('Propagating with Wavefront Errors applied...') - Ep, xp, yp = mcd.propagate(l=l,surface=-1, zoom=False, show=False, - print_progress=False, hx=hx, hy=hy) - # Do the Obscured Airy Pattern Fit for the perturbed PSF - xf, yf, _, _ = mcd.centroid_shift(method='Obscured Airy', initialGuess=(amplitude_guess, xe, ye, 0.0, 0.03)) - - # Resulting centroid shift due to the WFE - x_centroid_list.append(xf) - y_centroid_list.append(yf) - R = np.sqrt( (xf-xe)**2 + (yf-ye)**2 ) * mcd.platescale * 1e6 - print('xf: {}, yf: {}, shift (uas): {}'.format(xf,yf, R)) - - # Clear the memory by emptying the large arrays - mcd = None - wfe = None - E, x, y = None, None, None - Ep, xp, yp = None, None, None + # Load in the WFEs and set them to the corresponding ADC surfaces + print('\nSetting the wavefront errors...') + # null_wfe = np.zeros((3001,3001)) + mcd.ELT_FP.add_wavefront_error(surface_map_fp, xfp, yfp, wfe=wfe_fp) - # Reshaping the arrays - shape = (int(np.sqrt(len(x_expected_list))), int(np.sqrt(len(x_expected_list)))) - x_expected = np.asarray(x_expected_list).reshape(shape) - y_expected = np.asarray(y_expected_list).reshape(shape)#.T + # Propagate through the MICADO system + print('Propagating with Wavefront Errors applied...') + Ep, xp, yp = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=hx, hy=hy) + - # Saving the arrays - save_variable('../data/GD_IFP_obscured_PSF_xcentroid_expected_zern.pkl', x_expected) - save_variable('../data/GD_IFP_obscured_PSF_ycentroid_expected_zern.pkl', y_expected) + # Do the Obscured Airy Pattern Fit for the perturbed PSF + xf, yf, _, _ = mcd.centroid_shift(method='Obscured Airy', initialGuess=(amplitude_guess, xe, ye, 0.0, 0.03)) + + # Resulting centroid shift due to the WFE + x_centroid_list.append(xf) + y_centroid_list.append(yf) + R = np.sqrt( (xf-xe)**2 + (yf-ye)**2 ) * mcd.platescale * 1e6 + print('xf: {}, yf: {}, shift (uas): {}'.format(xf,yf, R)) + + # Clear the memory by emptying the large arrays + mcd = None + wfe = None + E, x, y = None, None, None + Ep, xp, yp = None, None, None - # Reshaping the arrays - shape = (int(np.sqrt(len(x_centroid_list))), int(np.sqrt(len(x_centroid_list)))) - x_centroid = np.asarray(x_centroid_list).reshape(shape) - y_centroid = np.asarray(y_centroid_list).reshape(shape) - - # Saving the arrays - save_variable('../data/GD_IFP_obscured_PSF_xcentroid_found_with_wfe_airy_zern.pkl', x_centroid) - save_variable('../data/GD_IFP_obscured_PSF_ycentroid_found_with_wfe_airy_zern.pkl', y_centroid) + # Reshaping the arrays + shape = (int(np.sqrt(len(x_expected_list))), int(np.sqrt(len(x_expected_list)))) + x_expected = np.asarray(x_expected_list).reshape(shape) + y_expected = np.asarray(y_expected_list).reshape(shape)#.T - plot_distortions() + # Saving the arrays + zoom_value = 97 + p + save_variable('../data/GD_IFP_zoom{:.0f}_obscured_PSF_xcentroid_expected_zern_psd.pkl'.format(zoom_value), x_expected) + save_variable('../data/GD_IFP_zoom{:.0f}_obscured_PSF_ycentroid_expected_zern_psd.pkl'.format(zoom_value), y_expected) + + # Reshaping the arrays + shape = (int(np.sqrt(len(x_centroid_list))), int(np.sqrt(len(x_centroid_list)))) + x_centroid = np.asarray(x_centroid_list).reshape(shape) + y_centroid = np.asarray(y_centroid_list).reshape(shape) + + # Saving the arrays + save_variable('../data/GD_IFP_zoom{:.0f}_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl'.format(zoom_value), x_centroid) + save_variable('../data/GD_IFP_zoom{:.0f}_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl'.format(zoom_value), y_centroid) + + plot_distortions(x_expected, y_expected, x_centroid, y_centroid) diff --git a/scripts/SPIE/Thesis/geometric_distortions_fp_wfe_continuous.py b/scripts/SPIE/Thesis/geometric_distortions_fp_wfe_continuous.py new file mode 100644 index 0000000..0ff6178 --- /dev/null +++ b/scripts/SPIE/Thesis/geometric_distortions_fp_wfe_continuous.py @@ -0,0 +1,214 @@ +import numpy as np +import sys +# sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +# sys.path.insert(0, '/home/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Documents/Published/FourierOptics GitHub release/') +from fourierOptics import * +from pdb import set_trace + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to generate the geometric distortions due to WFEs + off the ADC prisms. The distortions are produced by measuring + the offset of the PSF centroid with and without WFE as a + function of field position. + Script made for the SPIE manuscript. +''' + +def plot_distortions(): + # Coordinates without WFEs + xe = load_variable('../data/GD_IFP_obscured_PSF_xcentroid_expected_fs.pkl') + ye = load_variable('../data/GD_IFP_obscured_PSF_ycentroid_expected_fs.pkl') + + # Coordinates with WFEs + xf = load_variable('../data/GD_IFP_obscured_PSF_xcentroid_found_with_wfe_airy_fs.pkl') + yf = load_variable('../data/GD_IFP_obscured_PSF_ycentroid_found_with_wfe_airy_fs.pkl') + + mcd = ParaxialMICADO(n=3001, m=3001, dpix=50, l=1.49e-3) + + n = xe.shape[0] + X = xe.flatten() + Y = ye.flatten() + dX = xf - xe + dY = yf - ye + dx = dX.flatten() + dy = dY.flatten() + + # Magnitude of the distortions + R = np.sqrt(dx**2 + dy**2) * mcd.platescale * 1e6 + + # Dummy figure, used for the colorbar colors in the quiver plot + im = plt.imshow(np.reshape(R, (n,n))) + plt.close() + + fig = plt.figure() + ax = fig.add_subplot() + for i in range(len(xe)): + ax.plot(xe[:,i], ye[:,i], 'k-', alpha=0.2, linewidth=.5) + ax.plot(xe[i,:], ye[i,:], 'k-', alpha=0.2, linewidth=.5) + ax.quiver(X,Y,dx,dy, angles='xy', edgecolor='k', facecolor='None', linewidth=.5, units='dots') + ax.quiver(X,Y,dx,dy,R, alpha=.5, units='dots') + ax.set_aspect('equal', adjustable='box') + ax.set_xlabel('x (mm)') + ax.set_ylabel('y (mm)') + fig.colorbar(im, ax=ax, label='PSF centroid shift ($\\mu$as)') + plt.savefig('../figures/UPDATED_GD_focal_plane_fs.pdf', bbox_inches='tight', transparent=True) + plt.savefig('../figures/UPDATED_GD_focal_plane_fs.png', bbox_inches='tight', transparent=True) + plt.tight_layout() + plt.show() + + +if __name__ == '__main__': + # plt.rc('text', usetex=True) + plt.rc('font', family='serif') + plt.rc('axes', labelsize='large') + + + # Setup of simulation + l = 1.49e-3 + mcd = ParaxialMICADO(n=3001, m=3001, dpix=50) + # mcd.plot_ray_trace() + + # Get x3 and y3 (the coordinates to apply the ADC WFE to): + # The focal plane area is about 200x200 mm, so we apply the wavefront error to an area + # of 200x200 mm. + xfp, yfp = np.meshgrid(np.linspace(-150, 150, mcd.n), np.linspace(-150, 150, mcd.m)) + # _, xfp, yfp = mcd.propagate(l=l,surface=1, zoom=False, show=False, + # print_progress=False) + + # Generate a wavefront error for the intermediate focal plane. + print('\nSetting the wavefront errors...') + wfe_fp = WavefrontError(mcd.n,mcd.m) + wfe_rms = 15000e-6 + #PSD Wavefront error + + # Loading the phase generated using a Fourier Series + fs_array = load_variable('wfe_fp_phase_array_{}_{}.pkl'.format(mcd.n, mcd.m)) + + # Retrieving the coefficients for the Fourier series + np.random.seed(2183) + coeff = np.random.uniform(-1,1,(101,101)) + wfe_fp.coefficients = coeff + fs_array = fourier_series(xfp, yfp, coeff, Px=30, Py=30, xs=10, ys=10) + + # Setting the initial full window phase array to generate the full window surface height map. + wfe_fp.phase = fs_array + surface_map_fp = wfe_fp.generate_height_map(wfe_rms, 2, 1) + # plot_real_valued_results(surface_map_fp*1e6, xfp, yfp, logscale=False, fig_size=(8,6.4), grid=False, show=True, + # colormap='jet', labeltext='Surface height (nm)') + + # save_variable('../data/MCD_IFP_wfe_object_PSD_50nm.pkl', wfe_fp) + # wfe_fp = load_variable('../data/MCD_IFP_wfe_object_PSD_50nm.pkl') + # surface_map_fp = wfe_fp.data + # plt.show() + + + Elist, Eplist = [], [] + xlist, ylist = [], [] + xplist, yplist = [], [] + x_expected_list, y_expected_list = [], [] + x_centroid_list, y_centroid_list = [], [] + + # Set field points + hxList, hyList = np.mgrid[-1:1:11j,-1:1:11j]-(1/300) + # hxList, hyList = np.mgrid[-.1/11:.1/11:11j,-.1/11:.1/11:11j]-(1/300) # Zoomed in to central area + hxList = hxList.flatten() + hyList = hyList.flatten() + # hxList, hyList = [0.0], [0.0] + + # List of performed simulations + # Defined by PSD + ## GD_IFP : -1:1 + ## GD_IFP_zoom : -1/11 :1/11 + ## GD_IFP_zoom2 : -.5/11:.5/11 + ## GD_IFP_zoom3 : -.1/11:.1/11 + + # Defined by zernikes + ## GD_IFP_obscured_PSF_xcentroid_expected_zern : 50 nm rms 19 zernike modes : -1 : 1 + ## To do; check PSF shape at center field with zern wfe + ## To do; check interpolated WFE at any zoomed in field point. + + # Propagate through the field positions without any WFEs + for i in range(len(hxList)): + hx = hxList[i] + hy = hyList[i] + print('Data point: i={:.0f} of i_max={:.0f}'.format(i, len(hxList))) + print('hx: {} - hy: {}'.format(hx,hy) ) + + # Initializing MICADO paraxial model + mcd = ParaxialMICADO(n=3001, m=3001, dpix=50, entrance_pupil='obscured') + + # Propagate through the MICADO system + E, x, y = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=hx, hy=hy) + + # Find approximation for the centroid and amplitude as a first + # guess for the Obscured Airy pattern fit. + amplitude_guess = np.max(np.abs(E))**2 + xbc, ybc, _, _ = mcd.centroid_shift(method='Barycenter') + print('Barycenter guess: ', xbc, ybc) + + # Do the Obscured Airy Pattern Fit for the unperturbed PSF + xe, ye, _, _ = mcd.centroid_shift(method='Obscured Airy', initialGuess=(amplitude_guess, xbc, ybc, 0.0, 0.03)) + print('Airy result: ', xe, ye) + x_expected_list.append(xe) + y_expected_list.append(ye) + + + + # Load in the WFEs and set them to the corresponding ADC surfaces + print('\nSetting the wavefront errors...') + # null_wfe = np.zeros((3001,3001)) + mcd.ELT_FP.add_wavefront_error(surface_map_fp, xfp, yfp, wfe=wfe_fp) + + + # Propagate through the MICADO system + print('Propagating with Wavefront Errors applied...') + Ep, xp, yp = mcd.propagate(l=l,surface=-1, zoom=False, show=True, + print_progress=False, hx=hx, hy=hy) + + + # Do the Obscured Airy Pattern Fit for the perturbed PSF + xf, yf, _, _ = mcd.centroid_shift(method='Obscured Airy', initialGuess=(amplitude_guess, xe, ye, 0.0, 0.03)) + + # Resulting centroid shift due to the WFE + x_centroid_list.append(xf) + y_centroid_list.append(yf) + R = np.sqrt( (xf-xe)**2 + (yf-ye)**2 ) * mcd.platescale * 1e6 + print('xf: {}, yf: {}, shift (uas): {}'.format(xf,yf, R)) + + data_tuple = (i, hx, hy, xe, ye, xf, yf) + data_name = 'gd_point_{:03d}_hx{}_hy{}.pkl'.format(i, hx, hy) + save_variable('../data/GD_IFP_intermediate_results/{}'.format(data_name), data_tuple) + + # Clear the memory by emptying the large arrays + mcd = None + wfe = None + E, x, y = None, None, None + Ep, xp, yp = None, None, None + + + # Reshaping the arrays + shape = (int(np.sqrt(len(x_expected_list))), int(np.sqrt(len(x_expected_list)))) + x_expected = np.asarray(x_expected_list).reshape(shape) + y_expected = np.asarray(y_expected_list).reshape(shape)#.T + + # Saving the arrays + save_variable('../data/GD_IFP_obscured_PSF_xcentroid_expected_fs.pkl', x_expected) + save_variable('../data/GD_IFP_obscured_PSF_ycentroid_expected_fs.pkl', y_expected) + + + # Reshaping the arrays + shape = (int(np.sqrt(len(x_centroid_list))), int(np.sqrt(len(x_centroid_list)))) + x_centroid = np.asarray(x_centroid_list).reshape(shape) + y_centroid = np.asarray(y_centroid_list).reshape(shape) + + # Saving the arrays + save_variable('../data/GD_IFP_obscured_PSF_xcentroid_found_with_wfe_airy_fs.pkl', x_centroid) + save_variable('../data/GD_IFP_obscured_PSF_ycentroid_found_with_wfe_airy_fs.pkl', y_centroid) + + plot_distortions() + + + diff --git a/scripts/SPIE/Thesis/plotting_geometric_distortions_IFP.py b/scripts/SPIE/Thesis/plotting_geometric_distortions_IFP.py index 40ee220..21142f3 100644 --- a/scripts/SPIE/Thesis/plotting_geometric_distortions_IFP.py +++ b/scripts/SPIE/Thesis/plotting_geometric_distortions_IFP.py @@ -91,9 +91,10 @@ def convert_detector_size(x, y, old=True): def plot_distortions(xe, ye, xp, yp, lims=False, add_detectors=True, save=False, figtitle='geom_distortions', show=True, cmap='magma'): - vmin, vmax = 0.1, 100 + vmin, vmax = 0.01, 100 mcd = ParaxialMICADO(n=3001, m=3001, dpix=50, l=1.49e-3) n = xe.shape[0] + no_of_grids = int(n/121) X = xe.flatten() Y = ye.flatten() dX = xp - xe @@ -103,7 +104,7 @@ def plot_distortions(xe, ye, xp, yp, lims=False, add_detectors=True, # Magnitude of the distortions R = np.sqrt(dx**2 + dy**2) * mcd.platescale * 1e6 - print('RMS of the PSF shifts: {:.4f} μas.'.format(np.sqrt(np.mean((R)**2)))) + print('RMS of the PSF shifts: {:.4f} μas.'.format(np.sqrt(np.mean((R[:121])**2)))) # Dummy figure, used for the colorbar colors in the quiver plot # im = plt.imshow(np.reshape(R, (n,n))) @@ -113,104 +114,160 @@ def plot_distortions(xe, ye, xp, yp, lims=False, add_detectors=True, plt.close() # Initiate the figure - # fig = plt.figure(figsize=(9.6,4.8)) - # ax = fig.add_subplot(121, aspect='equal') - # ax2 = fig.add_subplot(122, aspect='equal') - fig, axes = plt.subplots(ncols=2, - nrows=1, - gridspec_kw={'width_ratios': [1,1]}, - figsize=(9.6,4.8)) - ax, ax2 = axes - ax.set_aspect('equal') - ax2.set_aspect('equal') + fig = plt.figure(figsize=(12,4.8)) + gs = GridSpec(ncols=3, nrows=2, width_ratios=[2,1,1], height_ratios=[1,1])#left=0.05, right=0.48, wspace=0.05) + ax = fig.add_subplot(gs[:, :1]) + ax2 = fig.add_subplot(gs[0, 1]) + ax3 = fig.add_subplot(gs[0, 2]) + ax4 = fig.add_subplot(gs[1, 1]) + ax5 = fig.add_subplot(gs[1, 2]) + + axilist = [ax, ax2, ax3, ax4, ax5] + for axi in axilist: + axi.set_aspect('equal') + # Add the detector layout if so desired. if add_detectors: - add_detector_layout((fig,ax), lw=.8, alpha=1, show_label=False) + add_detector_layout((fig,ax), lw=1, alpha=1, show_label=False) # Plot the vector plot xe = np.asarray(xe[:121]).reshape((11,11)) ye = np.asarray(ye[:121]).reshape((11,11)) - # for i in range(len(xe)): - # ax.plot(xe[:,i], ye[:,i], 'k-', alpha=0.2, linewidth=.5) - # ax.plot(xe[i,:], ye[i,:], 'k-', alpha=0.2, linewidth=.5) - # ax.quiver(X,Y,dx,dy,R, alpha=1, units='dots', zorder=300) - - for k in range(4): - print() + + # Now we iterate over the different iterations of the GD simulations + print('Total number of grids: {}'.format(no_of_grids)) + for k in range(no_of_grids): + + # Selecting the different iterations (11x11 grid, but flattened) ind_end = (k+1) * 121 ind_start = k * 121 dpix = np.diff(Y)[ind_start:ind_end][1] - print(k, ind_start, ind_end, dpix) + # Coordinates of the grid in this iteration xplot = np.asarray(X[ind_start:ind_end]).reshape((11,11)) yplot = np.asarray(Y[ind_start:ind_end]).reshape((11,11)) yplot, xplot = np.mgrid[yplot.min()+dpix/2:yplot.max()-dpix/2:12j, xplot.min()+dpix/2:xplot.max()-dpix/2:12j] - + # Magnitude of the distortion at coordinates xplot, yplot + # Used to color the background Rtemp = np.asarray(R[ind_start:ind_end]).reshape((11,11)) - ax.pcolor(-1*yplot.T, -1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) - # ax2.pcolor(-1*xplot.T, -1*yplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=1, vmax=100)) - shiftx = X[ind_start:ind_end][61]#np.abs(dpix)/2 + shiftx = X[ind_start:ind_end][61] shifty = Y[ind_start:ind_end][60] - print(Y[ind_start:ind_end]) - print('shiftx:', shiftx) - print('shifty:', shifty) - # ax2.pcolor(k*shifty-1*yplot.T, k*shiftx-1*xplot.T, Rtemp.T*10, cmap=cmap, norm=LogNorm(vmin=1, vmax=100)) - if k == 1: ax2.pcolor(shifty-1*yplot.T, shiftx-1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) - elif k == 2: ax2.pcolor(2*shifty-1*yplot.T, 2*shiftx-1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) - elif k == 3: ax2.pcolor(2*shifty-1*yplot.T, 2*shiftx-1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) - if k == 0: + + if k in [0,1,2,3,7]: # Plot the center part pcolor + ax.pcolor(1*shifty-1*yplot.T, 1*shiftx-1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) + # if k > 3: # Plot the zoom and megazoom + ax2.pcolor(2*shifty-1*yplot.T, 2*shiftx-1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) + ax3.pcolor(2*shifty-1*yplot.T, 2*shiftx-1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) + if k in [0,4,5,6,8]: # Plot the left bottom part pcolor + if k!=0: + ax.pcolor(yplot.T, xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) + ax4.pcolor(2*shifty-1*yplot.T, 2*shiftx-1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) + ax5.pcolor(2*shifty-1*yplot.T, 2*shiftx-1*xplot.T, Rtemp.T, cmap=cmap, norm=LogNorm(vmin=vmin, vmax=vmax)) + + # Quiver plots + # Full field plot + if k == 0: ax.quiver(X[ind_start:ind_end],Y[ind_start:ind_end], dx[ind_start:ind_end],dy[ind_start:ind_end], angles='xy', edgecolor='k', facecolor='k', - linewidth=.8, zorder=300, alpha=1, scale=1.4e-5, + linewidth=.8, zorder=300, alpha=1, scale=1.2e-5, units='dots') - elif k in [2,3]: + # Center plot zoom + if k in [0,1,2]: ax2.quiver(X[ind_start:ind_end],Y[ind_start:ind_end], dx[ind_start:ind_end],dy[ind_start:ind_end], angles='xy', edgecolor='k', facecolor='k', - linewidth=.8, zorder=300, alpha=1, scale=0.2e-5, + linewidth=.8, zorder=300, alpha=1, scale=7e-7, units='dots') - # ax.pcolor(xplot * (1 - dpix/105), yplot * (1 - dpix/105), Rtemp, cmap='jet', norm=LogNorm(vmin=1, vmax=100)) - - ax.set_xlabel('x (mm)', fontsize=16) - ax.set_ylabel('y (mm)', fontsize=16) - ax.tick_params(axis='both', direction='in', labelsize=14, zorder=-1, - top=True, right=True) - ax2.set_xlabel('x (mm)', fontsize=16) - ax2.set_ylabel('y (mm)', fontsize=16) - ax2.tick_params(axis='both', direction='in', labelsize=14, zorder=-1, - top=True, right=True) + + # Center plot megazoom + elif k in [2,3,7]: + ax3.quiver(X[ind_start:ind_end],Y[ind_start:ind_end], + dx[ind_start:ind_end],dy[ind_start:ind_end], + angles='xy', edgecolor='k', facecolor='k', + linewidth=.8, zorder=300, alpha=1, scale=1e-6, + units='dots') + + # Left bottom zoom + if k in [4,5]: + ax4.quiver(X[ind_start:ind_end],Y[ind_start:ind_end], + dx[ind_start:ind_end],dy[ind_start:ind_end], + angles='xy', edgecolor='k', facecolor='k', + linewidth=.8, zorder=300, alpha=1, scale=1.2e-5, + units='dots') + + if k in [5,6,8]: + ax5.quiver(X[ind_start:ind_end],Y[ind_start:ind_end], + dx[ind_start:ind_end],dy[ind_start:ind_end], + angles='xy', edgecolor='k', facecolor='k', + linewidth=.8, zorder=300, alpha=1, scale=7e-6, + units='dots') + + # Adjusting how the axes look in all subfigures + for axi in [ax, ax2, ax4]: + axi.set_ylabel('y (mm)', fontsize=16) + for axi in [ax, ax4, ax5]: + axi.set_xlabel('x (mm)', fontsize=16) + for axi in axilist: + axi.tick_params(axis='both', direction='in', labelsize=14, zorder=-1, + top=True, right=True) - - cb = fig.colorbar(im, ax=axes, fraction=0.025, pad=0.04) - cb.set_label(label='PSF centroid shift (μas)', size=14) - cb.ax.tick_params(labelsize=14) - - lbls = 10**np.asarray(cb.ax.get_yticks()) - print(lbls) - newlbls = ['{:.0f}'.format(x) for x in 10**np.asarray(cb.ax.get_yticks())] - newlbls[1] = '{:.1f}'.format(lbls[1]) - newlbls[0] = '{:.1f}'.format(lbls[0]) - cb.ax.set_yticklabels(newlbls) - print(cb.ax.get_yticks()) + # Setting some limits to the subplots if not lims: ax.set_xlim([-125, 115]) ax.set_ylim([-115, 115]) elif lims: ax.set_xlim(lims[0]) ax.set_ylim(lims[1]) - ax2.set_xlim([-2.01,2.75]) - ax2.set_ylim([-2.01,2.75]) - ax.indicate_inset_zoom(ax2, edgecolor="black") + ax2.set_xlim([-1+0.3513,1+0.3513]) + ax2.set_ylim([-1+0.3513,1+0.3513]) + ax3.set_xlim([-0.01+0.3513, 0.01+0.3513]) + ax3.set_ylim([-0.01+0.3513, 0.01+0.3513]) + # ax4.set_xlim([-2.01-84,2.75-84]) + # ax4.set_ylim([-2.01-84,2.75-84]) + ax4.set_xlim([-83.6-1,-83.6+1]) + ax4.set_ylim([-83.6-1,-83.6+1]) + ax5.set_xlim([-0.01-83.6070, 0.01-83.6070]) + ax5.set_ylim([-0.01-83.6070, 0.01-83.6070]) + + # For the center inset + ax.indicate_inset_zoom(ax2, edgecolor='black') + ax2.indicate_inset_zoom(ax3, edgecolor='black') + + # For the off-center inset + ax.indicate_inset_zoom(ax4, edgecolor='black') + ax4.indicate_inset_zoom(ax5, edgecolor='black') + + # after plotting the data, format the labels + ax5.set_xticklabels(['{:.2f}'.format(x) for x in ax5.get_xticks()]) + ax5.set_yticklabels(['{:.2f}'.format(x) for x in ax5.get_yticks()]) + fig.subplots_adjust(top=1.0, - bottom=0.0, - left=0.1, - right=0.8, - hspace=0.135, - wspace=0.3) + bottom=0.08, + left=0.075, + right=0.88, + hspace=0.145, + wspace=0.45) + + # Setting and modifying the colorbar + fig.subplots_adjust(right=0.88) + cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.85]) + cb = fig.colorbar(im, cax=cbar_ax) + cb.set_label(label='PSF centroid shift (μas)', size=14) + cb.ax.tick_params(labelsize=14) + + lbls = 10**np.asarray(cb.ax.get_yticks()) + + newlbls = ['{:.0f}'.format(x) for x in 10**np.asarray(cb.ax.get_yticks())] + newlbls[3] = '{:.1f}'.format(lbls[3]) + newlbls[2] = '{:.1f}'.format(lbls[2]) + newlbls[1] = '{:.2f}'.format(lbls[1]) + newlbls[0] = '{:.2f}'.format(lbls[0]) + cb.ax.set_yticklabels(newlbls) + # plt.tight_layout() if save: plt.savefig('{}.pdf'.format(figtitle), @@ -228,51 +285,127 @@ def plot_distortions(xe, ye, xp, yp, lims=False, add_detectors=True, if __name__ == '__main__': ps = ParaxialMICADO(n=3001, m=3001, dpix=50, l=1.49e-3).platescale - print(ps) plt.rc('font', family='serif') plt.rc('axes', labelsize='large') # convert_detector_size() - # Plotting the original distortions + # Plotting the original distortions, -1/11 :1/11 data_dir = '../data/GD Focal Plane WFEs/' xe_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_xcentroid_expected_zern_psd.pkl') ye_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_ycentroid_expected_zern_psd.pkl') xp_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') yp_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') - # Zoom 1: - xe_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_xcentroid_expected_zern_psd.pkl') - ye_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_ycentroid_expected_zern_psd.pkl') - xp_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') - yp_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + # # Zoom 1: -1/11 :1/11 + # xe_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_xcentroid_expected_zern_psd.pkl') + # ye_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_ycentroid_expected_zern_psd.pkl') + # xp_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + # yp_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') - # Zoom 2 - xe_og2 = load_variable(data_dir + 'GD_IFP_zoom2_obscured_PSF_xcentroid_expected_zern_psd.pkl') - ye_og2 = load_variable(data_dir + 'GD_IFP_zoom2_obscured_PSF_ycentroid_expected_zern_psd.pkl') - xp_og2 = load_variable(data_dir + 'GD_IFP_zoom2_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') - yp_og2 = load_variable(data_dir + 'GD_IFP_zoom2_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + # # Zoom 2: -.5/11:.5/11 + # xe_og2 = load_variable(data_dir + 'GD_IFP_zoom2_obscured_PSF_xcentroid_expected_zern_psd.pkl') + # ye_og2 = load_variable(data_dir + 'GD_IFP_zoom2_obscured_PSF_ycentroid_expected_zern_psd.pkl') + # xp_og2 = load_variable(data_dir + 'GD_IFP_zoom2_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + # yp_og2 = load_variable(data_dir + 'GD_IFP_zoom2_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') - # Zoom 3 - xe_og3 = load_variable(data_dir + 'GD_IFP_zoom3_obscured_PSF_xcentroid_expected_zern_psd.pkl') - ye_og3 = load_variable(data_dir + 'GD_IFP_zoom3_obscured_PSF_ycentroid_expected_zern_psd.pkl') - xp_og3 = load_variable(data_dir + 'GD_IFP_zoom3_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') - yp_og3 = load_variable(data_dir + 'GD_IFP_zoom3_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + # # Zoom 3: -.1/11:.1/11 + # xe_og3 = load_variable(data_dir + 'GD_IFP_zoom3_obscured_PSF_xcentroid_expected_zern_psd.pkl') + # ye_og3 = load_variable(data_dir + 'GD_IFP_zoom3_obscured_PSF_ycentroid_expected_zern_psd.pkl') + # xp_og3 = load_variable(data_dir + 'GD_IFP_zoom3_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + # yp_og3 = load_variable(data_dir + 'GD_IFP_zoom3_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') - xe = np.append(xe_og, [xe_og1, xe_og2, xe_og3]) - ye = np.append(ye_og, [ye_og1, ye_og2, ye_og3]) - xp = np.append(xp_og, [xp_og1, xp_og2, xp_og3]) - yp = np.append(yp_og, [yp_og1, yp_og2, yp_og3]) + # Zoom 97: -1/11: 1/11 + xe_og1 = load_variable(data_dir + 'GD_IFP_zoom97_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og1 = load_variable(data_dir + 'GD_IFP_zoom97_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og1 = load_variable(data_dir + 'GD_IFP_zoom97_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og1 = load_variable(data_dir + 'GD_IFP_zoom97_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') - # xe = np.append(xe_og2, xe_og3) - # ye = np.append(ye_og2, ye_og3) - # xp = np.append(xp_og2, xp_og3) - # yp = np.append(yp_og2, yp_og3) + # Zoom 98: -1/121 : 1/121 + xe_og2 = load_variable(data_dir + 'GD_IFP_zoom98_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og2 = load_variable(data_dir + 'GD_IFP_zoom98_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og2 = load_variable(data_dir + 'GD_IFP_zoom98_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og2 = load_variable(data_dir + 'GD_IFP_zoom98_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # Zoom 99: -1/11**3 : -1/11**3 + xe_og3 = load_variable(data_dir + 'GD_IFP_zoom99_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og3 = load_variable(data_dir + 'GD_IFP_zoom99_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og3 = load_variable(data_dir + 'GD_IFP_zoom99_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og3 = load_variable(data_dir + 'GD_IFP_zoom99_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') - figtitle = '../figures adapted for thesis/Focal Plane/GD_IFP_zern' - plot_distortions(xe, ye, xp, yp, show=True, save=True, figtitle=figtitle, + # Zoom7: -.01/11:.01/11 + xe_og4 = load_variable(data_dir + 'GD_IFP_zoom92_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og4 = load_variable(data_dir + 'GD_IFP_zoom92_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og4 = load_variable(data_dir + 'GD_IFP_zoom92_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og4 = load_variable(data_dir + 'GD_IFP_zoom92_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # Zoom8: -.05/11:.05/11 - not plotted + # Zoom93: -.004/11:.004/11 + xe_og5 = load_variable(data_dir + 'GD_IFP_zoom93_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og5 = load_variable(data_dir + 'GD_IFP_zoom93_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og5 = load_variable(data_dir + 'GD_IFP_zoom93_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og5 = load_variable(data_dir + 'GD_IFP_zoom93_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # # Zoom91: -.01/11:.01/11 but hx and hy shifted to an area near the corner of the image, + # # where distortion is supposed to be large + # xe_og6 = load_variable(data_dir + 'GD_IFP_zoom91_obscured_PSF_xcentroid_expected_zern_psd.pkl') + # ye_og6 = load_variable(data_dir + 'GD_IFP_zoom91_obscured_PSF_ycentroid_expected_zern_psd.pkl') + # xp_og6 = load_variable(data_dir + 'GD_IFP_zoom91_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + # yp_og6 = load_variable(data_dir + 'GD_IFP_zoom91_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # Zoom94: -1/11:1/11 but hx and hy shifted to an area near the corner of the image, + # where distortion is supposed to be large + xe_og6 = load_variable(data_dir + 'GD_IFP_zoom94_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og6 = load_variable(data_dir + 'GD_IFP_zoom94_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og6 = load_variable(data_dir + 'GD_IFP_zoom94_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og6 = load_variable(data_dir + 'GD_IFP_zoom94_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # Zoom95: -1/121:1/121 but hx and hy shifted to an area near the corner of the image, + # where distortion is supposed to be large + xe_og7 = load_variable(data_dir + 'GD_IFP_zoom95_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og7 = load_variable(data_dir + 'GD_IFP_zoom95_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og7 = load_variable(data_dir + 'GD_IFP_zoom95_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og7 = load_variable(data_dir + 'GD_IFP_zoom95_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # Zoom96: -1/11**3:1/11**3 but hx and hy shifted to an area near the corner of the image, + # where distortion is supposed to be large + xe_og8 = load_variable(data_dir + 'GD_IFP_zoom96_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og8 = load_variable(data_dir + 'GD_IFP_zoom96_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og8 = load_variable(data_dir + 'GD_IFP_zoom96_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og8 = load_variable(data_dir + 'GD_IFP_zoom96_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # Zoom 100: - 1/11**4:1/11**4, phase interpolation, center of the field + xe_og9 = load_variable(data_dir + 'GD_IFP_zoom100_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og9 = load_variable(data_dir + 'GD_IFP_zoom100_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og9 = load_variable(data_dir + 'GD_IFP_zoom100_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og9 = load_variable(data_dir + 'GD_IFP_zoom100_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # Zoom 101: - 1/11**4:1/11**4, phase interpolation, center of the field, with shift to left bottom corner + xe_og10 = load_variable(data_dir + 'GD_IFP_zoom101_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye_og10 = load_variable(data_dir + 'GD_IFP_zoom101_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp_og10 = load_variable(data_dir + 'GD_IFP_zoom101_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp_og10 = load_variable(data_dir + 'GD_IFP_zoom101_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + + + + # Combining the data points + # xe = np.append(xe_og, [xe_og1, xe_og2, xe_og3, xe_og4, xe_og5, xe_og6, xe_og7, xe_og8]) + # ye = np.append(ye_og, [ye_og1, ye_og2, ye_og3, ye_og4, ye_og5, ye_og6, ye_og7, ye_og8]) + # xp = np.append(xp_og, [xp_og1, xp_og2, xp_og3, xp_og4, xp_og5, xp_og6, xp_og7, xp_og8]) + # yp = np.append(yp_og, [yp_og1, yp_og2, yp_og3, yp_og4, yp_og5, yp_og6, yp_og7, yp_og8]) + xe = np.append(xe_og, [xe_og1, xe_og2, xe_og3, xe_og6, xe_og7, xe_og8, xe_og9, xe_og10]) + ye = np.append(ye_og, [ye_og1, ye_og2, ye_og3, ye_og6, ye_og7, ye_og8, ye_og9, ye_og10]) + xp = np.append(xp_og, [xp_og1, xp_og2, xp_og3, xp_og6, xp_og7, xp_og8, xp_og9, xp_og10]) + yp = np.append(yp_og, [yp_og1, yp_og2, yp_og3, yp_og6, yp_og7, yp_og8, yp_og9, yp_og10]) + + + + figtitle = '../figures adapted for thesis/Focal Plane/GD_IFP_psd' + plot_distortions(xe, ye, xp, yp, show=False, save=True, figtitle=figtitle, lims=[[-120,120],[-120,120]]) + plt.close() diff --git a/scripts/SPIE/Thesis/plotting_geometric_distortions_IFP_zern.py b/scripts/SPIE/Thesis/plotting_geometric_distortions_IFP_zern.py index a068cbd..f3c83f4 100644 --- a/scripts/SPIE/Thesis/plotting_geometric_distortions_IFP_zern.py +++ b/scripts/SPIE/Thesis/plotting_geometric_distortions_IFP_zern.py @@ -199,7 +199,7 @@ def plot_distortions2(xe, ye, xp, yp, lims=False, add_detectors=True, # ax.plot(xe[i,:], ye[i,:], 'k-', alpha=0.2, linewidth=.5) # ax.quiver(X,Y,dx,dy,R, alpha=1, units='dots', zorder=300) - for k in range(1): + for k in range(3): print() ind_end = (k+1) * 121 ind_start = k * 121 @@ -300,30 +300,30 @@ if __name__ == '__main__': # convert_detector_size() # Plotting the original distortions data_dir = '../data/GD Focal Plane WFEs/' - xe_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_xcentroid_expected_zern.pkl') - ye_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_ycentroid_expected_zern.pkl') - xp_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_xcentroid_found_with_wfe_airy_zern.pkl') - yp_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_ycentroid_found_with_wfe_airy_zern.pkl') + xe_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_xcentroid_expected_psd.pkl') + ye_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_ycentroid_expected_psd.pkl') + xp_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_xcentroid_found_with_wfe_airy_psd.pkl') + yp_og = load_variable(data_dir + 'GD_IFP_obscured_PSF_ycentroid_found_with_wfe_airy_psd.pkl') # # Zoom 1: - # xe_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_xcentroid_expected_zern.pkl') - # ye_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_ycentroid_expected_zern.pkl') - # xp_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_xcentroid_found_with_wfe_airy_zern.pkl') - # yp_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_ycentroid_found_with_wfe_airy_zern.pkl') + xe_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_xcentroid_expected_psd.pkl') + ye_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_ycentroid_expected_psd.pkl') + xp_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_xcentroid_found_with_wfe_airy_psd.pkl') + yp_og1 = load_variable(data_dir + 'GD_IFP_zoom_obscured_PSF_ycentroid_found_with_wfe_airy_psd.pkl') - xe = xe_og#np.append(xe_og, xe_og1) - ye = ye_og#np.append(ye_og, ye_og1) - xp = xp_og#np.append(xp_og, xp_og1) - yp = yp_og#np.append(yp_og, yp_og1) + xe = np.append(xe_og, xe_og1) + ye = np.append(ye_og, ye_og1) + xp = np.append(xp_og, xp_og1) + yp = np.append(yp_og, yp_og1) - # xe = np.append(xe_og2, xe_og3) - # ye = np.append(ye_og2, ye_og3) - # xp = np.append(xp_og2, xp_og3) - # yp = np.append(yp_og2, yp_og3) + xe = np.append(xe_og2, xe_og3) + ye = np.append(ye_og2, ye_og3) + xp = np.append(xp_og2, xp_og3) + yp = np.append(yp_og2, yp_og3) - figtitle = '../figures adapted for thesis/Focal Plane/GD_IFP_zern' - plot_distortions(xe, ye, xp, yp, show=True, save=False, figtitle=figtitle, + figtitle = '../figures adapted for thesis/Focal Plane/GD_IFP_zern2' + plot_distortions2(xe, ye, xp, yp, show=True, save=False, figtitle=figtitle, lims=[[-120,120],[-120,120]]) diff --git a/scripts/SPIE/Thesis/retrieving_average_tilt.py b/scripts/SPIE/Thesis/retrieving_average_tilt.py new file mode 100644 index 0000000..76446ac --- /dev/null +++ b/scripts/SPIE/Thesis/retrieving_average_tilt.py @@ -0,0 +1,86 @@ +import numpy as np +import matplotlib.pyplot as plt +from numpy.polynomial import polynomial +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Documents/Published/FourierOptics GitHub release/') +from fourierOptics import * +from pdb import set_trace + +def polyfit2d(x, y, z, deg): + x = np.asarray(x) + y = np.asarray(y) + z = np.asarray(z) + deg = np.asarray(deg) + vander = polynomial.polyvander2d(x, y, deg) + vander = vander.reshape((-1,vander.shape[-1])) # Wat gebeurt hier? reshape(-1) lijkt niets te doen? + z = z.reshape((vander.shape[0],)) + c = np.linalg.lstsq(vander, z, rcond=None)[0] + return c.reshape(deg+1) + + + +if __name__ == '__main__': + plt.rc('text', usetex=True) + plt.rc('font', family='serif') + + # Initializing MICADO paraxial model + mcd = ParaxialMICADO(n=3001, m=3001, dpix=50) + wfe = WavefrontError(mcd.n, mcd.m) + + # Setup + l = 1.49e-3 + + # Set field points + hxList = [0] + hyList = [0] + + wfe_p1 = load_variable('../data/compressed items/ADC_tWFE_p1_zern_20nm_psd_15nm.pkl') + wfe_p2 = load_variable('../data/compressed items/ADC_tWFE_p2_zern_20nm_psd_15nm.pkl') + wfe_p3 = load_variable('../data/compressed items/ADC_tWFE_p3_zern_20nm_psd_15nm.pkl') + wfe_p4 = load_variable('../data/compressed items/ADC_tWFE_p4_zern_20nm_psd_15nm.pkl') + + # wfe_p1 = load_variable('../data/compressed items/ADC_tWFE_p1_psd_25nm.pkl') + # wfe_p2 = load_variable('../data/compressed items/ADC_tWFE_p2_psd_25nm.pkl') + # wfe_p3 = load_variable('../data/compressed items/ADC_tWFE_p3_psd_25nm.pkl') + # wfe_p4 = load_variable('../data/compressed items/ADC_tWFE_p4_psd_25nm.pkl') + + + + prop_surf = 4 + E, x, y = mcd.propagate(l=l,surface=prop_surf, zoom=False, show=False, + print_progress=True, hx=0, hy=0) + # Save the resulting 'perfect' PSF for later + tx_combined = 0 + ty_combined = 0 + + print(mcd.ELT_ADC.p1.CA) + + for wfe_surf in [wfe_p1, wfe_p2, wfe_p3, wfe_p4]: + + # Don't show anything outside the clear aperture. + mask = np.sqrt(x**2 + y**2) > 40#mcd.ELT_ADC.p1.CA/2 + # Plot P1 + phasescreenMasked = wfe_surf.copy() + wfe_p1[mask] = np.nan + phasescreenMasked[mask] = 0 + + # Bepaling van de gemiddelde tilt voor deze coordinaten, i.e. een lineaire fit. + pf1 = polyfit2d(x.real[~mask], y.real[~mask], phasescreenMasked[~mask], [1,1]) + + # Coefficienten van de fit vertaalt naar graden tilt + o1x = np.arctan(pf1[1,0]) + o1y = np.arctan(pf1[0,1]) + print(o1x, o1y) + o1x_deg = np.degrees(o1x) + + o1y_deg = np.degrees(o1y) + tx_combined += o1x + ty_combined += o1y + print('tx: {} rad ({} deg = {:.2f} arcsec)'.format(o1x, o1x_deg, o1x_deg*3600)) + print('ty: {} rad ({} deg = {:.2f} arcsec)'.format(o1y, o1y_deg, o1y_deg*3600)) + + print('Combined results:\n-----------------') + print('tx: {} rad --> distortion of {} μas'.format(tx_combined, 1601.688*tx_combined*0.267*1e6)) + print('ty: {} rad --> distortion of {} μas'.format(ty_combined, 1601.688*ty_combined*0.267*1e6)) + tr_combined = np.sqrt(tx_combined**2 + ty_combined**2) + print('tr: {} rad --> distortion of {} μas'.format(tr_combined, 1601.688*tr_combined*0.267*1e6)) \ No newline at end of file diff --git a/scripts/SPIE/Thesis/wfe_fp_phase_array_101_101.pkl b/scripts/SPIE/Thesis/wfe_fp_phase_array_101_101.pkl new file mode 100644 index 0000000..6006002 Binary files /dev/null and b/scripts/SPIE/Thesis/wfe_fp_phase_array_101_101.pkl differ diff --git a/scripts/SPIE/Thesis/wfe_fp_phase_array_3001_3001.pkl b/scripts/SPIE/Thesis/wfe_fp_phase_array_3001_3001.pkl new file mode 100644 index 0000000..aa17633 Binary files /dev/null and b/scripts/SPIE/Thesis/wfe_fp_phase_array_3001_3001.pkl differ diff --git a/scripts/SPIE/Thesis/wfe_fp_phase_array_coefficients_101_101.pkl b/scripts/SPIE/Thesis/wfe_fp_phase_array_coefficients_101_101.pkl new file mode 100644 index 0000000..5213b03 Binary files /dev/null and b/scripts/SPIE/Thesis/wfe_fp_phase_array_coefficients_101_101.pkl differ diff --git a/scripts/SPIE/Thesis/wfe_fp_phase_array_coefficients_3001_3001.pkl b/scripts/SPIE/Thesis/wfe_fp_phase_array_coefficients_3001_3001.pkl new file mode 100644 index 0000000..1f06466 Binary files /dev/null and b/scripts/SPIE/Thesis/wfe_fp_phase_array_coefficients_3001_3001.pkl differ diff --git a/scripts/SPIE/figures adapted for thesis/Focal Plane/GD_IFP_psd.pdf b/scripts/SPIE/figures adapted for thesis/Focal Plane/GD_IFP_psd.pdf index 08b49ab..839e75c 100644 Binary files a/scripts/SPIE/figures adapted for thesis/Focal Plane/GD_IFP_psd.pdf and b/scripts/SPIE/figures adapted for thesis/Focal Plane/GD_IFP_psd.pdf differ diff --git a/scripts/SPIE/figures adapted for thesis/Focal Plane/GD_IFP_psd.png b/scripts/SPIE/figures adapted for thesis/Focal Plane/GD_IFP_psd.png index 25cafc5..c799037 100644 Binary files a/scripts/SPIE/figures adapted for thesis/Focal Plane/GD_IFP_psd.png and b/scripts/SPIE/figures adapted for thesis/Focal Plane/GD_IFP_psd.png differ