CODEWITHSUNDEEP
pythonphysicspropagationfresnel
harshit tiwari
harshit tiwari

Reputation: 31

Optical field propagation with angular spectrum method, and verification

I have written a code for gaussian beam propagation in free-space by angular spectrum method. To verify whether the simulation is correct or not, I am propagating the gaussian beam in free-space, and verify the resulted field's beam width with the theoretical formulas. The following is my code, I don't understand what's wrong in it, the beam width is not matching with theoretical beam width, I would really appreciate any guidance. Thank you

Here is my code:

wvl   = 1550e-9             # wavelength
wz    = 0.05                # beam waist
gdim  = 1                   # spatial extent of the grid;
resol = 512                 # grid resolution
I     = 1                   # intensity amplitude
zr    = (np.pi * wz**2)/wvl # rayleigh range
k     = 2*np.pi/wvl

# Theoretical beam width at reciever (Lz km) # Gaussian 
wzt     = lambda Lz, wg, lmda: wg * np.sqrt(1 + ((Lz * lmda)/(np.pi * wg**2))**2)

wz_     = wzt(0, wz, wvl)
z       = zr

# define grid
dx      =  np.sqrt((wvl*z)/resol) #(wvl * np.sqrt(z**2 + gdim**2))/gdim
gdim    =  dx * resol
x, y    =  np.meshgrid(np.arange(-gdim/2, gdim/2, dx),
                     np.arange(-gdim/2, gdim/2, dx))

r    = np.sqrt(x**2 + y**2)

def gaussian(wz, r, I):
    """
    wz: Beam width at z=0;
    r : Radial coordinates
    I : Intensity distribution
    """
    Fin = I*np.exp(-2*(r/wz)**2)
    return Fin
gbeam = gaussian(wz, r, I)
def ft2(g, delta=None):
    """
    Computes the 2D Fourier transform of g with scaling.
    
    Parameters:
    g : numpy.ndarray
        Input 2D array to transform.
    delta : float
        Spacing in the spatial domain.
    
    Returns:
    numpy.ndarray
        The scaled 2D Fourier transform of g.
    """
    G = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(g)))
    return G

def ang_spec_prop(e, z, resol, dx, wvl):
    k = (2*np.pi/wvl)
    
    fxx = np.fft.fftshift(np.fft.fftfreq(resol, dx))
    fxx, fyy = np.meshgrid(fxx, fxx)
    
    alfa = k ** 2 - (4 * (np.pi ** 2)) * (fxx ** 2 + fyy ** 2)
    
    tmp = np.sqrt(np.abs(alfa))
    kz  = np.where(alfa >= 0, tmp, 1j*tmp)
     
    h1   = np.exp(z * kz * 1j)
    e_ = ft2(ft2(e)*h1)    
    return e_
prop_as = ang_spec_prop(gbeam, z, gbeam.shape[0], dx, wvl)

# Calculate the beam width
def anacal_gbw(ub, r):
    # Find the intensity where it falls to 1/e^2 of its maximum
    ub       = np.abs(ub)**2
    r_flat   = r.flatten()
    imax     = ub.max()/(np.exp(1)**2)
    in_flat  = ub.flatten()
    wg_calc  = r_flat[np.abs(in_flat-imax).argmin()]
    return wg_calc
ana_w = anacal_gbw(prop_as, r) 
the_w = wzt(z, wz, wvl)

With these parameters, the_w (theoretical beam width) = 0.0707; but ana_w (analytical beam width) = 0.0792

Upvotes: 0

Views: 33

Answers (1)

Neil Butcher
Neil Butcher

Reputation: 563

your gaussian function seems to have a factor of 2 in, try this instead:

def gaussian(wz, r, I):
    """
    wz: Beam width at z=0;
    r : Radial coordinates
    I : Intensity distribution
    """
    Fin = I*np.exp(-(r/wz)**2)
    return Fin

With this change the calculated value comes out as 0.07060

If you also increase the grid resolution to 1024 it becomes 0.07077

Are you sure that 2 supposed to be there?

Upvotes: 0

Related Questions