Added Atmospheric Seeing

fourierOptics/__init__.py

@@ -8,3 +8,4 @@ from .psfFitter import *
 from .wavefrontErrors import *
 from .misc import *
 from .micado import *
+from .atmosphere import *

fourierOptics/atmosphere.py

@@ -0,0 +1,135 @@
+import numpy as np
+import scipy.special as spc
+from .misc import *
+
+'''
+	Various functions to simulate a seeing limited PSF. These functions are
+	quite simple and not suited to extensive atmospheric simulations. Nonetheless,
+	they are useful to inspect the impact of atmospheric seeing.
+
+	Author: 	Joost van den Born
+	Contact: 	born@astron.nl
+'''
+
+def seeing(wvl, r0):
+	'''
+		Calculates the seeing in arcseconds at the given wavelength
+		
+		Parameters
+		----------
+		wvl 	- 	Float or Array
+			Wavelegth of light in millimeters
+		r0 		- 	Float
+			Fried length in meter. See also fried_r0()
+
+		Returns
+		-------
+		sng 	- 	Float or Array
+			Seeing in arcseconds
+	'''
+	wvl *= 1e-3 # convert wvl from mm to meter
+	sng = 3600*180/np.pi*0.976*wvl/r0
+	return sng
+
+
+def fried_r0(wvl, seeing):
+	'''
+		Returns Fried's parameter value in meters
+		
+		Parameters
+		----------
+		wvl 	- 	Float or Array
+			Wavelegth of light in millimeters
+		seeing 	- 	Float or Array
+			Seeing in arcseconds
+
+		Returns
+		-------
+		r0 		- 	Float
+			Fried length in meter. See also fried_r0()
+	'''
+	wvl *= 1e-3 # convert wvl from mm to meter
+	r0 = 0.976 * wvl / seeing * 360/(2*np.pi) * 3600
+	return r0
+
+
+def airmass(z):
+	'''
+		Converts zenith angle in degrees to airmass
+	'''
+	return 1/np.cos(np.deg2rad(z))
+
+
+def spectrum(k, r0, L0=25):        
+	'''
+		Computes the atmospheric phase power spectrum density from the
+		spatial frequency f. Based on a Von Karman spectrum.
+		
+		Parameters
+		----------
+		k 		- 	Array
+			k-space spatial frequencies.
+		r0 		- 	Float
+			Fried length in meter. See also fried_r0().
+		L0 		- 	Float (default=25)
+			Maximum correlation length in meters. Typical value is 25 meter.
+
+		Returns
+		spectrum 	- Array
+			The Phase Power Spectrum Density values for the given frequencies
+	'''                        
+	cte = (24*spc.gamma(6/5)/5)**(5/6)*(spc.gamma(11/6)**2./(2*np.pi**(11/3)))
+	spectrum = r0**(-5/3)*cte*(k**2 + 1/L0**2)**(-11/6)
+	return spectrum 
+
+
+def add_seeing(wvl, sng, field_tuple):
+	'''
+		Computes the seeing limited PSF for a supplied amplitude
+		distribution.
+		
+		Parameters
+		----------
+		wvl 		-	Float
+			Wavelength in mm corresponding to the supplied PSF
+		sng 		- 	Float
+			Atmospheric seeing in arcseconds.
+		field_tuple - Tuple 
+			Contains the PSF information (amplitude array, x coordinates, y coordinates)
+		
+		Returns:
+		psf_seeing 	- Array
+			2D array containing the new seeing limited PSF information.
+
+	'''
+	r01   = fried_r0(wvl, seeing=sng) * 1e3
+	A_psf, x_psf, y_psf = field_tuple
+	n,m = np.shape(x_psf)
+	fx = x_psf / (3200*wvl)
+	fy = y_psf / (3200*wvl)
+	dk = find_pixel_scale(fx)
+
+	# von-Karman spectrum of the atmosphere for the given frequencies in k-space
+	psd   = spectrum(np.hypot(fx.real, fy.real), r0=r01, L0=25*1e3)
+	psd[int(n/2),int(m/2)] = 0
+
+	# Compute the atmospheric OTF
+	# According to Roddier 1981, the atmospheric OTF can be derived from the phase structure function
+	# assuming that the phase is stationary accross the pupil and behaves as a Gaussian variable
+	Bphi = np.fft.fft2(np.fft.ifftshift(psd))*dk**2
+	Dphi = np.fft.ifftshift((2 * (Bphi[0,0] - Bphi))) # ifftshift in Bphi en Dphi!
+	otf_atm = np.exp(-0.5*Dphi)
+	
+	# Compute the OTF of the supplied PSF
+	# We can derive the OTF from the inverse Fourier transform of the PSF
+	psf = np.abs(A_psf)**2
+	otf_psf = np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(psf)))
+	otf_psf = otf_psf/otf_psf.max()
+
+	# Compute the seeing-limited PSF
+	# Multiplication of OTF and then Fourier transform
+	otf  = otf_psf * otf_atm
+	psf_seeing = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(otf)))
+	return psf_seeing
+
+

