I want to simulate how a Gaussian beam would look like at the receiver plane when propagated through an atmospheric turbulence. For this I am using AO package. Since I didn't see any function for multipropagation i.e. choosing multiple phase screens between the distance and propagating partially. Taking reference from the book "Numerical Simulation of optical wave propagation with examples in Matlab", i modified the source code for angularSpectrum().
My code:
# -*- coding: utf-8 -*-
"""
@author:
"""
import numpy as np
import matplotlib.pyplot as plt
from aotools import phasescreen, opticalpropagation, atmos_conversions, fouriertransform
class beam:
def init(self, sigma, Eo, z, x, lambda_, r0, N, delta, L0, l0, num):
self.sigma = sigma
self.Eo = Eo
self.z = z
self.x = x
self.y = x
self.X,self.Y = np.meshgrid(self.x, self.y)
self.lambda_ = lambda_
self.r0 = r0
self.N = N
self.delta = delta
self.L0 = L0
self.l0 =l0
self.num = num
def electric_field(self): #generate a Gaussian beam
return self.Eo * np.exp(-0.5*(self.X **2 + self.Y**2)/self.sigma**2)
def PS(self): #generate phase screens
phase_arr = []
for i in range(self.num):
phase_arr.append(phasescreen.ft_sh_phase_screen(self.r0, self.N, self.delta, self.L0, self.l0, seed = None))
return np.array(phase_arr)
def angularSpectrum(self, inputSpacing, outputSpacing, t): #modified this from source code
"""
Propogates light complex amplitude using an angular spectrum algorithm
Parameters:
inputComplexAmp (ndarray): Complex array of input complex amplitude
wvl (float): Wavelength of light to propagate
inputSpacing (float): The spacing between points on the input array in metres
outputSpacing (float): The desired spacing between points on the output array in metres
z (float): Distance to propagate in metres
Returns:
ndarray: propagated complex amplitude
"""
wvl = self.lambda_
z = self.z
inputComplexAmp = self.electric_field()
# If propagation distance is 0, don't bother
# if self.z==0:
# return inputComplexAmp
mag = float(outputSpacing)/inputSpacing
N = self.N
k = 2*np.pi/wvl #optical wavevector
(x1,y1) = np.meshgrid(inputSpacing*np.arange(-N/2,N/2), #receiver plane meshgrid
inputSpacing*np.arange(-N/2,N/2))
r1sq = (x1**2 + y1**2) + 1e-10
Q1 = np.exp( 1j * k/2. * (1-mag)/delta * r1sq) #not changed
Uin = inputComplexAmp *Q1 * t[0]
w = N * 0.47
sg = np.exp(-(self.X**2 + self.Y**2)**8/(w**16))
for idx in range(self.num):
#Spatial Frequencies (of source plane)
df1 = 1. / (N*inputSpacing)
fX,fY = np.meshgrid(df1*np.arange(-N/2,N/2),
df1*np.arange(-N/2,N/2))
fsq = fX**2 + fY**2
#Scaling Param
Q2 = np.exp(-1j * np.pi**2 * 2 * z/mag/k*fsq)
Uin = sg * t[idx]*fouriertransform.ift2(
Q2 * fouriertransform.ft2(Uin/mag,inputSpacing), df1)
#Observation Plane Co-ords
x2,y2 = np.meshgrid( outputSpacing*np.arange(-N/2,N/2),
outputSpacing*np.arange(-N/2,N/2) )
r2sq = x2**2 + y2**2
#Quadratic phase factors
Q3 = np.exp(1j * k/2. * (mag-1)/(mag*z) * r2sq)
#Compute propagated field
outputComplexAmp = Q3 * Uin*Q1
return outputComplexAmp
def Fresnelpropagation(self): #optical propagation
phase_arr = self.PS()
Uout = self.angularSpectrum(self.delta, self.delta, phase_arr)
return np.array(Uout)
def avg_real(self): #take average of few iterations
av_m = 100
Uout = self.Fresnelpropagation()
Uav = np.zeros((Uout.shape))*1j
for i in range(av_m):
Uav += Uout
return Uav/av_m
def plot_(Uin,G):
# Assuming you have defined m, X, Y, and G appropriately for each propagating distance
m = G.shape[0]
fig, ax = plt.subplots(2, m+1) # Create n subplots in a single row
ax[0, 0].imshow(np.abs(Uin), cmap = 'viridis')
ax[1, 0].imshow(np.imag(Uin), cmap = 'viridis') # Display the absolute values of G[i] in each su
for i in range(1,m+1):
ax[0, i].imshow(np.abs(G[i-1]), cmap = 'viridis')
ax[1, i].imshow(np.imag(G[i-1]), cmap = 'viridis') # Display the absolute values of G[i] in each subplot
plt.show() # Display the plot
def func(cn2_,N): #try for a few cn2 values
count = 1
for cn2 in cn2_:
x = np.linspace(-D/2,D/2,N)
print(N)
r0 = atmos_conversions.cn2_to_r0(cn2, lambda_)
print(D/r0)
z = 1e2 #distance between two phase screens
#convert cn2 to ro
# sigma, Eo, z, x, lambda_, r0, N, delta, L0, l0, num
g = beam(sigma, Eo, z, x, lambda_, r0, N, delta, L0, l0, num)
Uin = g.electric_field()
Uav.append(g.avg_real(count))
count += 1
print(np.array(Uav).shape[0])
plot_(Uin, np.array(Uav))
plt.imshow(gaussian.electric_field())
D = 0.3 #receiver size
N = 100 #phase screen grid size
num = 10 #numner of phase screen
Eo = 1
delta = D/N #size of each pixel in a phase screen grid
L0 = 1 #the outer scale
l0 = 0.05 #the inner scale
lambda_ = 106410(-9) #1064nm
sigma = 0.01 #for beam width
#different values of cn2[![output][2]][2]
cn2_ = [110(-19), 1*10(-13), 110*(-10)]
Uav = []
func(cn2_,N)
I don't know if i did it correctly but the output I see is weird. I'd appreciate some help with understanding the output, if the code is correct. P.S the first plot is the source and the rest are on the basis of increasing atmospheric turbulence.
when running code for multiple times, i get different results. I do understand that the turbulence is random but shouldn't the output be somewhat similar based on the cn2 values? Output:
