From 04eb56d529a4730184b54c68980569e47b4c6123 Mon Sep 17 00:00:00 2001 From: yosata Date: Mon, 2 Nov 2020 17:16:33 +0100 Subject: [PATCH] Added off-axis propagation, SPIE scripts --- .gitignore | 4 + .../data/__pycache__/__init__.cpython-37.pyc | Bin 189 -> 0 bytes .../__pycache__/glassCatalog.cpython-37.pyc | Bin 590 -> 0 bytes fourierOptics/fourier.py | 30 +- fourierOptics/micado.py | 405 +++++++++++++++++- fourierOptics/misc.py | 248 ++++++++++- fourierOptics/psfFitter.py | 99 +++++ fourierOptics/wavefrontErrors.py | 1 - scripts/SPIE/centroid_shift_analysis.py | 190 ++++++++ scripts/SPIE/generate_ADC_footprint.py | 94 ++++ scripts/SPIE/generate_ADC_input_field.py | 35 ++ scripts/SPIE/generate_ray_trace_plot.py | 21 + scripts/SPIE/generate_strehl_relation.py | 159 +++++++ scripts/SPIE/generate_wfe.py | 88 ++++ scripts/SPIE/geometric_distortions.py | 162 +++++++ scripts/SPIE/plot_centroid_shift.py | 78 ++++ scripts/SPIE/plot_geometric_distortions.py | 65 +++ scripts/SPIE/plot_strehl.py | 79 ++++ scripts/SPIE/plot_wfe_analysis.py | 97 +++++ .../wavefront_error_frequency_analysis.py | 252 +++++++++++ 20 files changed, 2068 insertions(+), 39 deletions(-) delete mode 100644 fourierOptics/data/__pycache__/__init__.cpython-37.pyc delete mode 100644 fourierOptics/data/__pycache__/glassCatalog.cpython-37.pyc create mode 100644 scripts/SPIE/centroid_shift_analysis.py create mode 100644 scripts/SPIE/generate_ADC_footprint.py create mode 100644 scripts/SPIE/generate_ADC_input_field.py create mode 100644 scripts/SPIE/generate_ray_trace_plot.py create mode 100644 scripts/SPIE/generate_strehl_relation.py create mode 100644 scripts/SPIE/generate_wfe.py create mode 100644 scripts/SPIE/geometric_distortions.py create mode 100644 scripts/SPIE/plot_centroid_shift.py create mode 100644 scripts/SPIE/plot_geometric_distortions.py create mode 100644 scripts/SPIE/plot_strehl.py create mode 100644 scripts/SPIE/plot_wfe_analysis.py create mode 100644 scripts/SPIE/wavefront_error_frequency_analysis.py diff --git a/.gitignore b/.gitignore index cefc713..e844d52 100644 --- a/.gitignore +++ b/.gitignore @@ -103,3 +103,7 @@ dmypy.json # Pyre type checker .pyre/ + +# Exclude SPIE/Data en SPIE/Figures folders +scripts/SPIE/data +scripts/SPIE/figures diff --git a/fourierOptics/data/__pycache__/__init__.cpython-37.pyc b/fourierOptics/data/__pycache__/__init__.cpython-37.pyc deleted file mode 100644 index 7308dd841f5dda151ecb8c6e966eb432893bc5c6..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 189 zcmZ?b<>g`k0?7#_aUl9Jh=2h`Aj1KOi&=m~3PUi1CZpdg`kg1alCW1E;57#@Q-FaVNq0OH~RAd$k5!kEI8!kogA!kW&O!j{6G!jZz6 z!j-}ePpqd9a z<@P>&9O;l=RGHNC&+VCg`mP7*+dnpKv`?4y&UD``eak+b=UajB|ATw&(>JdDwOwt; z8~gM<&qC*3J*DiBE~2%ixu)%_efsat73#eCE%xb!?3-S_ICRQBowt40WiPQ6_UV%k z=M`P~_trk$>*A@s2Yy65q;F7LC3IlNC;RjSA>#|Nww3nj(~<<6ZK8pSQ{@*HYyOyJ zpYC@kW^X{@3;Xo^ruK3&UT=ghdc4xskS_pFb)Ua?PacE85b z^>>why31Rx{b%}^9nv#bsek5Qp5>s)TE*(8>*Q@7#S)z9ZxqGkW>m!Ou z_&9o)8~SN7MzJR+mL%rnr>|rvVg$wunD~{YA6lGRRIHzspI=;}?^2Xskd$AcAD&v2 ztsju#qVJfOm{VDtTC5*XS(1^T2jQ1wCMV|Tm!uY#C?qNr6y;~7CYR`^<(C#^rWW}Z r0A-8yQ-D_Kr{^RV7dt~N)hnpH#bJ}1pHiBWYR3*r8yrA_j}Za^2Fci$ diff --git a/fourierOptics/fourier.py b/fourierOptics/fourier.py index 1f41227..01f9b21 100644 --- a/fourierOptics/fourier.py +++ b/fourierOptics/fourier.py @@ -14,7 +14,8 @@ class FourierOptics: def plane_wave_propagation(self, A, xin, yin, dz=None, oversample=None, show=False, - save=False, scaleInputArray=False, outputCoordinates=None): + save=False, scaleInputArray=False, outputCoordinates=None, + tx=0, ty=0): ''' Propagation through the Fresnel approximation, which is valid in many different situations. A separate code for the Fraunhofer propagation is not striclty necessary, as the exponent @@ -31,6 +32,7 @@ class FourierOptics: scaleInputArray - [tuple] Rescale the resolution of the input array using interpolation. outputCoordinates - [list] list the corners of the desired output coordinate plane. Order: [xmin, xmax, ymin, ymax] + tx, ty - [float] direction of propagation (input in degrees) Returns: E_out - [array] (Complex) Output field distribution @@ -72,8 +74,8 @@ class FourierOptics: _FFTW_Obj() # E-field at z=dz is the fourier transform of A(z=0)*e^(-ikz) - _fx_l = np.fft.fftfreq(_n, d=_dpix) - _fy_l = np.fft.fftfreq(_m, d=_dpix) + _fx_l = np.fft.fftfreq(_n, d=_dpix) + np.tan(np.deg2rad(tx))/self.l + _fy_l = np.fft.fftfreq(_m, d=_dpix) + np.tan(np.deg2rad(ty))/self.l _fx, _fy = np.meshgrid(_fx_l, _fy_l) _kz = self.phase_exponent(_fx, _fy, self.l) # origin at start of array _phaseFactor = np.exp(-1j * _kz * dz) # origin at start of array @@ -92,7 +94,6 @@ class FourierOptics: _E_xyz /= np.sqrt(_IntensityNormFactor) # Normalization back to initial energy. # Coordinates of E_xy0 - # _y, _x = np.mgrid[-_n / 2 * _dpix : _n / 2 * _dpix : 1j*_n, -_m / 2 * _dpix : _m / 2 * _dpix : 1j * _m] _y, _x = yin.real, xin.real if outputCoordinates is not None: @@ -109,19 +110,19 @@ class FourierOptics: _E_xyz = _E_out_ip_real + _E_out_ip_imag*1j _y, _x = np.mgrid[outputCoordinates[0]:outputCoordinates[1]:1j*_n, outputCoordinates[2]:outputCoordinates[3]:1j * _m] - _E_xyz_extent = [np.min(_x.real), np.max(_x.real), np.min(_y.real), np.max(_x.real)] + _E_xyz_extent = [np.min(_x.real), np.max(_x.real), np.min(_y.real), np.max(_y.real)] if show or save: fig = plt.figure(figsize=(8,3)) ax = fig.add_subplot(121) - im = ax.imshow(np.abs(_E_xyz), cmap='hot', extent=_E_xyz_extent)#, vmax=1, vmin=0) + im = ax.imshow(np.abs(_E_xyz), cmap='hot', extent=_E_xyz_extent, origin='lower')#, vmax=1, vmin=0) ax.set_xlabel('x [mm]') ax.set_ylabel('y [mm]') fig.colorbar(im, ax=ax, label='Normalized intensity') ax2 = fig.add_subplot(122) - im2 = ax2.imshow(np.angle(_E_xyz), cmap='hot', extent=_E_xyz_extent) + im2 = ax2.imshow(np.angle(_E_xyz), cmap='hot', extent=_E_xyz_extent, origin='lower') ax2.set_xlabel('x [mm]') ax2.set_ylabel('y [mm]') fig.colorbar(im2, ax=ax2, label='Phase') @@ -137,7 +138,8 @@ class FourierOptics: def fresnel_propagation(self, A, xin, yin, dz=None, oversample=None, show=False, save=False, - scaleInputArray=False, outputCoordinates=None, apply_interpolation=True): + scaleInputArray=False, outputCoordinates=None, apply_interpolation=True, + interpol_type='linear', tx=0, ty=0): ''' Propagation through the Fresnel approximation, which is valid in many different situations. A separate code for the Fraunhofer propagation is not striclty necessary, as the exponent @@ -227,8 +229,8 @@ class FourierOptics: # Evaluate frequencies: _fx = np.fft.fftshift(np.fft.fftfreq(np.shape(_xi)[0], d=_dpix)) _fy = np.fft.fftshift(np.fft.fftfreq(np.shape(_xi)[1], d=_dpix)) - _x = _fx * self.l * dz - _y = _fy * self.l * dz + _x = _fx * self.l * dz + dz * np.tan(np.deg2rad(tx)) + _y = _fy * self.l * dz + dz * np.tan(np.deg2rad(ty)) _x2, _y2 = np.meshgrid(_x, _y) @@ -242,8 +244,8 @@ class FourierOptics: # Performing interpolation to get to the required final plane, done for real and imaginary part. if outputCoordinates and apply_interpolation: - _f_ip_real = interpolate.interp2d(_x.real, _y.real, _E_out.real, kind='linear') - _f_ip_imag = interpolate.interp2d(_x.real, _y.real, _E_out.imag, kind='linear') + _f_ip_real = interpolate.interp2d(_x.real, _y.real, _E_out.real, kind=interpol_type) + _f_ip_imag = interpolate.interp2d(_x.real, _y.real, _E_out.imag, kind=interpol_type) _E_out_ip_real = _f_ip_real(_xip, _yip) _E_out_ip_imag = _f_ip_imag(_xip, _yip) _E_out_ip = _E_out_ip_real + _E_out_ip_imag*1j @@ -264,14 +266,14 @@ class FourierOptics: # Plotting the results _fig = plt.figure(figsize=(8,3)) _ax = _fig.add_subplot(121) - _im = _ax.imshow(np.abs(_E_out_ip), cmap='hot', extent=[_xip.min(), _xip.max(), _yip.min(), _yip.max()]) + _im = _ax.imshow(np.abs(_E_out_ip), cmap='hot', extent=[_xip.min(), _xip.max(), _yip.min(), _yip.max()], origin='lower') _ax.set_xlim([_xip.min(), _xip.max()]) _ax.set_ylim([_yip.min(), _yip.max()]) _ax.set_xlabel('x [mm]') _ax.set_ylabel('y [mm]') _fig.colorbar(_im, ax=_ax, label='Amplitude') _ax2 = _fig.add_subplot(122) - _im2 = _ax2.imshow(np.angle(_E_out_ip), cmap='hot', extent=[_xip.min(), _xip.max(), _yip.min(), _yip.max()]) + _im2 = _ax2.imshow(np.angle(_E_out_ip), cmap='hot', extent=[_xip.min(), _xip.max(), _yip.min(), _yip.max()], origin='lower') _ax2.set_xlim([_xip.min(), _xip.max()]) _ax2.set_ylim([_yip.min(), _yip.max()]) _ax2.set_xlabel('x [mm]') diff --git a/fourierOptics/micado.py b/fourierOptics/micado.py index 4c69960..312990d 100644 --- a/fourierOptics/micado.py +++ b/fourierOptics/micado.py @@ -12,12 +12,20 @@ from .misc import * class ParaxialMICADO: - def __init__(self, l=1.35e-3, n=1001, m=1001, dpix=160): + def __init__(self, l=1.35e-3, n=1001, m=1001, dpix=160, entrance_pupil='ELT'): self.l = l self.dpix = dpix self.n, self.m = n, m self.x, self.y = np.meshgrid(np.linspace(-n/2*dpix, n/2*dpix, n), np.linspace(-m/2*dpix, m/2*dpix, m)) - self.ELT_MASK = Profiles().ELT_pupil(n,m,dpix).astype('complex128') + if entrance_pupil == 'ELT': + self.ELT_MASK = Profiles().ELT_pupil(n,m,dpix).astype('complex128') + elif entrance_pupil == 'obscured': + r1 = 38542/2 + r2 = 10700/2 + self.ELT_MASK = Profiles().obscured_top_hat(r1, r2, n, m, dpix).astype('complex128') + elif entrance_pupil == 'open': + r1 = 38542/2 + self.ELT_MASK = Profiles().top_hat(r1, n, m, dpix).astype('complex128') # Built optical system elements self.ELT_MIRROR = ParabolicMirror(dx=0, dy=0, name='ELT-paraxial', CA=38542, f=6.84022e5) @@ -69,48 +77,71 @@ class ParaxialMICADO: self.last_PSF = (PSF, x, y) - def centroid_shift(self, method='Coupling'): + def centroid_shift(self, method='Coupling', reference_PSF=None, **kwargs): ''' Do a centroid shift calculation for the given method. This is a quick way of doing it. It only returns the centroid shift (and uncertainties) and no other parameters. ''' assert (self.last_PSF is not None), 'PSF necessary for fitting.' - - PSF, x, y = self.last_PSF + if 'initialGuess' in kwargs: + guess_value = kwargs.get('initialGuess') + + if reference_PSF is None: + PSF, x, y = self.last_PSF + elif reference_PSF is not None: + PSF, x, y = reference_PSF fitter = PSFFitter(x.real,y.real,np.abs(PSF)**2) if method == 'Coupling': assert (self.clean_PSF is not None), 'Perfect PSF necessary for comparison.' Eclean, xclean, yclean = self.clean_PSF - fitted_PSF, popt, var, coupling = fitter.fit_custom_field(Eclean, PSF, return_coupling=True) + fitted_PSF, popt, var, coupling = fitter.fit_custom_field(Eclean, PSF, return_coupling=True, **kwargs) x0, y0 = np.asarray(popt).T sx0, sy0 = np.asarray(np.sqrt(var)).T return x0, y0, sx0, sy0, coupling elif method == 'Airy': - fitted_PSF, popt, pcov = fitter.fit_airy_PSF(initialGuess=(0.03**2, 0.0, 0.0, 0.0, 0.03)) - _, x0, y0, _, _ = popt_airy + if 'initialGuess' not in kwargs: + guess_value = (0.03**2, 0.0, 0.0, 0.0, 0.03) + fitted_PSF, popt, pcov = fitter.fit_airy_PSF(initialGuess=guess_value) + _, x0, y0, _, _ = popt + pcov_list = [ np.sqrt(pcov[i][i]) for i in range(len(pcov)) ] + _, sx0, sy0, _, _ = np.asarray(pcov_list).T + + elif method == 'Obscured Airy': + if 'initialGuess' not in kwargs: + guess_value = (0.03**2, 0.0, 0.0, 0.0, 0.03) + fitted_PSF, popt, pcov = fitter.fit_fixed_obscuration_airy_PSF(initialGuess=guess_value) + _, x0, y0, _, _ = popt pcov_list = [ np.sqrt(pcov[i][i]) for i in range(len(pcov)) ] _, sx0, sy0, _, _ = np.asarray(pcov_list).T elif method == 'Gaussian': - fitted_PSF, popt, pcov = fitter.fit_gaussian_PSF(initialGuess=(0.04, 0, 0, 0.01, 0.01, 0, 0)) + if 'initialGuess' not in kwargs: + guess_value = (0.04, 0, 0, 0.01, 0.01, 0, 0) + fitted_PSF, popt, pcov = fitter.fit_gaussian_PSF(initialGuess=guess_value) _, x0, y0, _, _, _, _ = popt pcov_list = [ np.sqrt(pcov[i][i]) for i in range(len(pcov)) ] _, sx0, sy0, _, _, _, _ = np.asarray(pcov_list).T elif method == 'Moffat': - fitted_PSF, popt, pcov = fitter.fit_moffat_PSF(initialGuess=(0.04, 0, 0, 0.01, 0.01, 0, 0, 4)) + if 'initialGuess' not in kwargs: + guess_value = (0.04, 0, 0, 0.01, 0.01, 0, 0, 4) + fitted_PSF, popt, pcov = fitter.fit_moffat_PSF(initialGuess=guess_value) _, x0, y0, _, _, _, _, _ = popt pcov_list = [ np.sqrt(pcov[i][i]) for i in range(len(pcov)) ] _, sx0, sy0, _, _, _, _, _ = np.asarray(pcov_list).T + elif method == 'Barycenter': + x0, y0, sx0, sy0 = fitter.fit_barycenter(amplitude=False) + return x0, y0, sx0, sy0 - def propagate(self, l=None, surface=-1, zoom=False, load_up_to_adc=None, coordinates=None, - show=False, image_output=None, save=None, apply_interpolation=True, print_progress=True): + def propagate(self, l=None, surface=-1, zoom=False, hx=0, hy=0, load_up_to_adc=None, + coordinates=None, show=False, image_output=None, save=None, apply_interpolation=True, + print_progress=True, defocus=None): ''' Propagate the input wavefront through the MICADO paraxial model. Optionally the propagation starts at the ADC, if @@ -140,6 +171,7 @@ class ParaxialMICADO: output_coordinates = self.get_image_coordinates(image_output, coordinatesList, len(surfaces), surface) saveList = self.get_save_image(save, len(surfaces), surface) apply_ip = self.get_image_interpolation(apply_interpolation, zoom, len(surfaces), surface) + tilt_x, tilt_y = self.get_tilt_angles(hx, hy) # Set up the FourierOptics object for propagation fou = FourierOptics(l=l) @@ -196,58 +228,63 @@ class ParaxialMICADO: fzoom = self.ELT_MIRROR.f Eout, x, y = fou.fresnel_propagation(Ein, x, y, dz=self.ELT_MIRROR.f, oversample=None, show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s], - apply_interpolation=apply_ip[s]) + apply_interpolation=apply_ip[s], tx=tilt_x[s], ty=tilt_y[s]) elif surfaces[s] == self.ELT_FP: if print_progress: print('Propagating from ELT-IFP to MCD-COL') Ein = Eout 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]) + apply_interpolation=apply_ip[s], tx=tilt_x[s], ty=tilt_y[s]) elif surfaces[s] == self.ELT_COL: if print_progress: print('Propagating from MCD-COL to MCD-ADC-P1') + x = x.real + y = y.real Ein = Eout * self.ELT_COL.phase_delay(x,y,l, paraxial=True) Eout, x, y = fou.plane_wave_propagation(Ein, x, y, dz=self.MCD.sequentialDistance[3], oversample=None, - show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s]) + show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s], tx=tilt_x[s], ty=tilt_y[s]) elif surfaces[s] == self.ELT_ADC and p == 1: if print_progress: print('Propagating from MCD-ADC-P1 to MCD-ADC-P2') Ein = Eout * self.ELT_ADC.get_phase_delay(x,y,l, theta=None, paraxial=True, prism=1) Eout, x, y = fou.plane_wave_propagation(Ein, x, y, dz=self.MCD.sequentialDistance[4], oversample=None, - show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s]) + show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s], tx=tilt_x[s], ty=tilt_y[s]) p += 1 elif surfaces[s] == self.ELT_ADC and p == 2: if print_progress: print('Propagating from MCD-ADC-P2 to MCD-ADC-P3') Ein = Eout * self.ELT_ADC.get_phase_delay(x,y,l, theta=None, paraxial=True, prism=2) Eout, x, y = fou.plane_wave_propagation(Ein, x, y, dz=self.MCD.sequentialDistance[5], oversample=None, - show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s]) + show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s], tx=tilt_x[s], ty=tilt_y[s]) p += 1 elif surfaces[s] == self.ELT_ADC and p == 3: if print_progress: print('Propagating from MCD-ADC-P3 to MCD-ADC-P4') Ein = Eout * self.ELT_ADC.get_phase_delay(x,y,l, theta=None, paraxial=True, prism=3) Eout, x, y = fou.plane_wave_propagation(Ein, x, y, dz=self.MCD.sequentialDistance[6], oversample=None, - show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s]) + show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s], tx=tilt_x[s], ty=tilt_y[s]) p += 1 elif surfaces[s] == self.ELT_ADC and p == 4: if print_progress: print('Propagating from MCD-ADC-P4 to MCD-CAM') Ein = Eout * self.ELT_ADC.get_phase_delay(x,y,l, theta=None, paraxial=True, prism=4) Eout, x, y = fou.plane_wave_propagation(Ein, x, y, dz=self.MCD.sequentialDistance[7], oversample=None, - show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s]) + show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s], tx=tilt_x[s], ty=tilt_y[s]) p += 1 elif surfaces[s] == self.ELT_DET: if print_progress: print('Propagating from MCD-CAM to MCD-DET') Ein = Eout * self.ELT_CAM.phase_delay(x,y,l, paraxial=True) + if defocus is None: + defocus = 0 if s == surface and zoom: xzoom, yzoom = x.real, y.real Einzoom = Ein - fzoom = self.ELT_CAM.f - Eout, x, y = fou.fresnel_propagation(Ein, x, y, dz=self.ELT_CAM.f, oversample=None, - show=show[s], save=saveList[s], outputCoordinates=output_coordinates[s], apply_interpolation=apply_ip[s]) + fzoom = self.ELT_CAM.f + defocus + Eout, x, y = fou.fresnel_propagation(Ein, x, y, dz=self.ELT_CAM.f+defocus, 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]) s += 1 if zoom: @@ -259,13 +296,15 @@ class ParaxialMICADO: newCoordinates = np.asarray(coordinatesList[ind]) + np.array([xmax, xmax, ymax, ymax]) newCoordinates = newCoordinates.real.tolist() Eout, x, y = fou.fresnel_propagation(Einzoom, xzoom, yzoom, dz=fzoom, oversample=None, - show=show_final_zoom, save=saveList[s-1], outputCoordinates=newCoordinates, apply_interpolation=True) + show=show_final_zoom, save=saveList[s-1], outputCoordinates=newCoordinates, + apply_interpolation=True, interpol_type='linear', tx=tilt_x[s-1], ty=tilt_y[s-1]) if print_progress: print('Propagation finished!') self.last_PSF = (Eout, x, y) return Eout, x, y + def get_surfaces(self): _surfaces = [ self.ELT_MASK, @@ -346,10 +385,330 @@ class ParaxialMICADO: return _apply_ip def set_plate_scale(self): + ''' + Returns the platescale of this paraxial MICADO + model in arcseconds/mm. + ''' _ps = (3600*180/np.pi) / (self.ELT_MIRROR.CA * 53.4466) self.platescale = _ps return _ps + def get_field_angle(self,hx,hy): + fx = hx * 0.0075 + fy = hy * 0.0075 + return fx, fy + + def get_tilt_angles(self, hx, hy): + ''' + Returns the propagation angles at each surface + in degrees. E.g. near the ADC for hx=1, tx=3.61 + degrees. + ''' + _x,_y,_tx,_ty,_z = self.perform_ray_trace_without_adc(hx,hy,0,0) + _tx[0] = _tx[1] + _ty[0] = _ty[1] + tx = np.asarray([_tx[i] for i in [0,1,2,3,5,7,9,11,12,13,13]]).flatten() + ty = np.asarray([_ty[i] for i in [0,1,2,3,5,7,9,11,12,13,13]]).flatten() + return tx, ty + + + def perform_ray_trace(self,hx,hy,px,py): + ''' + Ray trace the paraxial model of ELT-MICADO + + Parameters + ---------- + hx,hy - float or array + Normalized field angles + px,py - float or array + Normalized pupil angles + + Returns + ------- + xout - float or array + x-coordinates on all surfaces + yout - float or array + y-coordinates on all surfaces + zout + z-coordinates of all surfaces + + ''' + xout = [] + yout = [] + + # Show simplified 2D ray trace model + if type(px)==int or type(px)==float: + px = np.array([px]) + if type(py)==int or type(py)==float: + py = np.array([py]) + if type(hx)==int or type(hx)==float: + hx = np.array([hx]) + if type(hy)==int or type(hy)==float: + hy = np.array([hy]) + + xp0 = np.asarray(px) * self.ELT_MIRROR.CA/2 + yp0 = np.asarray(py) * self.ELT_MIRROR.CA/2 + + # for xp in range(len(xp0)): + positions_and_angles = [] + + # Calculate the field angles + txin, tyin = self.get_field_angle(hx, hy) + txin = np.deg2rad(txin) + tyin = np.deg2rad(tyin) + + #Start of system (perhaps at ELT_MASK?) + xp, yp, tx, ty = xp0, yp0, txin, tyin + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # Mask + + # ELT M1 + xp, yp, tx, ty = get_paraxial_lens_rtm(xp, yp, tx, ty, self.ELT_MIRROR.f, self.ELT_MIRROR.f) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # ELT M1 + xp, yp, tx, ty = get_free_space_propagation_rtm(xp, yp, tx, ty, self.ELT_MIRROR.f) + + # ELT-FP + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # ELT FP + xp, yp, tx, ty = get_free_space_propagation_rtm(xp, yp, tx, ty, self.ELT_COL.f) + + # ELT-MCD-COL + xp, yp, tx, ty = get_paraxial_lens_rtm(xp, yp, tx, ty, self.ELT_COL.f, self.ELT_COL.f) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # ELT MCD COL + xp, yp, tx, ty = get_free_space_propagation_rtm(xp, yp, tx, ty, 1209.4718) + + # ELT-MCD-ADC + # First prism: + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx, ty = get_prism_rtm( xp,yp,tx,ty, + thickness=self.ELT_ADC.p1.thickness, + a1x=np.deg2rad(self.ELT_ADC.p1.tilt1x), + a2x=np.deg2rad(self.ELT_ADC.p1.tilt2x), + a1y=np.deg2rad(self.ELT_ADC.p1.tilt1y), + a2y=np.deg2rad(self.ELT_ADC.p1.tilt2y), + n_glass=self.ELT_ADC.p1.find_refractive_index(self.l, unit='millimeter')) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,1) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + + # Second prism: + xp, yp, tx,ty = get_prism_rtm( xp,yp,tx,ty, + thickness=self.ELT_ADC.p2.thickness, + a1x=np.deg2rad(self.ELT_ADC.p2.tilt1x), + a2x=np.deg2rad(self.ELT_ADC.p2.tilt2x), + a1y=np.deg2rad(self.ELT_ADC.p2.tilt1y), + a2y=np.deg2rad(self.ELT_ADC.p2.tilt2y), + n_glass=self.ELT_ADC.p2.find_refractive_index(self.l, unit='millimeter')) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,8) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + + # Third prism + xp, yp, tx,ty = get_prism_rtm( xp,yp,tx,ty, + thickness=self.ELT_ADC.p3.thickness, + a1x=np.deg2rad(self.ELT_ADC.p3.tilt1x), + a2x=np.deg2rad(self.ELT_ADC.p3.tilt2x), + a1y=np.deg2rad(self.ELT_ADC.p3.tilt1y), + a2y=np.deg2rad(self.ELT_ADC.p3.tilt2y), + n_glass=self.ELT_ADC.p3.find_refractive_index(self.l, unit='millimeter')) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,1) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + + # Fourth prism + xp, yp, tx,ty = get_prism_rtm( xp,yp,tx,ty, + thickness=self.ELT_ADC.p4.thickness, + a1x=np.deg2rad(self.ELT_ADC.p4.tilt1x), + a2x=np.deg2rad(self.ELT_ADC.p4.tilt2x), + a1y=np.deg2rad(self.ELT_ADC.p4.tilt1y), + a2y=np.deg2rad(self.ELT_ADC.p4.tilt2y), + n_glass=self.ELT_ADC.p4.find_refractive_index(self.l, unit='millimeter')) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,50) + + # ELT-MCD-CAM: + xp, yp, tx, ty = get_paraxial_lens_rtm(xp, yp, tx, ty, self.ELT_CAM.f, self.ELT_CAM.f) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # After ELT-CAM + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,self.ELT_CAM.f) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # Image + + x_out = np.asarray(positions_and_angles)[:,0] + y_out = np.asarray(positions_and_angles)[:,1] + tx_out = np.asarray(positions_and_angles)[:,2] + ty_out = np.asarray(positions_and_angles)[:,3] + z_out = np.asarray([0, # ELT Mask + 0, # ELT M1 + 6.84022e5, # ELT-FP + 6.84022e5+1418.94, # ELT-MCD-COL + 6.84022e5+1418.94+1209.4718, # ADC + 6.84022e5+1418.94+1209.4718 + 15, # ADC p1 + 6.84022e5+1418.94+1209.4718 + 15 + 1, # Between p1 and p2 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15, # ADC p2 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8, # Between the prisms + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15, # ADC p3 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15 + 1, # Between p3 and p4 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15 + 1 + 15, # ADC p4 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15 + 1 + 15 + 50, # MCD-CAM + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15 + 1 + 15 + 50 + 4278.32701 # Focal plane + ]) + return x_out, y_out, tx_out, ty_out, z_out + + + def perform_ray_trace_without_adc(self,hx,hy,px,py): + ''' + Ray trace the paraxial model of ELT-MICADO + without an ADC in the optical path. + + Parameters + ---------- + hx,hy - float or array + Normalized field angles + px,py - float or array + Normalized pupil angles + + Returns + ------- + xout - float or array + x-coordinates on all surfaces + yout - float or array + y-coordinates on all surfaces + zout + z-coordinates of all surfaces + + ''' + xout = [] + yout = [] + + # Show simplified 2D ray trace model + if type(px)==int or type(px)==float: + px = np.array([px]) + if type(py)==int or type(py)==float: + py = np.array([py]) + if type(hx)==int or type(hx)==float: + hx = np.array([hx]) + if type(hy)==int or type(hy)==float: + hy = np.array([hy]) + + xp0 = np.asarray(px) * self.ELT_MIRROR.CA/2 + yp0 = np.asarray(py) * self.ELT_MIRROR.CA/2 + + # for xp in range(len(xp0)): + positions_and_angles = [] + # l_ind = self.l * 1e3 # wavelength from mm in microns + + # Calculate the field angles + txin, tyin = self.get_field_angle(hx, hy) + txin = np.deg2rad(txin) + tyin = np.deg2rad(tyin) + + #Start of system (perhaps at ELT_MASK?) + xp, yp, tx, ty = xp0, yp0, txin, tyin + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # Mask + + # ELT M1 + xp, yp, tx, ty = get_paraxial_lens_rtm(xp, yp, tx, ty, self.ELT_MIRROR.f, self.ELT_MIRROR.f) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # ELT M1 + xp, yp, tx, ty = get_free_space_propagation_rtm(xp, yp, tx, ty, self.ELT_MIRROR.f) + + # ELT-FP + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # ELT FP + xp, yp, tx, ty = get_free_space_propagation_rtm(xp, yp, tx, ty, self.ELT_COL.f) + + # ELT-MCD-COL + xp, yp, tx, ty = get_paraxial_lens_rtm(xp, yp, tx, ty, self.ELT_COL.f, self.ELT_COL.f) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # ELT MCD COL + xp, yp, tx, ty = get_free_space_propagation_rtm(xp, yp, tx, ty, 1209.4718) + + # ELT-MCD-ADC + # First prism: + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,self.ELT_ADC.p1.thickness) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,1) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + + # Second prism: + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,self.ELT_ADC.p2.thickness) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,8) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + + # Third prism + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,self.ELT_ADC.p3.thickness) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,1) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + + # Fourth prism + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,self.ELT_ADC.p4.thickness) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,50) + + # ELT-MCD-CAM: + xp, yp, tx, ty = get_paraxial_lens_rtm(xp, yp, tx, ty, self.ELT_CAM.f, self.ELT_CAM.f) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # After ELT-CAM + xp, yp, tx,ty = get_free_space_propagation_rtm(xp,yp,tx,ty,self.ELT_CAM.f) + positions_and_angles.append([xp,yp,np.degrees(tx),np.degrees(ty)]) # Image + + x_out = np.asarray(positions_and_angles)[:,0] + y_out = np.asarray(positions_and_angles)[:,1] + tx_out = np.asarray(positions_and_angles)[:,2] + ty_out = np.asarray(positions_and_angles)[:,3] + z_out = np.asarray([0, # ELT Mask + 0, # ELT M1 + 6.84022e5, # ELT-FP + 6.84022e5+1418.94, # ELT-MCD-COL + 6.84022e5+1418.94+1209.4718, # ADC + 6.84022e5+1418.94+1209.4718 + 15, # ADC p1 + 6.84022e5+1418.94+1209.4718 + 15 + 1, # Between p1 and p2 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15, # ADC p2 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8, # Between the prisms + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15, # ADC p3 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15 + 1, # Between p3 and p4 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15 + 1 + 15, # ADC p4 + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15 + 1 + 15 + 50, # MCD-CAM + 6.84022e5+1418.94+1209.4718 + 15 + 1 + 15 + 8 + 15 + 1 + 15 + 50 + 4278.32701 # Focal plane + ]) + return x_out, y_out, tx_out, ty_out, z_out + + + def plot_ray_trace(self, no_adc=False): + ''' + Function that plots a diagram of the ray traced paraxial + model of the MICADO instrument. Note that the prisms of + the ADC do not show their true shape. + ''' + + px = np.linspace(-1,1,21) + py = np.zeros(px.shape) + hx = np.asarray([0]*21) + hy = np.zeros(hx.shape) + if no_adc: + x0, y0, _, _, z0 = self.perform_ray_trace_without_adc(hx,hy,px,py) + hx = np.asarray([1]*21) + x1, y1, _, _, z1 = self.perform_ray_trace_without_adc(hx,hy,px,py) + hx = np.asarray([-1]*21) + x2, y2, _, _, z2 = self.perform_ray_trace_without_adc(hx,hy,px,py) + else: + x0, y0, _, _, z0 = self.perform_ray_trace(hx,hy,px,py) + hx = np.asarray([1]*21) + x1, y1, _, _, z1 = self.perform_ray_trace(hx,hy,px,py) + hx = np.asarray([-1]*21) + x2, y2, _, _, z2 = self.perform_ray_trace(hx,hy,px,py) + + fig = plt.figure(figsize=(15,6)) + ax = fig.add_subplot(111) + for i in range(len(px)): + ax.plot(z0, x0[:,i], 'ro') + ax.plot(z0, x0[:,i], 'r-') + ax.plot(z1, x1[:,i], 'bo') + ax.plot(z1, x1[:,i], 'b-') + ax.plot(z2, x2[:,i], 'go') + ax.plot(z2, x2[:,i], 'g-') + ax.set_xlim([6.84e5,6.915e5]) + ax.set_ylim([-290,290]) + plt.grid() + plt.show() + + class ADC_surfaces: diff --git a/fourierOptics/misc.py b/fourierOptics/misc.py index 162391a..dbb8989 100644 --- a/fourierOptics/misc.py +++ b/fourierOptics/misc.py @@ -1,5 +1,6 @@ import numpy as np import pickle +import matplotlib.pyplot as plt ''' Miscellaneous functions that do not belong specifically to any of the @@ -74,6 +75,58 @@ def perform_extrapolation(x, y, n, m): xnew, ynew = np.meshgrid(xnew, ynew) return xnew, ynew +def plot_results(E, x, y, logscale=False, fig_size=(8,3),grid=False, show=True, colormap='hot'): + x = x.real + y = y.real + _fig = plt.figure(figsize=fig_size) + _ax = _fig.add_subplot(121) + if logscale: + from matplotlib.colors import LogNorm + _im = _ax.imshow(np.abs(E), cmap=colormap, extent=[x.min(), x.max(), y.min(), y.max()], origin='lower', norm=LogNorm()) + else: + _im = _ax.imshow(np.abs(E), cmap=colormap, extent=[x.min(), x.max(), y.min(), y.max()], origin='lower') + if grid: + _ax.grid(True) + _ax.set_xlim([x.min(), x.max()]) + _ax.set_ylim([y.min(), y.max()]) + _ax.set_xlabel('x [mm]') + _ax.set_ylabel('y [mm]') + _fig.colorbar(_im, ax=_ax, label='Amplitude') + _ax2 = _fig.add_subplot(122) + _im2 = _ax2.imshow(np.angle(E), cmap=colormap, extent=[x.min(), x.max(), y.min(), y.max()], origin='lower') + _ax2.set_xlim([x.min(), x.max()]) + _ax2.set_ylim([y.min(), y.max()]) + _ax2.set_xlabel('x [mm]') + _ax2.set_ylabel('y [mm]') + _fig.colorbar(_im2, ax=_ax2, label='Phase') + plt.tight_layout() + if show: + plt.show() + return _fig, _ax, _ax2 + + +def plot_real_valued_results(E, x, y, logscale=False, fig_size=(8,8), grid=False, show=True, colormap='hot', labeltext=None): + x = x.real + y = y.real + _fig = plt.figure(figsize=fig_size) + _ax = _fig.add_subplot(111) + if logscale: + from matplotlib.colors import LogNorm + _im = _ax.imshow(E, cmap=colormap, extent=[x.min(), x.max(), y.min(), y.max()], origin='lower', norm=LogNorm()) + else: + _im = _ax.imshow(E, cmap=colormap, extent=[x.min(), x.max(), y.min(), y.max()], origin='lower') + if grid: + _ax.grid(True) + _ax.set_xlim([x.min(), x.max()]) + _ax.set_ylim([y.min(), y.max()]) + _ax.set_xlabel('x [mm]') + _ax.set_ylabel('y [mm]') + _fig.colorbar(_im, ax=_ax, label=labeltext) + plt.tight_layout() + if show: + plt.show() + return _fig, _ax + def save_variable(filename, variable): ''' @@ -90,4 +143,197 @@ def load_variable(filename): ''' with open(filename, 'rb') as f: variable = pickle.load(f) - return variable \ No newline at end of file + return variable + + +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. + + ''' + dx = (xp - x) * magnification + dy = (yp - y) * magnification + xpn = x + dx + ypn = y + dy + + fig = plt.figure(figsize=(10,10)) + ax = fig.add_subplot(111, aspect=1) + for i in range(len(xp)): + ax.plot(x[:,i], y[:,i], 'k-', alpha=0.2) + ax.plot(x[i,:], y[i,:], 'k-', alpha=0.2) + ax.plot(xpn[:,i], ypn[:,i], 'k-', alpha=1) + ax.plot(xpn[i,:], ypn[i,:], 'k-', alpha=1) + ax.plot(xpn, ypn, 'rx') + ax.set_xlabel('x [mm]') + ax.set_ylabel('y [mm]') + return fig, ax + + +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) + ''' + 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 + + xp1 = x_in / (1 + np.tan(a1x) * np.tan(theta_x_in)) + yp1 = y_in / (1 + np.tan(a1y) * np.tan(theta_y_in)) + + tp1x = -1 * np.arcsin( 1 / n_glass * np.sin(theta_x_in - a1x) ) + a1x + xp2 = xp1 - (thickness + xp1 * (np.tan(a1x) - np.tan(a2x))) * (np.cos(a2x) * np.sin(tp1x)) / np.cos(tp1x - a2x) + + tp1y = -1 * np.arcsin( 1 / n_glass * np.sin(theta_y_in - a1y) ) + a1y + yp2 = yp1 - (thickness + yp1 * (np.tan(a1y) - np.tan(a2y))) * (np.cos(a2y) * np.sin(tp1y)) / np.cos(tp1y - a2y) + + x_out = xp2 * (1 + np.tan(a2x) * np.tan(theta_x_out)) + y_out = yp2 * (1 + np.tan(a2y) * np.tan(theta_y_out)) + + return x_out, y_out, theta_x_out, theta_y_out + + +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) + ''' + x_out = x_in + y_out = y_in + theta_x_out = np.arctan( (x_in + -f_lens_x * np.tan(theta_x_in)) / -f_lens_x ) + theta_y_out = np.arctan( (y_in + -f_lens_y * np.tan(theta_y_in)) / -f_lens_y ) + return x_out, y_out, theta_x_out, theta_y_out + + +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) + ''' + 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 diff --git a/fourierOptics/psfFitter.py b/fourierOptics/psfFitter.py index 3f9aca6..6ca2886 100644 --- a/fourierOptics/psfFitter.py +++ b/fourierOptics/psfFitter.py @@ -16,6 +16,29 @@ class PSFFitter: self.y = y self.psfdata = psfdata + def fit_barycenter(self, xy=None, E=None, amplitude=True): + ''' + Calculates the weighted mean of the PSF. + This is very naive way to find the centroid position. + The function assumes the Efield is the amplitude and not (yet) + the intensity. + ''' + if xy is None and E is None: + x, y = self.x, self.y + E = self.psfdata + else: + x, y = xy + + if amplitude: + I = np.abs(E)**2 + else: + I = E + _x0 = np.sum(x * I) / np.sum(I) + _sx0 = np.sqrt(1 / (x.size-1) * np.sum( (x - _x0)**2 * I) ) + _y0 = np.sum(y * I) / np.sum(I) + _sy0 = np.sqrt(1 / (y.size-1) * np.sum( (y - _y0)**2 * I) ) + return _x0, _y0, _sx0, _sy0 + def gaussian_PSF(self, xy, amplitude, x0, y0, sigma_x, sigma_y, theta, offset): ''' @@ -131,6 +154,82 @@ class PSFFitter: return self.airy_PSF((self.x,self.y), *_popt).reshape(self.x.shape), _popt, _pcov + def obscured_airy_PSF(self, xy, amplitude, x0, y0, epsilon, offset, radius): + ''' + 2D Obscured Airy pattern with offset included + + Input: + xy - Combined x,y coordinates. x and y should be a meshgrid + amplitude - Amplitude of the profile, i.e. maximum of the profile + x0, y0 - Shift in x and y of the central peak + epsilon - Obscuration ratio, e.g. due to the secondary mirror + offset - Noise floor or sky background level + ''' + _x,_y = xy + _x0 = float(x0) + _y0 = float(y0) + _r = np.sqrt((_x-_x0)**2 + (_y-_y0)**2) + _Rz = 1.2196698912665045 + _inp = (np.pi * _r) / (radius/_Rz) + _J1 = j1(_inp) + _eJ1 = j1(epsilon * _inp) + _fraction = np.nan_to_num(2 * _J1 / _inp, nan=1.0) + _efraction = np.nan_to_num(2 * epsilon * _eJ1 / _inp, nan=1.0) + _a = offset + amplitude * (_fraction - _efraction)**2 + return _a.ravel() + + + def fit_obscured_airy_PSF(self, initialGuess=None): + ''' + Repackaging of scipy's curve_fit to fit an Airy pattern to the data. + No arguments are necessary, but an initial guess of the parameters + is recommended. + + Input: + initialGuess - Tuple, containing (amplitude, x0, y0, epsilon, offset, radius) + ''' + _popt, _pcov = opt.curve_fit(self.obscured_airy_PSF, (self.x,self.y), self.psfdata.flatten(), p0=initialGuess) + return self.obscured_airy_PSF((self.x,self.y), *_popt).reshape(self.x.shape), _popt, _pcov + + + def fixed_obscuration_airy_PSF(self, xy, amplitude, x0, y0, offset, radius): + ''' + 2D Obscured Airy pattern with offset included + + Input: + xy - Combined x,y coordinates. x and y should be a meshgrid + amplitude - Amplitude of the profile, i.e. maximum of the profile + x0, y0 - Shift in x and y of the central peak + offset - Noise floor or sky background level + ''' + _x,_y = xy + _epsilon = 0.2776 + _x0 = float(x0) + _y0 = float(y0) + _r = np.sqrt((_x-_x0)**2 + (_y-_y0)**2) + _Rz = 1.2196698912665045 + _inp = (np.pi * _r) / (radius/_Rz) + _J1 = j1(_inp) + _eJ1 = j1(_epsilon * _inp) + _fraction = np.nan_to_num(2 * _J1 / _inp, nan=1.0) + _efraction = np.nan_to_num(2 * _epsilon * _eJ1 / _inp, nan=1.0) + _a = offset + amplitude * (_fraction - _efraction)**2 + return _a.ravel() + + + def fit_fixed_obscuration_airy_PSF(self, initialGuess=None): + ''' + Repackaging of scipy's curve_fit to fit an Airy pattern to the data. + No arguments are necessary, but an initial guess of the parameters + is recommended. + + Input: + initialGuess - Tuple, containing (amplitude, x0, y0, epsilon, offset, radius) + ''' + _popt, _pcov = opt.curve_fit(self.fixed_obscuration_airy_PSF, (self.x,self.y), self.psfdata.flatten(), p0=initialGuess) + return self.fixed_obscuration_airy_PSF((self.x,self.y), *_popt).reshape(self.x.shape), _popt, _pcov + + def reference_PSF(self, p, xy, perturbedPSF): ''' Residual function used for the mimizer functions in diff --git a/fourierOptics/wavefrontErrors.py b/fourierOptics/wavefrontErrors.py index 3f3e25f..2f22c22 100644 --- a/fourierOptics/wavefrontErrors.py +++ b/fourierOptics/wavefrontErrors.py @@ -4,7 +4,6 @@ import types from math import ceil from scipy.special import factorial, j1 - class WavefrontError: ''' This class contains methods to describe wavefront errors; diff --git a/scripts/SPIE/centroid_shift_analysis.py b/scripts/SPIE/centroid_shift_analysis.py new file mode 100644 index 0000000..753f3c1 --- /dev/null +++ b/scripts/SPIE/centroid_shift_analysis.py @@ -0,0 +1,190 @@ +import numpy as np +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to analyse the shift of the PSF centroid as a function + of the RMS value of the WFE. Higher wavefront errors will + result in a larger shift of the centroid. + Script made for the SPIE manuscript. +''' + + +if __name__ == '__main__': + plt.rc('font', family='serif') + plt.rc('text', usetex=True) + plt.rc('axes', labelsize='large') + + # Setting some input parameters + dpix = 160 + n, m = 1001, 1001 + l = 1.49e-3 # wavelength + zern_rms = 25e-5 # RMS of the zernike component of the WFE + midspat_rms = 25e-5 # RMS of the PSD component of the WFE + rmsList = 10**np.linspace(-6, -3, 25) # 1 to 1000 nm WFE RMS in mm, equidistant in log-space + input_field_location = './data/ADC_input_field_1001x1001_1.49_onaxis_obscured.pkl' + + # Initiate ELT-MICADO and wavefront error objects + mcd = ParaxialMICADO(l=l, n=n, m=m, dpix=dpix, entrance_pupil='open') + wfe = WavefrontError(n,m) + + print ('Your resolution is {}x{} with a pixel size of {}'.format(n,m,dpix)) + + # Do the unperturbed propagation for a reference PSF + E, x, y = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + + # Generate a PSD surface + psd_surface = wfe.generate_height_map(midspat_rms, 2, 1) + + # Fit the PSD surface using zernike terms + # Generating a normalized grid for the Zernike evaluation + yz,xz = np.mgrid[-1:1:wfe.n*1j, -1:1:wfe.m*1j] + _r = np.sqrt(xz**2 + yz**2) + _mask = _r > 1 + coeffic = wfe.fit_zernike_coefficients(psd_surface, modes=66) + + # StrehlList_zern, StrehlList_psd = [], [] + rms_zern_list, rms_psd_list = [], [] + + # psd lists + popt_airy_list_psd, pcov_airy_list_psd = [], [] + popt_cust_list_psd, custom_var_list_psd, strehl_list_psd = [], [], [] + Exyzlist_psd = [] + + # zern lists + popt_airy_list_zern, pcov_airy_list_zern = [], [] + popt_cust_list_zern, custom_var_list_zern, strehl_list_zern = [], [], [] + Exyzlist_zern = [] + + # Now iterate over the wavefront rms values. + for rms in rmsList: + # Print results + print('RMS: {:.3f} nm'.format(rms*1e6)) + + + ######################################## + # WFE described by Zernike polynomials # + ######################################## + # Make a new perturbed wavefront error with the required RMS + coeffic = wfe.scale_zernike_coefficients(coeffic, rms) + hzern = wfe.zernike_cartesian(coeffic, xz, yz) + + # Add the new wavefront error + mcd.ELT_ADC.p1.add_wavefront_error(hzern) + + # Performing the propagation using Fresnel approximation to the focal plane. + Exyz_zern,xzern,yzern = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + + psfFit_zern = PSFFitter(xzern.real,yzern.real,np.abs(Exyz_zern)**2) + Exyzlist_zern.append(Exyz_zern) + + # Fit using airy pattern + amplitude_guess = np.max(np.abs(Exyz_zern))**2 + airyFitted, popt_airy, pcov_airy = psfFit_zern.fit_fixed_obscuration_airy_PSF(initialGuess=(amplitude_guess, 0.0, 0.0, 0.0, 0.03)) + popt_airy_list_zern.append(popt_airy) + pcov_airy_list_zern.append([ np.sqrt(pcov_airy[i][i]) for i in range(len(pcov_airy)) ]) + strehl_list_zern.append( np.abs(psfFit_zern.get_coupling(E, Exyz_zern))**2 ) + ######################################## + + + + ######################################## + # WFE described by PSD profile # + ######################################## + factor = rms / midspat_rms + h_psd = wfe.scale_rms(factor, data=psd_surface) + + # Add the new wavefront error + mcd.ELT_ADC.p1.add_wavefront_error(h_psd) + + # Performing the propagation using Fresnel approximation to the focal plane. + Exyz_psd,xpsd,ypsd = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + + psfFit_psd = PSFFitter(xpsd.real,ypsd.real,np.abs(Exyz_psd)**2) + + # Fit using coupling integral + Exyzlist_psd.append(Exyz_psd) + + # Fit using airy pattern + amplitude_guess = np.max(np.abs(Exyz_psd))**2 + airyFitted, popt_airy, pcov_airy = psfFit_psd.fit_fixed_obscuration_airy_PSF(initialGuess=(amplitude_guess, 0.0, 0.0, 0.0, 0.03)) + popt_airy_list_psd.append(popt_airy) + pcov_airy_list_psd.append([ np.sqrt(pcov_airy[i][i]) for i in range(len(pcov_airy)) ]) + strehl_list_psd.append( np.abs(psfFit_psd.get_coupling(E, Exyz_psd))**2 ) + ######################################## + + + # Find the RMS value of the wavefront errros + rms_zern = np.sqrt(np.mean(np.asarray(hzern[~_mask])**2)) + rms_psd = np.sqrt(np.mean(np.asarray(h_psd)**2)) + rms_zern_list.append(rms_zern) + rms_psd_list.append(rms_psd) + + print('Strehl: {:.3f} ... {:.3f}'.format(strehl_list_psd[-1], strehl_list_zern[-1])) + + + + # PSD + # Optimum fit parameters for all scaled surface errors + _, x0_airy_psd, y0_airy_psd, _, _ = np.asarray(popt_airy_list_psd).T + + # Errors in the optimum fit parameters for all scaled surface errors + _, x0_airy_psd_err, y0_airy_psd_err, _, _ = np.asarray(pcov_airy_list_psd).T + + # Distance from original centroid at (0,0) + r_airy_psd = np.sqrt(x0_airy_psd**2+y0_airy_psd**2) + r_err_psd_airy = np.sqrt( (x0_airy_psd**2*x0_airy_psd_err**2+y0_airy_psd**2*y0_airy_psd_err**2) / (x0_airy_psd**2+y0_airy_psd**2)) + + rms_psd_nm = np.asarray(rms_psd_list)*1e6 + r_airy_psd_mas = r_airy_psd * mcd.platescale * 1000 + r_err_psd_airy_mas = r_err_psd_airy * mcd.platescale * 1000 + + + # Zernikes + # Optimum fit parameters for all scaled surface errors + _, x0_airy_zern, y0_airy_zern, _, _ = np.asarray(popt_airy_list_zern).T + + # Errors in the optimum fit parameters for all scaled surface errors + _, x0_airy_zern_err, y0_airy_zern_err, _, _ = np.asarray(pcov_airy_list_zern).T + + # Distance from original centroid at (0,0) + r_airy_zern = np.sqrt(x0_airy_zern**2+y0_airy_zern**2) + r_err_zern_airy = np.sqrt( (x0_airy_zern**2*x0_airy_zern_err**2+y0_airy_zern**2*y0_airy_zern_err**2) / (x0_airy_zern**2+y0_airy_zern**2)) + + rms_zern_nm = np.asarray(rms_zern_list)*1e6 + r_airy_zern_mas = r_airy_zern * mcd.platescale * 1000 + r_err_zern_airy_mas = r_err_zern_airy * mcd.platescale * 1000 + + save_variable('./data/centroid_shift_zernikes.pkl', (rms_zern_nm, r_airy_zern_mas, r_err_zern_airy_mas, strehl_list_zern)) + save_variable('./data/centroid_shift_psd.pkl', (rms_psd_nm, r_airy_psd_mas, r_err_psd_airy_mas, strehl_list_psd)) + + + # Plot some of the point spread functions + fig = plt.figure(figsize=(9,3)) + ax1 = fig.add_subplot(131) + ax1.imshow(np.abs(Exyzlist_zern[0][450:550,450:550]), cmap='hot', interpolation='bicubic') + ax1.set_xticks([]) + ax1.set_yticks([]) + ax2 = fig.add_subplot(132) + ax2.imshow(np.abs(Exyzlist_zern[-5][450:550,450:550]), cmap='hot', interpolation='bicubic') + ax2.set_xticks([]) + ax2.set_yticks([]) + ax4 = fig.add_subplot(133) + ax4.imshow(np.abs(Exyzlist_psd[-5][450:550,450:550]), cmap='hot', interpolation='bicubic') + ax4.set_xticks([]) + ax4.set_yticks([]) + plt.tight_layout() + plt.savefig('./figures/PSF_shape.png', bbox_inches='tight', transparent=True) + plt.savefig('./figures/PSF_shape.pdf', bbox_inches='tight', transparent=True) + plt.show() + \ No newline at end of file diff --git a/scripts/SPIE/generate_ADC_footprint.py b/scripts/SPIE/generate_ADC_footprint.py new file mode 100644 index 0000000..b262fa5 --- /dev/null +++ b/scripts/SPIE/generate_ADC_footprint.py @@ -0,0 +1,94 @@ +import numpy as np +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to generate two plots of the footprint on MCD-CAM lens, + showing the impact of the WFEs on the amplitude and phase + near the exit pupil. + Script made for the SPIE manuscript. +''' + +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 = [-1,-1,0, 1,1] + hyList = [-1, 1,0,-1,1] + + # 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) ) + + # Propagate through the MICADO system + prop_surf = 4 + E, x, y = mcd.propagate(l=l,surface=prop_surf, zoom=False, show=False, + print_progress=False, hx=hx, hy=hy) + # Save the resulting 'perfect' PSF for later + if prop_surf == 4: + # Don't show anything outside the clear aperture. + mask = np.sqrt(x**2 + y**2) > mcd.ELT_ADC.p1.CA/2 + + # Plot P1 + footprint_p1 = E.copy() + footprint_p1[mask]=np.nan + # figp1, axp1, axp2 = plot_results(footprint_p1[850:2150,850:2150], x[850:2150,850:2150], y[850:2150,850:2150], + # fig_size=(12,5), show=False) + figp1, axp1 = plot_real_valued_results(np.abs(footprint_p1)[850:2150,850:2150], x[850:2150,850:2150], y[850:2150,850:2150], + labeltext='Amplitude', show=False, fig_size=None) + plt.savefig('./figures/ELT-MCD-ADC_footprint_no_wfe_{}_{}.pdf'.format(hx,hy), bbox_inches='tight', transparent=True) + plt.savefig('./figures/ELT-MCD-ADC_footprint_no_wfe_{}_{}.png'.format(hx,hy), bbox_inches='tight', transparent=True) + plt.show() + + + + # Load in the WFEs and set them to the corresponding ADC surfaces + wfe_p1 = load_variable('./data/ADC_tWFE_p1_psd_25nm.pkl') + wfe_p2 = load_variable('./data/ADC_tWFE_p2_psd_25nm.pkl') + wfe_p3 = load_variable('./data/ADC_tWFE_p3_psd_25nm.pkl') + wfe_p4 = load_variable('./data/ADC_tWFE_p4_psd_25nm.pkl') + + mcd.ELT_ADC.p1.add_wavefront_error(wfe_p1) + mcd.ELT_ADC.p2.add_wavefront_error(wfe_p2) + mcd.ELT_ADC.p3.add_wavefront_error(wfe_p3) + mcd.ELT_ADC.p4.add_wavefront_error(wfe_p4) + + + + # Iterating over the field positions. + print('\nPropagation with Wavefront Errors applied') + prop_surf = 4 + for i in range(len(hxList)): + hx = hxList[i] + hy = hyList[i] + print('hx: {} - hy: {}'.format(hx,hy) ) + + # Propagate through the MICADO system + Ep, xp, yp = mcd.propagate(l=l,surface=prop_surf+4, zoom=False, show=False, + print_progress=False, hx=hx, hy=hy) + + if prop_surf == 4: + # Don't show anything outside the clear aperture. + mask = np.sqrt(x**2 + y**2) > mcd.ELT_ADC.p1.CA/2 + + # Plot P1 + footprint_p1 = Ep.copy() + figp1, axp1, axp2 = plot_results(footprint_p1[1050:1350,1150:1450], x[1050:1350,1150:1450], y[1050:1350,1150:1450], + fig_size=(12,5), show=False) + plt.savefig('./figures/ELT-MCD-COL_footprint_with_wfe_{}_{}_zoom.pdf'.format(hx,hy), bbox_inches='tight', transparent=True) + plt.savefig('./figures/ELT-MCD-COL_footprint_with_wfe_{}_{}_zoom.png'.format(hx,hy), bbox_inches='tight', transparent=True) + plt.show() + diff --git a/scripts/SPIE/generate_ADC_input_field.py b/scripts/SPIE/generate_ADC_input_field.py new file mode 100644 index 0000000..3642aae --- /dev/null +++ b/scripts/SPIE/generate_ADC_input_field.py @@ -0,0 +1,35 @@ +import numpy as np +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to propagate an input field to a surface in the + ELT-MICADO system. Here we propagate to the first ADC + surface and save the result, saving us a lot of time + in later analyses. + Script made for the SPIE manuscript. +''' + +if __name__ == '__main__': + plt.rc('font', family='serif') + plt.rc('text', usetex=True) + + # Setting some input parameters + dpix = 160 + n, m = 1001, 1001 + l = 1.49e-3 # wavelength + input_field_location = './data/ADC_input_field_{}x{}_{:.2f}_onaxis_obscured.pkl'.format(n,m,l*1e3) + + + # Initiate ELT-MICADO and wavefront error objects + mcd = ParaxialMICADO(l=l, n=n, m=m, dpix=dpix, entrance_pupil='obscured') + + print ('Your resolution is {}x{} with a pixel size of {}'.format(n,m,dpix)) + + # Propagation up to the first ADC surface. + E, x, y = mcd.propagate(l=l,surface=3, zoom=False, show=True, + print_progress=True, hx=0, hy=0) + save_variable(input_field_location, (E,x,y)) \ No newline at end of file diff --git a/scripts/SPIE/generate_ray_trace_plot.py b/scripts/SPIE/generate_ray_trace_plot.py new file mode 100644 index 0000000..60face9 --- /dev/null +++ b/scripts/SPIE/generate_ray_trace_plot.py @@ -0,0 +1,21 @@ +import numpy as np +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to generate the ray traced plot of the paraxial + ELT-MICADO model. + Script made for the SPIE manuscript. +''' + +if __name__ == '__main__': + plt.rc('text', usetex=True) + plt.rc('font', family='serif') + + # Initializing MICADO paraxial model + mcd = ParaxialMICADO(l=1.49e-3, n=3001, m=3001, dpix=50) + mcd.plot_ray_trace() + \ No newline at end of file diff --git a/scripts/SPIE/generate_strehl_relation.py b/scripts/SPIE/generate_strehl_relation.py new file mode 100644 index 0000000..f58f73b --- /dev/null +++ b/scripts/SPIE/generate_strehl_relation.py @@ -0,0 +1,159 @@ +import numpy as np +import os +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to compare the Strehl ratio of 50 different instances + of a generated WFE. Higher spatial frequency perturbations + should follow the Marechal approximation more closely than + lower order shape errors only. We verify this using Power + Spectral Density and Zernike modes (66 coefficients). + Script made for the SPIE manuscript. +''' + + +def expected_strehl2(rms, l): + return np.exp( -(2*np.pi*rms/l)**2 ) + + +def append_to_existing_data(A, filename): + if os.path.exists(filename): + A_o = load_variable(filename) + A2 = np.append(A, A_o) + else: + A2 = A + return A2 + + +if __name__ == '__main__': + # Setting some input parameters + dpix = 160 + n, m = 1001, 1001 + l = 1.49e-3 # wavelength + zern_rms = 25e-5 # RMS of the zernike component of the WFE + midspat_rms = 25e-5 # RMS of the PSD component of the WFE + rmsList = np.linspace(1e-6, 1e-3, 25) # 1 to 1000 nm WFE RMS in mm + number_of_iterations = 50 + input_field_location = './data/ADC_input_field_1001x1001_1.49_onaxis_obscured.pkl' + + + # Initiate ELT-MICADO and wavefront error objects + mcd = ParaxialMICADO(l=l, n=n, m=m, dpix=dpix, entrance_pupil='open') + wfe = WavefrontError(n,m) + + print ('Your resolution is {}x{} with a pixel size of {}'.format(n,m,dpix)) + + # Propagate unperturbed field + E, x, y = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + + # Initiate psf fitter object to calculate the strehl ratios + psf = PSFFitter(x, y, E) + + # Create some empty lists for later + strehl_zern_complete, strehl_psd_complete = [], [] + rms_zern_list_complete, rms_psd_list_complete = [], [] + + + # Repeat the analysis 50 times + for i in range(number_of_iterations): + + print('Iteration {}/{}'.format(i+1,number_of_iterations)) + + # Generate a PSD surface + psd_surface = wfe.generate_height_map(midspat_rms, 2, 1) + + # Fit the PSD surface using zernike terms + # Generating a normalized grid for the Zernike evaluation + yz,xz = np.mgrid[-1:1:wfe.n*1j, -1:1:wfe.m*1j] + _r = np.sqrt(xz**2 + yz**2) + _mask = _r > 1 + coeffic = wfe.fit_zernike_coefficients(psd_surface, modes=66) + + StrehlList_zern, StrehlList_psd = [], [] + rms_zern_list, rms_psd_list = [], [] + + + # Compare results at different RMS values + for rms in rmsList: + + ######################################## + # WFE described by Zernike polynomials # + ######################################## + # Make a new perturbed wavefront error with the required RMS + coeffic = wfe.scale_zernike_coefficients(coeffic, rms) + hzern = wfe.zernike_cartesian(coeffic, xz, yz) + + # Add the new wavefront error + mcd.ELT_ADC.p1.add_wavefront_error(hzern) + + # Performing the propagation using Fresnel approximation to the focal plane. + Exyz_zern,xzern,yzern = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + Strehl_zern = np.abs(psf.get_coupling(E, Exyz_zern))**2 + StrehlList_zern.append(Strehl_zern) + ######################################## + + + ######################################## + # WFE described by PSD profile # + ######################################## + factor = rms / midspat_rms + h_psd = wfe.scale_rms(factor, data=psd_surface) + + # Add the new wavefront error + mcd.ELT_ADC.p1.add_wavefront_error(h_psd) + + # Performing the propagation using Fresnel approximation to the focal plane. + Exyz_psd,xpsd,ypsd = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + Strehl_psd = np.abs(psf.get_coupling(E, Exyz_psd))**2 + StrehlList_psd.append(Strehl_psd) + ######################################## + + # Find the RMS value of the wavefront errros + rms_zern = np.sqrt(np.mean(np.asarray(hzern[~_mask])**2)) + rms_psd = np.sqrt(np.mean(np.asarray(h_psd)**2)) + rms_zern_list.append(rms_zern) + rms_psd_list.append(rms_psd) + + # Print results + print('RMS: {:.3f}, S_zern: {:.3f}, S_psd: {:.3f}, \ + RMS_zern: {:.3f}, RMS_psd: {:.3f}'.format( + rms*1e6, Strehl_zern, Strehl_psd, rms_zern*1e6, rms_psd*1e6)) + + + # Put the results in the correct list + strehl_zern_complete.append(StrehlList_zern) + strehl_psd_complete.append(StrehlList_psd) + rms_zern_list_complete.append(rms_zern_list) + rms_psd_list_complete.append(rms_psd_list) + + rms_zern_list_complete = np.asarray(rms_zern_list_complete) + rms_psd_list_complete = np.asarray(rms_psd_list_complete) + + # Import old results if they are available and append them to the results. + rmsList_present = rmsList.copy() + rmsList = append_to_existing_data(rmsList, './data/SR_RMS_values_obscured.pkl') + strehl_zern_complete = append_to_existing_data(strehl_zern_complete, './data/SR_strehl_zern_obscured.pkl') + strehl_psd_complete = append_to_existing_data(strehl_psd_complete, './data/SR_strehl_psd_obscured.pkl') + rms_zern_list_complete = append_to_existing_data(rms_zern_list_complete, './data/SR_rms_zern_complete_obscured.pkl') + rms_psd_list_complete = append_to_existing_data(rms_psd_list_complete, './data/SR_rms_psd_complete_obscured.pkl') + + mean_strehl_zern = np.mean(strehl_zern_complete, axis=0) + err_strehl_zern = np.std(strehl_zern_complete, axis=0) + mean_strehl_psd = np.mean(strehl_psd_complete, axis=0) + err_strehl_psd = np.std(strehl_psd_complete, axis=0) + + save_variable('./data/SR_RMS_values_obscured.pkl', rmsList) + save_variable('./data/SR_strehl_zern_obscured.pkl', strehl_zern_complete) + save_variable('./data/SR_strehl_psd_obscured.pkl', strehl_psd_complete) + save_variable('./data/SR_rms_zern_complete_obscured.pkl', rms_zern_list_complete) + save_variable('./data/SR_rms_psd_complete_obscured.pkl', rms_psd_list_complete) diff --git a/scripts/SPIE/generate_wfe.py b/scripts/SPIE/generate_wfe.py new file mode 100644 index 0000000..f49edff --- /dev/null +++ b/scripts/SPIE/generate_wfe.py @@ -0,0 +1,88 @@ +import numpy as np +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +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 wavefront errors that will be used for + the ADC prisms during all subsequent analyses. Script made for + the SPIE manuscript. +''' + +def main(): + plt.rc('text', usetex=True) + plt.rc('font', family='serif') + # Initializing MICADO paraxial model + + l = 1.49e-3 + mcd = ParaxialMICADO(n=3001, m=3001, dpix=50) + wfe = WavefrontError(mcd.n, mcd.m) + + # Get x3 and y3 (the coordinates to apply the ADC WFE to): + _, x3, y3 = mcd.propagate(l=l,surface=4, zoom=False, show=False, + print_progress=False) + + # Apply a set of wavefront errors to each ADC surface. + print('\nSetting the wavefront errors...') + wfe_p1 = ADC_surfaces(x3,y3) + wfe_p2 = ADC_surfaces(x3,y3) + wfe_p3 = ADC_surfaces(x3,y3) + wfe_p4 = ADC_surfaces(x3,y3) + + wfe_p1.generate_ADC_wfe(zernRMS=20e-6, psdRMS=15e-6) + wfe_p2.generate_ADC_wfe(zernRMS=20e-6, psdRMS=15e-6) + wfe_p3.generate_ADC_wfe(zernRMS=20e-6, psdRMS=15e-6) + wfe_p4.generate_ADC_wfe(zernRMS=20e-6, psdRMS=15e-6) + + mask = np.sqrt(wfe_p1.xz**2 + wfe_p1.yz**2) > 1 + + # Plot P1 + + surf_p1 = wfe_p1.latest_surface.copy() + surf_p1[mask]=np.nan + figp1, axp1 = plot_real_valued_results(surf_p1[850:2150,850:2150]*1e6, x3[850:2150,850:2150], y3[850:2150,850:2150], + colormap='jet', labeltext='Surface height (nm)', fig_size=None, show=False) + plt.savefig('./figures/ADC_tWFE_p1_zern_20nm_psd_15nm.pdf', bbox_inches='tight', transparent=True) + plt.savefig('./figures/ADC_tWFE_p1_zern_20nm_psd_15nm.png', bbox_inches='tight', transparent=True) + plt.close() + + # Plot P2 + surf_p2 = wfe_p2.latest_surface.copy() + surf_p2[mask]=np.nan + figp2, axp2 = plot_real_valued_results(surf_p2[850:2150,850:2150]*1e6, x3[850:2150,850:2150], y3[850:2150,850:2150], + colormap='jet', labeltext='Surface height (nm)', fig_size=None, show=False) + plt.savefig('./figures/ADC_tWFE_p2_zern_20nm_psd_15nm.pdf', bbox_inches='tight', transparent=True) + plt.savefig('./figures/ADC_tWFE_p2_zern_20nm_psd_15nm.png', bbox_inches='tight', transparent=True) + plt.close() + + # Plot P3 + surf_p3 = wfe_p3.latest_surface.copy() + surf_p3[mask]=np.nan + figp3, axp3 = plot_real_valued_results(surf_p3[850:2150,850:2150]*1e6, x3[850:2150,850:2150], y3[850:2150,850:2150], + colormap='jet', labeltext='Surface height (nm)', fig_size=None, show=False) + plt.savefig('./figures/ADC_tWFE_p3_zern_20nm_psd_15nm.pdf', bbox_inches='tight', transparent=True) + plt.savefig('./figures/ADC_tWFE_p3_zern_20nm_psd_15nm.png', bbox_inches='tight', transparent=True) + plt.close() + + # Plot P4 + surf_p4 = wfe_p4.latest_surface.copy() + surf_p4[mask]=np.nan + figp4, axp4 = plot_real_valued_results(surf_p4[850:2150,850:2150]*1e6, x3[850:2150,850:2150], y3[850:2150,850:2150], + colormap='jet', labeltext='Surface height (nm)', fig_size=None, show=False) + plt.savefig('./figures/ADC_tWFE_p4_zern_20nm_psd_15nm.pdf', bbox_inches='tight', transparent=True) + plt.savefig('./figures/ADC_tWFE_p4_zern_20nm_psd_15nm.png', bbox_inches='tight', transparent=True) + plt.close() + + # Get the wavefront errors and save them for later use. + save_variable('./data/ADC_tWFE_p1_zern_20nm_psd_15nm.pkl', wfe_p1.latest_surface) + save_variable('./data/ADC_tWFE_p2_zern_20nm_psd_15nm.pkl', wfe_p2.latest_surface) + save_variable('./data/ADC_tWFE_p3_zern_20nm_psd_15nm.pkl', wfe_p3.latest_surface) + save_variable('./data/ADC_tWFE_p4_zern_20nm_psd_15nm.pkl', wfe_p4.latest_surface) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/scripts/SPIE/geometric_distortions.py b/scripts/SPIE/geometric_distortions.py new file mode 100644 index 0000000..6b7841d --- /dev/null +++ b/scripts/SPIE/geometric_distortions.py @@ -0,0 +1,162 @@ +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') +from fourierOptics import * + +''' + 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_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye = load_variable('./data/GD_obscured_PSF_ycentroid_expected_zern_psd.pkl') + + # Coordinates with WFEs + xf = load_variable('./data/GD_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yf = load_variable('./data/GD_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.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/GD_full_field_zern_psd.pdf', bbox_inches='tight', transparent=True) + plt.savefig('./figures/GD_full_field_zern_psd.png', bbox_inches='tight', transparent=True) + plt.tight_layout() + plt.show() + + +if __name__ == '__main__': + # plt.rc('text', usetex=True) + # plt.rc('font', family='serif') + + # Setup of simulation + l = 1.49e-3 + 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 = hxList.flatten() + hyList = hyList.flatten() + + # Loading in the WFEs + wfe_p1 = load_variable('./data/ADC_tWFE_p1_zern_20nm_psd_15nm.pkl') + wfe_p2 = load_variable('./data/ADC_tWFE_p2_zern_20nm_psd_15nm.pkl') + wfe_p3 = load_variable('./data/ADC_tWFE_p3_zern_20nm_psd_15nm.pkl') + wfe_p4 = load_variable('./data/ADC_tWFE_p4_zern_20nm_psd_15nm.pkl') + + # 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) ) + + # 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_ADC.p1.add_wavefront_error(wfe_p1) + mcd.ELT_ADC.p2.add_wavefront_error(wfe_p2) + mcd.ELT_ADC.p3.add_wavefront_error(wfe_p3) + mcd.ELT_ADC.p4.add_wavefront_error(wfe_p4) + + + # 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 + + + # 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_obscured_PSF_xcentroid_expected_zern_psd.pkl', x_expected) + save_variable('./data/GD_obscured_PSF_ycentroid_expected_zern_psd.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_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl', x_centroid) + save_variable('./data/GD_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl', y_centroid) + + plot_distortions() + + + diff --git a/scripts/SPIE/plot_centroid_shift.py b/scripts/SPIE/plot_centroid_shift.py new file mode 100644 index 0000000..3582c11 --- /dev/null +++ b/scripts/SPIE/plot_centroid_shift.py @@ -0,0 +1,78 @@ +import numpy as np +import os +import sys +import matplotlib.pyplot as plt +import matplotlib.cm as cm +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to plot the results of the centroid shift analysis. + Script made for the SPIE manuscript. +''' + + +def plot_colourline(x,y,c, ax): + ''' + Source: https://stackoverflow.com/questions/36505587/color-line-by-third-variable-python + Function to colour a line according to a thrid variable. + In our case, this will be the Strehl ratio + ''' + c = cm.plasma_r(c) + for i in np.arange(len(x)-1): + ax.plot([x[i],x[i+1]], [y[i],y[i+1]], c=c[i]) + return + + + + +if __name__ == '__main__': + plt.rc('font', family='serif') + plt.rc('text', usetex=True) + + # Load in the results from centroid_shift_analysis.py + rms_zern_nm, r_airy_zern_mas, r_err_zern_airy_mas, strehl_list_zern = load_variable('./data/centroid_shift_zernikes.pkl') + rms_psd_nm, r_airy_psd_mas, r_err_psd_airy_mas, strehl_list_psd = load_variable('./data/centroid_shift_psd.pkl') + + # Convert the errors to logarithmic space, assuming small errors. + log_err_zern_airy_mas = 0.434 * r_err_zern_airy_mas / r_airy_zern_mas + log_err_psd_airy_mas = 0.434 * r_err_psd_airy_mas / r_airy_psd_mas + + + # Plot the resulting WFE scaling relation + fig = plt.figure(figsize=(9,7)) + ax1 = fig.add_subplot(111) + ax1.errorbar((0.8138)*rms_zern_nm, r_airy_zern_mas, yerr=log_err_zern_airy_mas, color='k', + linestyle='none', capsize=3, + markeredgewidth=1) + ax1.errorbar((0.8138)*rms_psd_nm, r_airy_psd_mas, yerr=log_err_psd_airy_mas, color='k', + linestyle='none', capsize=3, + markeredgewidth=1) + ax1.scatter((0.8138)*rms_zern_nm, r_airy_zern_mas, c=strehl_list_zern, marker='v', label='Zernikes', cmap='plasma_r') + + im = ax1.scatter((0.8138)*rms_psd_nm, r_airy_psd_mas, c=strehl_list_psd, marker='o', label='PSD', cmap='plasma_r') + + # We colour the lines by their Strehl ratio + plot_colourline((0.8138)*rms_zern_nm, r_airy_zern_mas, strehl_list_zern, ax=ax1) + plot_colourline((0.8138)*rms_psd_nm, r_airy_psd_mas, strehl_list_psd, ax=ax1) + ax1.set_xscale('log') + ax1.set_yscale('log') + ax1.set_xlabel('$\\sigma_{\\mathrm{WFE}}$ (nm)', fontsize=18) + ax1.set_ylabel('Centroid shift (mas)', fontsize=18) + ax1.set_xlim([np.min(rms_zern_nm)*(0.8138), np.max(rms_zern_nm)*(0.8138)]) + ax1.legend(fontsize=18) + ax1.grid(True, alpha=0.3, which='major') + ax1.grid(True, alpha=0.1, which='minor') + ax1.tick_params(labelsize=18) + cb = fig.colorbar(im, ax=ax1) + cb.set_label(label='Strehl', size=18) + cb.ax.tick_params(labelsize=18) + plt.savefig('./figures/WFE_scaling_relation_obscured_airy.png', bbox_inches='tight', transparent=True) + plt.savefig('./figures/WFE_scaling_relation_obscured_airy.pdf', bbox_inches='tight', transparent=True) + plt.tight_layout() + plt.show() + + \ No newline at end of file diff --git a/scripts/SPIE/plot_geometric_distortions.py b/scripts/SPIE/plot_geometric_distortions.py new file mode 100644 index 0000000..00ebb41 --- /dev/null +++ b/scripts/SPIE/plot_geometric_distortions.py @@ -0,0 +1,65 @@ +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to plot the results of the geometric distortion analysis. + Script made for the SPIE manuscript. +''' + +if __name__ == '__main__': + plt.rc('font', family='serif') + plt.rc('text', usetex=True) + + # Full field + xe = load_variable('../SPIE/data/GD full field (Obscured Airy zern psd)/GD_obscured_PSF_xcentroid_expected_zern_psd.pkl') + ye = load_variable('../SPIE/data/GD full field (Obscured Airy zern psd)/GD_obscured_PSF_ycentroid_expected_zern_psd.pkl') + xp = load_variable('../SPIE/data/GD full field (Obscured Airy zern psd)/GD_obscured_PSF_xcentroid_found_with_wfe_airy_zern_psd.pkl') + yp = load_variable('../SPIE/data/GD full field (Obscured Airy zern psd)/GD_obscured_PSF_ycentroid_found_with_wfe_airy_zern_psd.pkl') + + # Shape of the arrays + n = xe.shape[0] + + # Reshaping of the arrays so that they can be used for the quiver plot. + X = xe.flatten() + Y = ye.flatten() + dX = xp-xe + dY = yp-ye + + dx = dX.flatten() + dy = dY.flatten() + + # Two values are comparatively (extremely large), they are set to 0 + # to be able to see the other values. + R = (np.sqrt(dx**2 + dy**2) * 100.13151872258219 * 1000) # Magnitude of the shift in micro arcseconds + print('RMS of the PSF shifts: {:.4f} uas.'.format(np.sqrt(np.mean((R)**2)))) + + + # Dummy figure + im = plt.imshow(np.reshape(R, (n,n))) + plt.close() + + # Actual plotting of the figure + fig = plt.figure(figsize=(9,7)) + 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=1) + ax.quiver(X,Y,dx,dy,R, alpha=1) + ax.set_aspect('equal', adjustable='box') + ax.set_xlabel('x (mm)', fontsize=18) + ax.set_ylabel('y (mm)', fontsize=18) + ax.set_xlim([-340,340]) + ax.set_ylim([-340,340]) + ax.tick_params(labelsize=18) + cb = fig.colorbar(im, ax=ax) + cb.set_label(label='PSF centroid shift ($\\mu$as)', size=18) + cb.ax.tick_params(labelsize=18) + plt.savefig('./figures/GD_full_field_obscured_Airy_psd_zern.pdf', bbox_inches='tight', transparent=True) + plt.savefig('./figures/GD_full_field_obscured_Airy_psd_zern.png', bbox_inches='tight', transparent=True) + plt.tight_layout() + plt.show() \ No newline at end of file diff --git a/scripts/SPIE/plot_strehl.py b/scripts/SPIE/plot_strehl.py new file mode 100644 index 0000000..83eb40a --- /dev/null +++ b/scripts/SPIE/plot_strehl.py @@ -0,0 +1,79 @@ +import numpy as np +import os +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to generate the plots of the relation between the + Strehl ratio and the RMS of the Wavefront Error. A comparison + with the Marechal approximation is made. + Script made for the SPIE manuscript. +''' + +def expected_strehl2(rms, l): + return np.exp( -(2*np.pi*rms/l)**2 ) + + +def append_to_existing_data(A, filename): + if os.path.exists(filename): + A_o = load_variable(filename) + A = np.append(A, A_o) + return A + + +def plot_some_stuff(): + fig = plt.figure(figsize=(9,7)) + ax = fig.add_subplot(111) + ax2 = ax.twiny() + + rmsValues = rmsList[:25] + for i in range(len(rms_zern_list_complete)): + ax.plot((0.8138)*rms_zern_list_complete[i]*1e6, strehl_zern_complete[i], 'y-', + alpha=0.1, linewidth=.5) + ax.plot((0.8138)*rms_psd_list_complete[i]*1e6, strehl_psd_complete[i], 'b-', + alpha=0.1, linewidth=.5) + + ax.plot(np.linspace(1e-6,1e-3,100)*1e6, expected_strehl2(np.linspace(1e-6,1e-3,100),l),'k--', + label='$S=\\exp[-(2\\pi\\sigma/\\lambda)^2]$') + ax.errorbar((0.8138)*rmsValues*1e6, mean_strehl_zern, yerr=err_strehl_zern, color='y', + linestyle='-', label='Zernikes', capsize=3, markeredgewidth=1)#, + # linewidth=1) + ax.errorbar((0.8138)*rmsValues*1e6, mean_strehl_psd, yerr=err_strehl_psd, color='b', + linestyle='-', label='PSD', capsize=3, markeredgewidth=1)#, + # linewidth=1) + ax.set_xlabel('$\\sigma_{\\mathrm{WFE}}$ (nm)', fontsize=18) + ax.set_ylabel('Strehl', fontsize=18) + ax.set_ylim([0,1]) + ax.set_xlim([np.min(rmsValues)*1e6*(0.8138), np.max(rmsValues)*1e6*(0.8138)]) + ax2.set_xlim([np.min(rmsValues*(0.8138)/l), np.max(rmsValues*(0.8138)/l)]) + ax2.set_xlabel('$\\sigma/\\lambda$', fontsize=18) + ax.legend(fontsize=18) + ax.grid(True, alpha=0.3) + ax.tick_params(labelsize=18) + ax2.tick_params(labelsize=18) + plt.tight_layout() + plt.savefig('./figures/SR_relation_WFE_on_ADC_p1_obscured.png', bbox_inches='tight', transparent=True) + plt.savefig('./figures/SR_relation_WFE_on_ADC_p1_obscured.pdf', bbox_inches='tight', transparent=True) + plt.show() + +if __name__ == '__main__': + plt.rc('font', family='serif') + plt.rc('text', usetex=True) + plt.rc('axes', labelsize='large') + + l = 1.49e-3 + rmsList = load_variable('./data/SR_RMS_values_obscured.pkl') + strehl_zern_complete = load_variable('./data/SR_strehl_zern_obscured.pkl') + strehl_psd_complete = load_variable('./data/SR_strehl_psd_obscured.pkl') + rms_zern_list_complete = load_variable('./data/SR_rms_zern_complete_obscured.pkl') + rms_psd_list_complete = load_variable('./data/SR_rms_psd_complete_obscured.pkl') + + mean_strehl_zern = np.mean(strehl_zern_complete, axis=0) + err_strehl_zern = np.std(strehl_zern_complete, axis=0) + mean_strehl_psd = np.mean(strehl_psd_complete, axis=0) + err_strehl_psd = np.std(strehl_psd_complete, axis=0) + + plot_some_stuff() diff --git a/scripts/SPIE/plot_wfe_analysis.py b/scripts/SPIE/plot_wfe_analysis.py new file mode 100644 index 0000000..21156a1 --- /dev/null +++ b/scripts/SPIE/plot_wfe_analysis.py @@ -0,0 +1,97 @@ +import numpy as np +import os +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to plot the results of the analysis where the shift + of the PSF centroid as a function of the RMS value of the + WFE was investigated. + Script made for the SPIE manuscript. +''' + + +if __name__ == '__main__': + plt.rc('font', family='serif') + plt.rc('text', usetex=True) + plt.rc('axes', labelsize='large') + + # Optical system constant + platescale = 100.13151872258219 + + ############################### + ### Lowpass filter analysis ### + ############################### + + # Load in the data + r_airy, r_err_airy = load_variable('./data/lowpass_frequency_filter_airy_distance.pkl') + x0_airy, y0_airy, x0_airy_err, y0_airy_err = load_variable('./data/lowpass_frequency_filter_airy_locations.pkl') + central_frequencies = np.asarray([0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 20, 25, + 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 120, 140, 160, 180, 200, + 220, 250, 280, 310, 340, 370, 400, 450, 500]) + + r_airy_mas = r_airy * platescale + r_err_airy_mas = r_err_airy * platescale + + # Plotting the centroid shift for all bandpasses + fig = plt.figure(figsize=(9,7)) + ax1 = fig.add_subplot(111) + ax1.scatter(central_frequencies, r_airy_mas, edgecolor='k', facecolors='none', alpha=0.7) + ax1.set_xscale('log') + ax1.set_xlabel('Cut-off Frequency of the lowpass filter ($\\rho/\\rho_c$)', fontsize=18) + ax1.set_ylabel('Centroid shift (mas)', fontsize=18) + ax1.tick_params(labelsize=18) + ax1.grid(True, alpha=0.3, which='major') + ax1.grid(True, alpha=0.1, which='minor') + ax1.set_xlim([0.3,550]) + plt.tight_layout() + plt.savefig('./figures/wfe_frequency_analysis/lowpass_spatial_freq_distance.png', bbox_inches='tight', transparent=True) + plt.savefig('./figures/wfe_frequency_analysis/lowpass_spatial_freq_distance.pdf', bbox_inches='tight', transparent=True) + + + + + + ################################ + ### Bandpass filter analysis ### + ################################ + + # Load in the data + r_airy, r_err_airy = load_variable('./data/bandpass_frequency_filter_airy_distance.pkl') + x0_airy, y0_airy, x0_airy_err, y0_airy_err = load_variable('./data/bandpass_frequency_filter_airy_locations.pkl') + central_frequencies = load_variable('./data/bandpass_frequency_filter_central_frequencies.pkl') + print(central_frequencies, len(central_frequencies), len(r_airy)) + r_airy_mas = r_airy * platescale + r_err_airy_mas = r_err_airy * platescale + + + # Fit a line in Logarithmic space + log_cf = np.log10(central_frequencies) + log_r = np.log10(r_airy_mas) + coefficients = np.polyfit(log_cf[2:], log_r[2:], 1) + best_fit_fnc = np.poly1d(coefficients) + best_fit_y = best_fit_fnc(log_cf) + + + # Plotting the centroid shift for all bandpasses + fig = plt.figure(figsize=(9,7)) + ax1 = fig.add_subplot(111) + ax1.scatter(central_frequencies, r_airy_mas, edgecolor='k', facecolors='none', alpha=0.7) + ax1.plot(central_frequencies[2:], 10**best_fit_y[2:], 'b-', alpha=0.5, label='Best fit: $y = 0.030 x^{-0.602}$') + ax1.set_xscale('log') + ax1.set_yscale('log') + ax1.set_xlabel('Central Frequency of the bandpass filter ($\\rho/\\rho_c$)', fontsize=18) + ax1.set_ylabel('Centroid shift (mas)', fontsize=18) + ax1.tick_params(labelsize=18) + ax1.grid(True, alpha=0.3, which='major') + ax1.grid(True, alpha=0.1, which='minor') + ax1.set_xlim([0.3,550]) + ax1.set_ylim([10**-6, 10**0]) + ax1.legend(fontsize=18) + plt.tight_layout() + plt.savefig('./figures/wfe_frequency_analysis/bandpass_spatial_freq_distance.png', bbox_inches='tight', transparent=True) + plt.savefig('./figures/wfe_frequency_analysis/bandpass_spatial_freq_distance.pdf', bbox_inches='tight', transparent=True) + plt.show() \ No newline at end of file diff --git a/scripts/SPIE/wavefront_error_frequency_analysis.py b/scripts/SPIE/wavefront_error_frequency_analysis.py new file mode 100644 index 0000000..8d123c4 --- /dev/null +++ b/scripts/SPIE/wavefront_error_frequency_analysis.py @@ -0,0 +1,252 @@ +import numpy as np +import os +import sys +sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics') +from fourierOptics import * + +''' + Author: Joost van den Born, PhD student at University of Groningen + Contact: born@astro.rug.nl + Description: Script to analyse the shift of the PSF centroid as a + function of the RMS value of the WFE. + Script made for the SPIE manuscript. +''' + + +if __name__ == '__main__': + plt.rc('font', family='serif') + plt.rc('axes', labelsize='large') + + # Pick a wavelength + l = 1.49e-3 + input_field_location = './data/ADC_input_field_1001x1001_1.49_onaxis_circular.pkl' + + # Initializing MICADO paraxial model + mcd = ParaxialMICADO(n=1001, m=1001, dpix=160, entrance_pupil='open') + + # Propagation of the unperturbed system. + Exyz_og, xp, yp = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + + # Generate the wavefront error that will be used. + wfe = WavefrontError(mcd.n,mcd.m) + wfe.generate_height_map(5e-5, 2, 1) # 50 nm RMS for standard, power slope 2, cut-off 1 + h_data = wfe.data + + # Plotting the WFE, PSD and image planes + oC = [np.min(xp), np.max(xp), np.min(yp), np.max(yp)] + + # Mask out areas where the ADC prism isn't + _, xadc, yadc = mcd.propagate(l=l,surface=3, zoom=False, show=False, + print_progress=True, hx=0, hy=0) + mask = np.sqrt(xadc**2 + yadc**2) > 62.5 + + # Setting up empty lists for fitting process + popt_airy_list, popt_cust_list = [], [] + pcov_airy_list, custom_var_list = [], [] + Exyzlist, PSDlist, surfacelist = [], [], [] + central_frequencies = [] + + # Bins in LogSpace + x = mcd.x + flist = np.asarray([0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 20, 25, + 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 120, 140, 160, 180, 200, + 220, 250, 280, 310, 340, 370, 400, 450, 500]) + plt.imshow(wfe.data*1e6, cmap='jet');plt.colorbar();plt.show() + + + ############################### + ### Lowpass filter analysis ### + ############################### + + # Do the analysis for each frequency in flist + for fmax in flist: + print('Lowpass filter - Freq.: {} - Sample {:02d}/{:02d}'.format(fmax, flist.tolist().index(fmax) + 1, len(flist) ) ) + if fmax == flist[0]: + wfe.set_lowpass_filter(fmax) + else: + i = flist.tolist().index(fmax) - 1 + wfe.set_lowpass_filter(flist[i]) + + scal = wfe.filteredData + + # Adding the wavefront error + mcd.ELT_ADC.p1.add_wavefront_error(scal) + + # Performing the propagation using Fresnel approximation to the focal plane. + Exyz_psd,xp,yp = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + + # Set up fitting object + psfFit_psd = PSFFitter(xp.real,yp.real,np.abs(Exyz_psd)**2) + + # Fit airy pattern and coupling integral + amplitude_guess = np.max(np.abs(Exyz_psd))**2 + airyFitted, popt_airy, pcov_airy = psfFit_psd.fit_airy_PSF(initialGuess=(amplitude_guess, 0.0, 0.0, 0.0, 0.03)) + popt_airy_list.append(popt_airy) + pcov_airy_list.append([ np.sqrt(pcov_airy[i][i]) for i in range(len(pcov_airy)) ]) + + # Other properties + original_PSD_curve = wfe.PSD2D[int((mcd.n-1)/2),int((mcd.n-1)/2):] + wfe.calculate_1D_PSD(filteredDataFlag=True) + PSDlist.append(wfe.filteredPSD2D[int((mcd.n-1)/2),int((mcd.n-1)/2):]) + surfacelist.append(wfe.filteredData) + + # Plot of the surface profile + fig2 = plt.figure(figsize=(12,3)) + axi1 = fig2.add_subplot(1, 3, 1) + wavefrontmap = wfe.filteredData.copy()*1e6 + wavefrontmap[mask] = np.nan + imi1 = axi1.imshow(wavefrontmap, cmap='jet', origin='lower', extent=[xadc.min(), xadc.max(), yadc.min(), yadc.max()]) + axi1.set_xlabel('x (mm)') + axi1.set_ylabel('y (mm)') + axi1.set_xlim([-65,65]) + axi1.set_ylim([-65,65]) + fig2.colorbar(imi1, ax=axi1, label='Surface height (nm)') + + # Plot of the PSD + axi2 = fig2.add_subplot(1, 3, 2) + axi2.loglog(original_PSD_curve, 'k-') + axi2.loglog(PSDlist[-1], 'y-') + axi2.grid(which='major', alpha=0.3) + axi2.grid(which='minor', alpha=0.1) + axi2.set_xlabel('Normalized Radial Frequency ($\\rho/\\rho_c$)') + axi2.set_ylabel('PSD (mm$^4$)') + + # Plot of the PSF + axi3 = fig2.add_subplot(1, 3, 3) + differential_amplitude = np.abs(Exyz_og) - np.abs(Exyz_psd) + if (differential_amplitude == 0).all(): + norm_par = None + else: + norm_par = LogNorm() + imi3 = axi3.imshow(differential_amplitude, extent=oC, cmap='gray', origin='lower', norm=norm_par) + fig2.colorbar(imi3, ax=axi3, label='Differential amplitude') + axi3.set_xlim([-1,1]) + axi3.set_ylim([-1,1]) + axi3.set_xlabel('x (mm)') + axi3.set_ylabel('y (mm)') + plt.tight_layout() + plt.savefig('./figures/wfe_frequency_analysis/lowpass_spatial_freq_{:.0f}.png'.format(flist.tolist().index(fmax) + 1), bbox_inches='tight', transparent=True) + plt.close() + + + # Optimum fit parameters for all scaled surface errors + amplitude_airy, x0_airy, y0_airy, R_airy, offset_airy = np.asarray(popt_airy_list).T + amplitude_airy_err, x0_airy_err, y0_airy_err, R_airy_err, offset_airy_err = np.asarray(pcov_airy_list).T + + # Distance from original centroid at (0,0) + r_airy = np.sqrt(x0_airy**2+y0_airy**2) + r_err_airy = np.sqrt( (x0_airy**2*x0_airy_err**2+y0_airy**2*y0_airy_err**2) / (x0_airy**2+y0_airy**2)) + r_airy_mas = r_airy * mcd.platescale * 1000 + r_err_airy_mas = r_err_airy * mcd.platescale * 1000 + + save_variable('./data/lowpass_frequency_filter_airy_distance.pkl', (r_airy, r_err_airy)) + save_variable('./data/lowpass_frequency_filter_airy_locations.pkl', (x0_airy, y0_airy, x0_airy_err, y0_airy_err)) + print('Platescale in mas/mm: ', mcd.platescale * 1000) + + + + ################################ + ### Bandpass filter analysis ### + ################################ + + # Setting up empty lists for fitting process + popt_airy_list, popt_cust_list = [], [] + pcov_airy_list, custom_var_list = [], [] + Exyzlist, PSDlist, surfacelist = [], [], [] + central_frequencies = [] + + # Do the analysis for each frequency in flist + for fmax in flist: + print('Bandpass - Freq.: {} - Sample {:02d}/{:02d}'.format(fmax, flist.tolist().index(fmax) + 1, len(flist) ) ) + if fmax == flist[0]: + # Use a lowpass filter only for the very first element in flist + wfe.set_lowpass_filter(fmax) + central_frequencies.append(fmax) + else: + # Set a bandpass between the current and previous frequency element + i = flist.tolist().index(fmax) - 1 + wfe.set_bandpass_filter(flist[i], fmax) + central_frequencies.append(np.mean([flist[i],fmax])) + + + # Adding the wavefront error + scal = wfe.filteredData + mcd.ELT_ADC.p1.add_wavefront_error(scal) + + # Performing the propagation using Fresnel approximation to the focal plane. + Exyz_psd,xp,yp = mcd.propagate(l=l,surface=-1, zoom=False, show=False, + print_progress=False, hx=0, hy=0, + load_up_to_adc=input_field_location) + + # Set up fitting object + psfFit_psd = PSFFitter(xp.real,yp.real,np.abs(Exyz_psd)**2) + + # Fit airy pattern and coupling integral + amplitude_guess = np.max(np.abs(Exyz_psd))**2 + airyFitted, popt_airy, pcov_airy = psfFit_psd.fit_airy_PSF(initialGuess=(amplitude_guess, 0.0, 0.0, 0.0, 0.03)) + popt_airy_list.append(popt_airy) + pcov_airy_list.append([ np.sqrt(pcov_airy[i][i]) for i in range(len(pcov_airy)) ]) + + # Other properties + original_PSD_curve = wfe.PSD2D[int((mcd.n-1)/2),int((mcd.n-1)/2):] + wfe.calculate_1D_PSD(filteredDataFlag=True) + PSDlist.append(wfe.filteredPSD2D[int((mcd.n-1)/2),int((mcd.n-1)/2):]) + surfacelist.append(wfe.filteredData) + + # Plot of the surface profile + fig2 = plt.figure(figsize=(12,3)) + axi1 = fig2.add_subplot(1, 3, 1) + wavefrontmap = wfe.filteredData.copy()*1e6 + wavefrontmap[mask] = np.nan + imi1 = axi1.imshow(wavefrontmap, cmap='jet', origin='lower', extent=[xadc.min(), xadc.max(), yadc.min(), yadc.max()]) + axi1.set_xlabel('x (mm)') + axi1.set_ylabel('y (mm)') + axi1.set_xlim([-65,65]) + axi1.set_ylim([-65,65]) + fig2.colorbar(imi1, ax=axi1, label='Surface height (nm)') + + # Plot of the PSD + axi2 = fig2.add_subplot(1, 3, 2) + axi2.loglog(original_PSD_curve, 'k-') + axi2.loglog(PSDlist[-1], 'y-') + axi2.grid(which='major', alpha=0.3) + axi2.grid(which='minor', alpha=0.1) + axi2.set_xlabel('Normalized Radial Frequency ($\\rho/\\rho_c$)') + axi2.set_ylabel('PSD (mm$^4$)') + + # Plot of the PSF + axi3 = fig2.add_subplot(1, 3, 3) + differential_amplitude = np.abs(Exyz_og) - np.abs(Exyz_psd) + if (differential_amplitude == 0).all(): + norm_par = None + else: + norm_par = LogNorm() + imi3 = axi3.imshow(differential_amplitude, extent=oC, cmap='gray', origin='lower', norm=norm_par) + fig2.colorbar(imi3, ax=axi3, label='Differential amplitude') + axi3.set_xlim([-1,1]) + axi3.set_ylim([-1,1]) + axi3.set_xlabel('x (mm)') + axi3.set_ylabel('y (mm)') + plt.tight_layout() + plt.savefig('./figures/wfe_frequency_analysis/bandpass_spatial_freq_{:.0f}.png'.format(flist.tolist().index(fmax) + 1), bbox_inches='tight', transparent=True) + plt.close() + + # Optimum fit parameters for all scaled surface errors + amplitude_airy, x0_airy, y0_airy, R_airy, offset_airy = np.asarray(popt_airy_list).T + amplitude_airy_err, x0_airy_err, y0_airy_err, R_airy_err, offset_airy_err = np.asarray(pcov_airy_list).T + + # Distance from original centroid at (0,0) + r_airy = np.sqrt(x0_airy**2+y0_airy**2) + r_err_airy = np.sqrt( (x0_airy**2*x0_airy_err**2+y0_airy**2*y0_airy_err**2) / (x0_airy**2+y0_airy**2)) + r_airy_mas = r_airy * mcd.platescale * 1000 + r_err_airy_mas = r_err_airy * mcd.platescale * 1000 + + save_variable('./data/bandpass_frequency_filter_airy_distance.pkl', (r_airy, r_err_airy)) + save_variable('./data/bandpass_frequency_filter_airy_locations.pkl', (x0_airy, y0_airy, x0_airy_err, y0_airy_err)) + save_variable('./data/bandpass_frequency_filter_central_frequencies.pkl', central_frequencies) + print('Platescale in mas/mm: ', mcd.platescale * 1000) +