fourierOptics/fourier.py

@@ -6,7 +6,7 @@ from scipy.optimize import leastsq
 from matplotlib.ticker import FormatStrFormatter
 from matplotlib.colors import LogNorm
 import pyfftw
-
+from .misc import find_pixel_scale
 
 class FourierOptics:
 	def __init__(self, l):
@@ -61,9 +61,7 @@ class FourierOptics:
 		_n, _m = np.shape(A)
 		
 		# Calculate the pixel size from the xi, eta plane
-		_dpix = np.mean(np.diff(xin.real))
-		if _dpix < 1e-99:
-			_dpix = np.mean(np.diff(xin.real))
+		_dpix = find_pixel_scale(xin)
 
 		# A at z=0 is then
 		_E_xy0 		= pyfftw.empty_aligned((_n,_m), dtype='complex128')
@@ -181,11 +179,6 @@ class FourierOptics:
 			A_out_ip 		= A_out_ip_real + A_out_ip_imag*1j
 			A 				= np.rint(A_out_ip_real)
 
-		# if oversample:
-			# _Aos 		= np.zeros((np.shape(A)[0]*oversample, np.shape(A)[1]*oversample))
-			# _Aos[int(len(_Aos)/2-len(A)/2): int(len(_Aos)/2+len(A)/2), int(len(_Aos)/2-len(A)/2): int(len(_Aos)/2+len(A)/2)] = A
-			# A 			= _Aos
-
 		_n, _m 			= np.shape(A)
 		_k 				= 2 * np.pi / self.l
 		
@@ -194,9 +187,7 @@ class FourierOptics:
 		
 		
 		# Calculate the pixel size from the xi, eta plane
-		_dpix = np.mean(np.diff(_xi))
-		if _dpix < 1e-99:
-			_dpix = np.mean(np.diff(_xi.T))
+		_dpix = find_pixel_scale(_xi)
 
 		# Gaussian test: use 8*w0 radius for FFT:
 		if outputCoordinates == 'auto':

fourierOptics/misc.py

@@ -28,6 +28,7 @@ def perform_zero_padding(A, xs, ys):
 	y_margin = int((ys - A.shape[1])/2)
 	return np.pad(A, (x_margin, y_margin))
 
+
 def find_nearest(array, value):
 	'''
 		For a value it finds the nearest element in an array
@@ -39,6 +40,7 @@ def find_nearest(array, value):
 	idx = np.unravel_index(diff, array.shape)
 	return idx
 
+
 def find_maximum(array):
 	'''
 		Finds the index of the maximum of the given array
@@ -46,6 +48,16 @@ def find_maximum(array):
 	idx = np.unravel_index(array.argmax(), array.shape)
 	return idx
 
+def find_pixel_scale(x):
+	''' 
+		Calculates the pixel size from the supplied array. 
+	'''
+	_dpix = np.mean(np.diff(x.real))
+	if _dpix < 1e-99:
+		_dpix = np.mean(np.diff((x.T).real))
+	return _dpix
+
+
 def select_subsection(x, y , outputCoordinates):
 	'''
 		Returns the indices of a subsection of an x,y meshgrid
@@ -56,6 +68,7 @@ def select_subsection(x, y , outputCoordinates):
 	select_4 = find_nearest(y, outputCoordinates[3])[0]
 	return [select_1, select_2, select_3, select_4]
 
+
 def perform_extrapolation(x, y, n, m):
 	'''
 		Performs linear extrapolation of a meshgrid to a new 
