From 914910e684a7436e44d075784c12823d58b7bee2 Mon Sep 17 00:00:00 2001 From: yosata Date: Wed, 3 Feb 2021 10:16:54 +0100 Subject: [PATCH] Added Atmospheric Seeing --- fourierOptics/__init__.py | 3 +- fourierOptics/atmosphere.py | 135 ++++++++++++++ fourierOptics/fourier.py | 15 +- fourierOptics/misc.py | 82 ++++++++- fourierOptics/optics.py | 13 +- scripts/SPIE/generate_ADC_footprint.py | 4 +- .../gratama_seeing_limited_dispersion_test.py | 172 ++++++++++++++++++ 7 files changed, 404 insertions(+), 20 deletions(-) create mode 100644 fourierOptics/atmosphere.py create mode 100644 scripts/gratama_seeing_limited_dispersion_test.py diff --git a/fourierOptics/__init__.py b/fourierOptics/__init__.py index 1a9800b..14e3a7d 100644 --- a/fourierOptics/__init__.py +++ b/fourierOptics/__init__.py @@ -7,4 +7,5 @@ from .profiles import * from .psfFitter import * from .wavefrontErrors import * from .misc import * -from .micado import * \ No newline at end of file +from .micado import * +from .atmosphere import * \ No newline at end of file diff --git a/fourierOptics/atmosphere.py b/fourierOptics/atmosphere.py new file mode 100644 index 0000000..ba74c9f --- /dev/null +++ b/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 + + diff --git a/fourierOptics/fourier.py b/fourierOptics/fourier.py index 01f9b21..3c21a64 100644 --- a/fourierOptics/fourier.py +++ b/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': diff --git a/fourierOptics/misc.py b/fourierOptics/misc.py index dbb8989..8c211e8 100644 --- a/fourierOptics/misc.py +++ b/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) diff --git a/fourierOptics/optics.py b/fourierOptics/optics.py index fca9642..10bcfea 100644 --- a/fourierOptics/optics.py +++ b/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) diff --git a/scripts/SPIE/generate_ADC_footprint.py b/scripts/SPIE/generate_ADC_footprint.py index b262fa5..47efc91 100644 --- a/scripts/SPIE/generate_ADC_footprint.py +++ b/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() diff --git a/scripts/gratama_seeing_limited_dispersion_test.py b/scripts/gratama_seeing_limited_dispersion_test.py new file mode 100644 index 0000000..81a7cbe --- /dev/null +++ b/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)) \ No newline at end of file