Added off-axis propagation, SPIE scripts

.gitignore

@@ -103,3 +103,7 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+# Exclude SPIE/Data en SPIE/Figures folders
+scripts/SPIE/data
+scripts/SPIE/figures

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]')

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))	
+		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.'
+		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:

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):
 	'''
@@ -91,3 +144,196 @@ def load_variable(filename):
 	with open(filename, 'rb') as f:
 		variable = pickle.load(f)
 	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

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 

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;

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()
+	

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()
+

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))

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()
+	

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)

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()

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()
+	
+
+

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()
+
+	

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()

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()

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()

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)
+