@@ -75,7 +88,41 @@ 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'):
+
+def plot_results(E, x, y, logscale=False, fig_size=(8,3), grid=False, show=True, 
+				colormap='hot'):
+	'''
+		Quick function to plot complex valued image data, for example the
+		amplitude field of a PSF. The absolute value of the complex values
+		will be used. 
+		If the amplitude values can be smaller than zero, use 
+		plot_real_valued_results().
+
+		Parameters
+		----------
+		E 		-	2D Array
+			Amplitude or intensity of the field to be shown.
+		x 		- 	Array
+			X-coordinates of the field to be plotted
+		y 		- 	Array
+			Y-coordinates of the field to be plotted
+		logscale - 	Bool (default=False)
+			If True, imshow() will use a logarithmic color scale
+		fig_size - 	Tuple (default=(8,8))
+			Size of the figure
+		grid 	- 	Bool (default=False)
+			If True, the figure will show gridlines.
+		show 	-	Bool (default=True)
+			If True, the figure will be shown before continuing.
+		colormap -  String (default='hot')
+			String, defines the colormap imshow() will use.
+		labeltext - String (default=None)
+			Defines the label to be shown along the colorbar
+
+		Returns
+		-------
+		_fig, _ax - Matplotlib figure and axis object.
+	'''
 	x = x.real
 	y = y.real
 	_fig = plt.figure(figsize=fig_size)
@@ -105,7 +152,38 @@ def plot_results(E, x, y, logscale=False, fig_size=(8,3),grid=False, show=True,
 	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):
+def plot_real_valued_results(E, x, y, logscale=False, fig_size=(8,8), grid=False, show=True, 
+							colormap='hot', labeltext=None):
+	'''
+		Quick function to plot real valued image data, for example the
+		intensity field of a PSF.
+
+		Parameters
+		----------
+		E 		-	2D Array
+			Amplitude or intensity of the field to be shown.
+		x 		- 	Array
+			X-coordinates of the field to be plotted
+		y 		- 	Array
+			Y-coordinates of the field to be plotted
+		logscale - 	Bool (default=False)
+			If True, imshow() will use a logarithmic color scale
+		fig_size - 	Tuple (default=(8,8))
+			Size of the figure
+		grid 	- 	Bool (default=False)
+			If True, the figure will show gridlines.
+		show 	-	Bool (default=True)
+			If True, the figure will be shown before continuing.
+		colormap -  String (default='hot')
+			String, defines the colormap imshow() will use.
+		labeltext - String (default=None)
+			Defines the label to be shown along the colorbar
+
+		Returns
+		-------
+		_fig, _ax - Matplotlib figure and axis object.
+	'''
+
 	x = x.real
 	y = y.real
 	_fig = plt.figure(figsize=fig_size)

fourierOptics/optics.py

@@ -1,5 +1,6 @@
 import numpy as np
 from scipy.optimize import fsolve
+from scipy import ndimage
 from .data import glassCatalog
 
 class LensObject:
@@ -708,6 +709,13 @@ class Prism:
 	def set_prism_orientation(self, theta):
 		self.tilt1x, self.tilt1y = self.tilt_rotation(self.tilt1x_OG, self.tilt1y_OG, theta)
 		self.tilt2x, self.tilt2y = self.tilt_rotation(self.tilt2x_OG, self.tilt2y_OG, theta)
+		if self.wavefrontError:
+			if self.sag_combined is not None:
+				self.sag_combined = ndimage.rotate(self.sag_combined, theta, reshape=False)
+			elif self.sag1 is not None:
+				self.sag1 = ndimage.rotate(self.sag1, theta, reshape=False)
+			elif self.sag2 is not None:
+				self.sag2 = ndimage.rotate(self.sag2, theta, reshape=False)
 		self.thetaPrism = theta
 	
 
@@ -718,8 +726,7 @@ class Prism:
 		return tx_new, ty_new
 
 
-	def get_surface_sag(self, x, y, theta):
-		# self.set_prism_orientation(theta)
+	def get_surface_sag(self, x, y):
 		self.first_surface_sag = y * np.tan(np.deg2rad(self.tilt1x)) + x * np.tan(np.deg2rad(self.tilt1y))
 		self.second_surface_sag = y * np.tan(np.deg2rad(self.tilt2x)) + x * np.tan(np.deg2rad(self.tilt2y))
 		
