Atmospherical propagation of a gaussian beam.

[Aziadocs]

Hi, my name is Fernando and I am from the University of Porto. My goal is to simulate the propagation of a gaussian beam through turbulent atmosphere, and maybe aoTools can do this, right? From what I can see, I must propagate using Fresnel Integral, for example for 2500m, then multiply by e^(j*kolmogorov phase screen), and then propagate 2500m again to complete 5000km. To add more turbulence I will need to rise the number of phase screens right?

I am having some troubles, even propagating a beam 1m...

1. First I created my gaussian:
   "Create gaussian profile" dim=0.3 #m hdim=dim/2 step=0.0001 lmbda0=1.550e-6 #m x = numpy.arange(-hdim,hdim, step) y = numpy.arange(-hdim,hdim, step) xx, yy = numpy.meshgrid(x, y, sparse=True) sigma=0.05 z = 1/(2*numpy.pi) * numpy.exp(-(xx**2+yy**2)/(2*sigma*sigma))

2. Then propagated 1 meter:
   hlinha=opticalpropagation.oneStepFresnel(z,lmbda0, step, 1)

and if I plot this I get something crazy like this

I don't know what I missed... must d1 value d1 be the step?
Other thing, why there are no cross terms in the Huygens-Fresnel-Kirchoff integral inside the oneStepFresnel function?
Why shouldn't be C = ft2(Uin *numpy.exp(-1j * k/(2*z) * (x1**2 + y1**2 - 2*x1*x2 - 2*y1*y2)), d1) instead?

Formula from Svelto :
$u(x, y, z)=\frac{j}{\lambda L} \iint_{s} u\left(x_{1}, y_{1}, z_{1}\right) \exp -j k\left[\frac{\left(x-x_{1}\right)^{2}+\left(y-y_{1}\right)^{2}}{2 L}\right] d x_{1} d y_{1}$

Thank you a lot ^_^

[matthewtownson]

To be honest I am not sure how the 1 step fresnel propagation works, it is something I will have to look into. I can reproduce the results you show, so I don't think you are doing anything wrong.

I usually use the angular spectrum propagation function, when I try your example with this the results are more like what I would expect. It doesn't answer your question about 1 step fresnel, but if the angular spectrum propagation is suitable for your use then I would recommend using that for now.

[Aziadocs]

Hi, Matthew.
Something is odd for sure... Probably it is the fft2, I will check it later!
Now I tried with angularSpectrum and it gave me a gaussian for 1 m propagation (cool). But the initial peak of the gaussian was 0.12, and with 1m propagation lowered to 0.020 (a factor of 6....). I would appreciate if you could help me on this one :) It would be awesome if I could get some cool results to put on my thesis!
wavelength = 1550e-9 propagation_distance=5 wavefront=z pxl_scale=step propagated_screen=opticalpropagation.angularSpectrum(wavefront,wavelength, pxl_scale, pxl_scale,propagation_distance) scintillation_pattern=numpy.abs(propagated_screen)**2 a = plt.axes(projection='3d') a.plot_surface(xx, yy, scintillation_pattern, cmap="rainbow");

[matthewtownson]

I am not sure about the peak of the Gaussian, I would expect it to reduce as the beam propagates but don't know how much off the top of my head. It is something that might be worth checking. If you sum up the intensity of the beam before and after propagation they are the same, so the light is not being lost in the propagation. I would guess it is an effect to do with the size of the aperture you are propagating and the beam spreading quickly

[matthewtownson]

Hi @Aziadocs have you had a chance to check this is behaving as you would expect?

[Aziadocs]

Hello, Matthew.
Yes, I did. The method involves a really large matrix for my parameters so I twisted it a little bit. I was getting boundary artifacts and stuff like that.
Now I am getting into my university cluster so I can get the memory to do those calculations.
thank you for checking!
Fernando Maia

[KlenM]

Hi, @Aziadocs, i would recommend using a GPU rather than a CPU to simulate propagation with large matrices, because it can speed up the simulation time by several orders. (I just ran a simulation for a random channel and got a 40x speedup).
This is easy to do by simply replaced numpy (np) module in the code with cupy. (Or another library, jax.numpy)

If you don't have a GPU for this, Google Colab provides free GPU usage for 5-12 hours a day.

If you have really large matrices (>8k*8k), Google Colab provides a free-use of TPU, which can speed up simulation by several times compared to GPU. But I haven't tested this.

[matthewtownson]

Great to hear you're making progress.

And @NYurin is correct that these methods are not the most efficient if you are doing many runs at large scale, and they can be accelerated relatively easily. Accelerating parts of aotools is something we are considering, but for now we are erring on the side of caution on adding new dependencies.

I think this issue can be closed now unless you have anything else related?

[tozbilgin]

To be honest I am not sure how the 1 step fresnel propagation works, it is something I will have to look into. I can reproduce the results you show, so I don't think you are doing anything wrong.

I usually use the angular spectrum propagation function, when I try your example with this the results are more like what I would expect. It doesn't answer your question about 1 step fresnel, but if the angular spectrum propagation is suitable for your use then I would recommend using that for now.

I have been working on Gaussian beam propagation based on this post. However I have a few questions.

First, how should I create a Gaussian wavefront, I use the complex beam parameter q(z) (Alternatively R(z) and w(z)).
Second, how to propagate the wavefront using opticalpropagation function of AOtools. I create Kolmogorov or von Karman phasescreens to include turbulence. On the other hand I have the Gaussian wavefront. I don't have any idea to combine them, do I just multiply them. I am confused :(

I am sharing the code I am working on.

import aotools
import numpy as np
from aotools import opticalpropagation
from aotools.turbulence import infinitephasescreen, phasescreen
import matplotlib.pylab as plt
from tqdm import tqdm, trange, tqdm_notebook


D = 0.3 #1. total size of the grid [m] 
stencil_length_factor = 32
nx_size = 64  #number of grid points per side
pxl_scale = D/nx_size # % grid spacing [m] wvl*z/W

x = np.arange(-D/2,D/2, pxl_scale) 
y = np.arange(-D/2,D/2, pxl_scale) 
xx, yy = np.meshgrid(x, y, sparse=True)

wavelength = 810e-9 # (m) 
k = 2*np.pi/wavelength #Wave number

cn2 = 1e-13
r0 = aotools.turbulence.atmos_conversions.cn2_to_r0(cn2=cn2, lamda=wavelength)

l0=0.001 # inner scale (m)
L0 = 100 / (nx_size/2) # outer scale (m)
j = 0+1j

L = 5e3 #Distance (m)

w0 = 2.5e-2 # inital beam radius (m)
teta0 = 10e-6 #beam divergence (rad)
wL = w0/2 + L*teta0 #final beam width at receiver (without turbulence) (m)
zR = np.pi*w0**2/wavelength #Rayleigh range

P = 0.1 #Total power (W)
IL = 2*P / (np.pi*wL)**2 #I(r,z) = 2*P / (pi*w(z)^2) * exp(-2*r^2 / w(z)^2)


def qz(z): #complex beam parameter
    return z + j*zR
gaussian_wavefront = 1/qz(L) * np.exp( -j*k*(xx**2+yy**2)/(2*qz(L)) ) # *np.exp(j*k*L)


#Theoretic calculation for comparison
sR2 = 1.23*cn2*k**(7/6)*L**(11/6) #Rytov variance
wRxL = wL *np.sqrt(1 + 1.33*sR2 *(2*L/ (k*wL**2))**(5/6)) #Atmospheric turbulence added

print("w0= ", w0, ", wRx_no_turb=",wL, ", wRx_with_turb= ", wRxL)





# Generate a large set of phase screens
N_scrns = 2

#wavefronts = np.zeros((N_scrns, nx_size, nx_size))
Is = np.zeros((N_scrns, nx_size, nx_size))
powerspec = np.zeros_like(Is)

for i in tqdm(range(N_scrns)):
    phase_screen = phasescreen.ft_phase_screen(r0, nx_size, pxl_scale, L0, l0, FFT=None, seed=None) #von Karman
    #phase_screen = infinitephasescreen.PhaseScreenVonKarman(nx_size, pxl_scale, r0, L0).scrn
    #phase_screen = infinitephasescreen.PhaseScreenKolmogorov(nx_size, pxl_scale, r0, L0, stencil_length_factor=stencil_length_factor).scrn    
    
    
    wavefront = gaussian_wavefront * np.exp(j*phase_screen)
    
    m=1
    propagated_wavefront = aotools.opticalpropagation.angularSpectrum(wavefront, wavelength, m*pxl_scale, m*pxl_scale, L)
    
    Is[i] = IL*np.abs(propagated_wavefront)**2

    
I=Is.mean(0)   

ax = plt.axes(projection='3d') 
ax.plot_surface(xx, yy, I, cmap="rainbow");

import matplotlib.colors as colors
X, Y = np.meshgrid(x, y)
fig,ax2=plt.subplots(1,1)
cp = ax2.contourf(X, Y, I, cmap="rainbow")
fig.colorbar(cp)
plt.show

[Aziadocs]

Yes, that's exactly how I did it! Just multiply them! Like eq (2) in DOI:10.1117/12.2033442
I manage to modify the angular spectrum method and I am getting some interesting results for 7km propagation. I was interested in measuring the beam deviation. Now I am writing my thesis oof

[tozbilgin]

Yes, that's exactly how I did it! Just multiply them! Like eq (2) in DOI:10.1117/12.2033442
I manage to modify the angular spectrum method and I am getting some interesting results for 7km propagation. I was interested in measuring the beam deviation. Now I am writing my thesis oof