@@ -766,7 +773,7 @@ class Prism:
 		_n = self.find_refractive_index(l)
 
 		# Retrieve the sag of the lens surfaces for coordinates x,y
-		_A1, _A2 	= self.get_surface_sag(x,y, theta)
+		_A1, _A2 	= self.get_surface_sag(x,y)
 
 		# Masking out the parts that fall outside the lens clear aperture
 		dist_from_center = np.sqrt(x**2 + y**2)

scripts/SPIE/generate_ADC_footprint.py

@@ -49,8 +49,8 @@ if __name__ == '__main__':
 			# 							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.savefig('./figures/ELT-MCD-ADC_footprint_no_wfe_{}_{}.pdf'.format(hx,hy), bbox_inches='tight', transparent=True)
+			plt.savefig('./figures/ELT-MCD-ELT-M1.png', bbox_inches='tight', transparent=True)	
 			plt.show()
 
 

scripts/gratama_seeing_limited_dispersion_test.py

@@ -0,0 +1,172 @@
+import numpy as np
+import sys
+sys.path.insert(0, '/Users/joost/Dropbox/Werk/PhD/Analyses/Python/Analytical/Fourier Optics')
+from fourierOptics import *
+import scipy.special as spc
+from pdb import set_trace
+'''
+def find_pixel_scale(x):
+	""" Calculate the pixel size from the x, y plane """
+	_dpix = np.mean(np.diff(x.real))
+	if _dpix < 1e-99:
+		_dpix = np.mean(np.diff((x.T).real))
+	return _dpix
+
+def seeing(wvl, r0):
+	"""
+		Seeing value in arcsec at wavelength wvl
+		wvl in millimeter
+		r0 in meter
+	"""
+	wvl *= 1e-3 # convert wvl from mm to meter
+	return 3600*180/np.pi*0.976*wvl/r0
+
+def fried_r0(wvl, seeing):
+	"""
+		Returns Fried's parameter value in meters
+		wvl in millimeter
+	"""
+	wvl *= 1e-3 # convert wvl from mm to meter
+	return 0.976*wvl/seeing*206264.8
+
+def airmass(z):
+	"""
+		Converts zenith angle in degrees to airmass
+	"""
+	return 1/np.cos(np.deg2rad(z))
+
+def spectrum(k, r0, L0=25):        
+	"""SPECTRUM Phase power spectrum density computes the phase power
+		spectrum density from the spatial frequency f and an
+		atmosphere object
+		
+		Normal value for L0 = 25.
+	"""                        
+	cte = (24*spc.gamma(6/5)/5)**(5/6)*(spc.gamma(11/6)**2./(2*np.pi**(11/3)))
+	return r0**(-5/3)*cte*(k**2 + 1/L0**2)**(-11/6)
+
+def add_seeing(wvl, sng, field_tuple):
+	"""
+		Computes the seeing limited PSF for a supplied amplitude
+		distribution.
+		
+		Parameters:
+		wvl - wavelength in mm
+		sng - seeing in arcseconds
+		field_tuple - tuple containing (amplitude array, x array, y array)
+		
+		Returns:
+		psf_seeing
+
+	"""
+	r01   = fried_r0(wvl, seeing=sng) * 1e3
+	A_psf, x_psf, y_psf = field_tuple
+	n,m = np.shape(x_psf)
+	fx = x_psf / (3200*wvl)
+	fy = y_psf / (3200*wvl)
+	dk = find_pixel_scale(fx)
+
+	# von-Karman spectrum of the atmosphere for the given frequencies in k-space
+	psd   = spectrum(np.hypot(fx.real, fy.real), r0=r01, L0=25*1e3)
+	psd[int(n/2),int(m/2)] = 0
+
+	# Compute the atmospheric OTF
+	# According to Roddier 1981, the atmospheric OTF can be derived from the phase structure function
+	# assuming that the phase is stationary accross the pupil and behaves as a Gaussian variable
+	Bphi = np.fft.fft2(np.fft.ifftshift(psd))*dk**2
+	Dphi = np.fft.ifftshift((2 * (Bphi[0,0] - Bphi))) # ifftshift in Bphi en Dphi!
+	otf_atm = np.exp(-0.5*Dphi)
+	
+	# Compute the OTF of the supplied PSF
+	# We can derive the OTF from the inverse Fourier transform of the PSF
+	psf = np.abs(A_psf)**2
+	otf_psf = np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(psf)))
+	otf_psf = otf_psf/otf_psf.max()
+
+	# Compute the seeing-limited PSF
+	# Multiplication of OTF and then Fourier transform
+	otf  = otf_psf * otf_atm
+	psf_seeing = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(otf)))
+	return psf_seeing
+'''
+
+
+
+if __name__ == '__main__':
+	l  = 0.5e-3      # mm
+	D_GT = 400          # mm
+	dpix = 1.6
+	n,m  = 1001, 1001
+	x,y  = np.meshgrid(np.linspace(-n/2*dpix, n/2*dpix, n), np.linspace(-m/2*dpix, m/2*dpix, m))
+	wfe = WavefrontError(n, m)
+
+	# Create the Bahtinov mask and set it on the entrance pupil (ELT-M1)
+	GT_MASK = Profiles().bahtinov_mask(D_GT/2, 186.53/2, n, m, dpix=1.6, N=30, oversample=1).astype('complex128')
+	GT_M1 = ParabolicMirror(dx=0, dy=0, name='GT_Paraxial', CA=400, f=3200)
+	GT_platescale = (3600*180/np.pi) / (GT_M1.CA * 8)
+
+	# Set up
+	Elist = []
+	lList = np.linspace(0.398e-3, 0.727e-3, 25).tolist()
+
+	# Dispersion over the band
+	dispMax = 1.5 / 3600 		# maximum dispersion in degrees
+	dispersion = np.linspace(0, dispMax, len(lList)).tolist()
+
+	# Seeing conditions
+	sng = 1.5
+
+
+	# Iterate over the wavelengths
+	for i in range(len(lList)):
+		
+		l = lList[i]
+		print(l)
+		# Generating the tilt to the input field.
+		d = dispersion[i]
+		field_tilt = wfe.generate_tilt(x,y,tx=0, ty=d)
+		
+		# Adding the tilt to the input field
+		Ein = GT_MASK * np.exp(-1j * 2*np.pi/l * field_tilt) * GT_M1.phase_delay(x,y,l, paraxial=True)
+
+		# Propagate the input field throug the paraxial model
+		# Output diffraction limited PSF
+		fou = FourierOptics(l)
+		oC = [-0.6, 0.6, -0.6, 0.6]
+		Eout, xout, yout = fou.fresnel_propagation(Ein, x, y, dz=GT_M1.f, show=False, 
+			outputCoordinates=oC, apply_interpolation=True)
+		psf_GT_DL = (Eout, xout, yout)
+
+		psf_seeing = add_seeing(l, sng, psf_GT_DL)
+
+		Elist.append(psf_seeing)
+		# Elist.append(Eout)
+
+	# Average over all the wavelengths to get the final field.
+	Eca = np.asarray([np.abs(x) for x in Elist]).mean(axis=0) 	# Amplitude 												# Intensity
+	Ecp = np.asarray([np.angle(x) for x in Elist]).mean(axis=0) # Phase
+
+	Eca = Eca/Eca.max() * 64000
+
+	# Adding background level and nosie, estimated from FITS data of the Gratama Telescope
+	Eca = Eca + np.random.normal(2000, 13.5, np.shape(Eca))
+
+	# Convert the scale to sky-angles using the plate scale
+	x = xout * GT_platescale
+	y = yout * GT_platescale
+	
+	# save_variable('GT_psf_with_dispersion.pkl', (Eca, x, y))
+
+	# Plotting the results
+	fig = plt.figure(figsize=(9,7))
+	ax = fig.add_subplot(111)
+	im = ax.imshow(Eca, cmap='hot', origin='lower', extent=[x.min(),x.max(),y.min(),y.max()], norm=LogNorm())
+	ax.set_xlim([x.min(),x.max()])
+	ax.set_ylim([y.min(),y.max()])
+	ax.set_xlabel('x [arcsec]')
+	ax.set_ylabel('y [arcsec]')
+	fig.colorbar(im, ax=ax, label='Amplitude')
+	plt.tight_layout()
+	plt.show()
+
+	plot_real_valued_results(np.sqrt(Eca), x, y, fig_size=(9,7))