Good luck with your thesis Aziadocs!

Yes I am multiplying the two by wavefront = gaussian_wavefront * np.exp(j*phase_screen) and propagate this wavefront with propagated_wavefront = aotools.opticalpropagation.angularSpectrum(wavefront, wavelength, pxl_scale, m*pxl_scale, L), however I cannot see any difference unless I extremely increase the distance value.

For example for a 7km system with wavelength=810nm, Cn2=1e-13 and initial beam radius w0=2.5cm the final beam radius should be (theoretically) 8.25 cm without turbulence and 67 cm with turbulence (all codes are below). However I don't see any difference.

I could not find the problem. Can someone help me or can you share some working code of yours?

  import aotools
  import numpy as np
  from aotools import opticalpropagation
  from aotools.turbulence import infinitephasescreen, phasescreen
  import matplotlib.pylab as plt
  from tqdm import tqdm, trange, tqdm_notebook
  
  
  D = .2 #1. total size of the grid [m] 
  stencil_length_factor = 32
  nx_size = 64  #number of grid points per side
  pxl_scale = D/nx_size # % grid spacing [m] wvl*z/W
  
  x = np.arange(-D/2,D/2, pxl_scale) 
  y = np.arange(-D/2,D/2, pxl_scale) 
  xx, yy = np.meshgrid(x, y, sparse=True)
  
  wavelength = 810e-9 # (m) 
  k = 2*np.pi/wavelength #Wave number
  
  cn2 = 1e-13
  r0 = aotools.turbulence.atmos_conversions.cn2_to_r0(cn2=cn2, lamda=wavelength)
  
  l0=0.001 # inner scale (m)
  L0 = 100 / (nx_size/2) # outer scale (m)
  j = 0+1j
  
  L = 7e3 #Distance (m)
  
  w0 = 2.5e-2 # inital beam radius (m)
  teta0 = 10e-6 #beam divergence (rad)
  wL = w0/2 + L*teta0 #final beam width at receiver (without turbulence) (m)
  zR = np.pi*w0**2/wavelength #Rayleigh range
  
  P = 0.1 #Total power (W)
  IL = 2*P / (np.pi*wL)**2 #I(r,z) = 2*P / (pi*w(z)^2) * exp(-2*r^2 / w(z)^2)
  
  
  def qz(z): #complex beam parameter
      return z + j*zR
  gaussian_wavefront = 1/qz(L) * np.exp( -j*k*(xx**2+yy**2)/(2*qz(L)) ) # *np.exp(j*k*L)
  
  
  
  
  
  #Theoretic calculation for comparison
  sR2 = 1.23*cn2*k**(7/6)*L**(11/6) #Rytov variance
  wRxL = wL *np.sqrt(1 + 1.33*sR2 *(2*L/ (k*wL**2))**(5/6)) #Atmospheric turbulence added
  
  print("w0= ", w0, ", wRx_no_turb=",wL, ", wRx_with_turb= ", wRxL)
  
  
  
  
  
  # Generate a large set of phase screens
  N_scrns = 2
  
  #wavefronts = np.zeros((N_scrns, nx_size, nx_size))
  Is = np.zeros((N_scrns, nx_size, nx_size))
  powerspec = np.zeros_like(Is)
  
  for i in tqdm(range(N_scrns)):
      phase_screen = phasescreen.ft_phase_screen(r0, nx_size, pxl_scale, L0, l0, FFT=None, seed=None) #von Karman
      #phase_screen = infinitephasescreen.PhaseScreenVonKarman(nx_size, pxl_scale, r0, L0).scrn
      #phase_screen = infinitephasescreen.PhaseScreenKolmogorov(nx_size, pxl_scale, r0, L0, stencil_length_factor=stencil_length_factor).scrn    
      
  
      wavefront = gaussian_wavefront * np.exp(j*phase_screen)
      
      m=1
      propagated_wavefront = aotools.opticalpropagation.angularSpectrum(wavefront, wavelength, pxl_scale, m*pxl_scale, L)
      
      Is[i] = IL*np.abs(propagated_wavefront)**2
  
      
  I=Is.mean(0)   
  
  ax = plt.axes(projection='3d') 
  ax.plot_surface(m*xx, m*yy, I, cmap="rainbow");
  
  import matplotlib.colors as colors
  X, Y = np.meshgrid(m*x, m*y)
  fig,ax2=plt.subplots(1,1)
  cp = ax2.contourf(X, Y, I, cmap="rainbow")
  fig.colorbar(cp)
  plt.show

[matthewtownson]

sorry for the long silence. I am not sure about your initial Gaussian beam, it looks like you put a Gaussian shape to the complex (phase) part, I think this should be applied to the amplitude and the phase part is flat.

[matthewtownson]

Has this been resolved? Either way, I think this is a discussion on using aotools rather than an issue within aotools itself. So I am moving it from here to discussions.