Tutorial: LightPipes for Python¶

Developers of LightPipes for Python:
Gleb Vdovin, Fred van Goor, Guyskk, Leonard Doyle, Flexible Optical B.V. (OKO Tech)

Package to model wave optics phenomena and beam propagation based on Scalar theory of Diffraction (Fresnel-Kirchoff)
Versions 2.0.+, pure Python implementation using numpy, scipy and pyFFTW packages
Works with Python 3.+ versions
Supported on Windows, Linux and Mac

Install as: pip install LightPipes
Documentation

Features: Propagation through lenses with coordinate transforms, simulation of any combination of Zernike aberrations, mode analysis in laser resonators, interferometers, inverse problems, waveguides and propagation in media with non-uniform distribution of refraction index

from LightPipes import *
import numpy as np
import matplotlib
%matplotlib inline
import pylab as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import gridspec
import time

print(LPversion)
LPtest()

2.0.5
Test OK

Start calculations¶

`Begin`¶

Let's begin the journey with our first LightPipes function!

All LightPipes calculations start with this built-in function which defines the field on a square grid.

Field = Begin(size,wavelength,N)

Recall that 'Shift+Tab' inside a function shows the arguments it takes

Seek more help with the 'help' command as help(Begin)

size= 10.0*mm #size of the square grid
wavelength= 500*nm
N= 3 #grid dimension

Field= Begin(size,wavelength,N) #A square field at this point

type(Field) #It's a LightPipes field object

LightPipes.field.Field

print(type(Field.field), Field.field.dtype, Field.field.shape)

<class 'numpy.ndarray'> complex128 (3, 3)

Field.field

array([[1.+0.j, 1.+0.j, 1.+0.j],
       [1.+0.j, 1.+0.j, 1.+0.j],
       [1.+0.j, 1.+0.j, 1.+0.j]])

So, we've defined a uniform field with intensity 1 and phase 0 everywhere

Apertures¶

Circular aperture, CircAperture
Rectangular aperture, RectAperture
Gaussian aperture, GaussAperture

Screens¶

Circular screen, CircScreen
Rectangular screen, RectScreen
Gaussian screen, GaussScreen

Screens are the inversions of their respective apertures

Screens and apertures:¶

size= 10.0*mm 
wavelength= 500*nm
N= 512 
Field= Begin(size, wavelength, N)

R= 3*mm #Aperture/Screen radius
x_shift= 2*mm; y_shift= 2*mm #Shift of the aperture/screen

Caperture= CircAperture(Field,R,x_shift,y_shift)
Iaperture= Intensity(Caperture)

Cscreen= CircScreen(Field,R,x_shift,y_shift)
Iscreen= Intensity(Cscreen)

fig= plt.figure(figsize=(10,6))
ax1= fig.add_subplot(121)
ax2= fig.add_subplot(122)
ax1.imshow(Iaperture); ax1.axis('off'); ax1.set_title('Circular aperture')
ax2.imshow(Iscreen); ax2.axis('off'); ax2.set_title('Circular screen')
plt.show()

Combination of screens and apertures:¶

size= 25.0*mm
wavelength= 500*nm
N= 512
Field = Begin(size, wavelength, N) 

R1= 5*mm; x_shift= 0*mm; y_shift = 0*mm; T= 0.8
F1= GaussAperture(R1,x_shift,y_shift,T,Field) #A 5mm Gaussian aperture

R2= 2*mm
F2= CircScreen(R2,x_shift,y_shift,F1) #A circular screen

R3= 4*mm
F3= CircAperture(R3,x_shift,y_shift,F2) #A circular screen

#Field intensities
I1= Intensity(F1)  
I2= Intensity(F2)
I3= Intensity(F3)

fig= plt.figure(figsize=(16,9))
ax1= fig.add_subplot(131)
ax2= fig.add_subplot(132)
ax3= fig.add_subplot(133)
ax1.imshow(I1,cmap='hot'); ax1.axis('off'); ax1.set_title('Gaussian aperture')
ax2.imshow(I2,cmap='hot'); ax2.axis('off'); ax2.set_title('+ Circular screen')
ax3.imshow(I3,cmap='hot'); ax3.axis('off'); ax3.set_title('+ Circular aperture')
plt.show()

t0 = time.time()
X= range(N); Y= range(N)
X, Y= np.meshgrid(X,Y) 
fig= plt.figure(figsize=(18,12))
ax= fig.gca(projection='3d')
ax.plot_surface(X, Y, I3,
                rstride=1,cstride=1,cmap='rainbow',linewidth = 0.0)
ax.set_zlabel('Intensity [a.u.]')
ax.set_title('Field intensity after the Gaussian aperture + circular screen + circular aperture')
plt.show()
print('Took', round(time.time()-t0,2), 'sec')

Took 36.29 sec

`Interpol`¶

Interpolates the field to a new grid size, grid dimension.

So with Interpol, we can reduce the grid dimension of the field to speed things up, especially when dealing with 3d surface plots.

# Field intensity after the Gaussian aperture (using Interpol)
NewGridDimension = int(N/4)
F1_interp = Interpol(size,NewGridDimension,0,0,0,1,F1)
I1_interp = Intensity(F1_interp)

t0 = time.time()
X= range(NewGridDimension)
Y= range(NewGridDimension)
X, Y= np.meshgrid(X,Y) 

fig = plt.figure(figsize=(18,6)) 
gs = gridspec.GridSpec(1, 2, width_ratios=[1,2]) 
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1], projection='3d')

### first subplot: cross section of the field intensity
x= []
for i in range(NewGridDimension):
    x.append((-size/2+i*size/NewGridDimension)/mm)
ax1.plot(x,I1_interp[int(NewGridDimension/2)])
ax1.set_xlabel('x [mm]')
ax1.set_ylabel('Intensity [a.u.]')

### second subplot: 3d surface plot of the field intensity
ax2.plot_surface(X,Y,I1_interp,rstride=1,cstride=1,cmap='rainbow',linewidth=0.0)
ax2.set_zlabel('Intensity [a.u.]')
ax2.set_title('Field intensity: Gaussian aperture')

plt.tight_layout()
plt.show()
print('Took', round(time.time()-t0,2), 'sec')

Took 2.62 sec

Free-space propagation¶

The Fresnel-Kirchoff formulation is a scalar theory of diffraction and enables us to predict diffraction effects under these conditions:

aperture is large compared to the wavelength of light, $a>>\lambda$
field is calculated at a distance much larger than the wavelength of light, $z>>\lambda$
$a$: characterstic size (radius) of the diffracting aperture
$\lambda$: wavelength
$z$: propagation distance

Fresnel and Fraunhofer diffraction¶

In Fresnel-Kirchhoff, the expanding wavefronts propagating away from the aperture are spherical

Fresnel and Fraunhofer diffractions are approximations of the Fresnel-Kirchhoff formalism

After a certain distance from the aperture and in the proximity of the optic axis, the wavefronts can be approximated as parabolic: Fresnel or near-field approximation

After a sufficiently large distance and in the proximity of the optic axis, the wavefronts can be further approximated as planar wavefronts: Fraunhofer or far-field approximation

Note that Fraunhofer is a special case of Fresnel, so Frensel approximation applies in the Fraunhofer regime too

Near and Far field:¶

Fresnel number, $F =a^2/\lambda z$ coarsely determines the near/far diffraction regimes:

When $F \geq 1 \Longrightarrow z \leq a^2/\lambda$ , this is referred to as near-field
Spherical wavefronts can be approximated as parabolic at this distance
Fresnel approximation applies

When $F << 1 \Longrightarrow z >> a^2/\lambda$, the is referred to as far-field
Spherical wavefronts can be approximated as planar at this distance
Fraunhofer approximations applies (Fresnel works too)

For example: With a 1 inch aperture and a 500nm beam, $a^2/\lambda$= 0.2 miles. So, the observation distance $z$ must be larger than 2 miles (an order of magnitude larger) for the Fraunhofer approximation to be valid!

Validity ranges:¶

Rayleigh-Sommerfeld: $a >> \lambda$ and $z \gtrsim \lambda$

Fresnel-Kirchoff: $a >> \lambda$ and $z >> \lambda$

Fresnel: $z >> a >> \lambda$

Fraunhofer: $z >> a^2/\lambda$

For Fresnel and Fraunhofer: observation plane close to the optic axis such that for any observation point on the plane, $r \approx z$

Typically a factor or 10x or more when $>>$

Fourier transform relations:¶

Fresnel diffraction: Field in the image plane is a 2D Fourier transform of the field at the aperture multiplied by a quadratic phase factor

Fraunhofer diffraction: Field in the image plane is a 2D Fourier transform of the field at the aperture

Lens: a Fourier transform operator¶

Field in the image plane is a 2D Fourier transform of the field at the pupil or back focal plane of the lens

It brings to focus the diffraction pattern that would otherwise be observed at infinity (Fraunhofer diffraction pattern)

Four methods available to model free-space propagation of light:

Forvard Propagate the field using FFT
Fresnel Propagate the field using convolution followed by FFT
Forward Propagate the field using direct integration
Steps Propagate the field using a Finite difference method

`Forvard(Fin,z)`¶

Propagates the field using FFT to solve Fresnel-Kirchoff integral
Simplest and fastest propagation command
Models light propagation inside a square waveguide with reflecting walls at the grid edges
To approximate a free space propagation, need to set the intensity near the grid edges to be negligible

def pltmethod(gridsize, lamda, gridN, aperture, Intensity, zval, method, title):
    "Generate propagation plot using a single method"
    #--------------------
    # Check diffraction regime    
    FresnelNumber= round(aperture**2/(lamda*zval),4)
    print('Fresnel number (F):', FresnelNumber)
    if FresnelNumber < 0.1:
        title= title + 'Far-field diffraction (F= ' + str(FresnelNumber) + ') modeled using ' + method
    elif FresnelNumber >= 1:
        title= title + 'Near-field diffraction (F= ' + str(FresnelNumber) + ') modeled using ' + method
    else:
        title= title + 'Diffraction (F= ' + str(FresnelNumber) + ') modeled using ' + method
    #--------------------
    #Line profiles
    yy= Intensity[int(gridN/2)]
    xx=np.arange(gridN)
    xx=(xx/gridN-1/2)*gridsize/mm

    #--------------------
    fig= plt.figure(figsize=(18,8))
    fig.suptitle(title, fontsize=18)
    ax1= fig.add_subplot(121)
    ax2= fig.add_subplot(122)
    
    ax1.imshow(Intensity,cmap='hot')
    ax1.axis('off')
    ax1.set_title('After propagating '  + str(zval/m) + ' m in z using ' + method)

    ax2.plot(xx,yy)
    ax2.set_xlabel('x [mm]')
    ax2.set_ylabel('Intensity [a.u.]')
    
    plt.tight_layout(pad=7)
    plt.show()

Propagation through a circular aperture:¶

size= 720.0*mm
wavelength= 0.5*um
N= 2304 
Field=Begin(size,wavelength,N)

title = 'Circular Aperture: '
R= 1.5*mm; x_shift= 0*mm; y_shift = 0*mm
Field_z0= CircAperture(Field,R,x_shift,y_shift)
I_z0= Intensity(Field_z0)

#Propagation distance z= [0.045m; 4.5m; 450m] ==> F= [100, 1, 0.1]
z= 450*m
   
#Propagate the field
t0= time.time()
F= Forvard(Field_z0,z)
I= Intensity(1,F)
print('Forvard: Took ', round(time.time()-t0,4), 'sec')

Forvard: Took  0.7611 sec

print('Airy disk first minima (in Fraunhofer regime) should be at',round((1.22*z*wavelength)/(2*R)/mm, 2),'mm')
pltmethod(size, wavelength, N, R, I, z, 'Forvard', title)

Airy disk first minima (in Fraunhofer regime) should be at 91.5 mm
Fresnel number (F): 0.01

`Fresnel(Fin,z)`¶

Propagate the field using Convolution
2-5x slower compared to Forvard and consumes more memory
Frensel diffraction integral is converted to convolution then FFT is performed
Unlike Forvard, no need to have intensities vanish at the edges, so can use a smaller grid size
Fresnel results invalid if $z$ is comparable or less than the aperture at which the field is diffracted. Use Forvard or Steps that solve the full Fresnel-Kirchoff equation in these cases

size= 480.0*mm
wavelength= 0.5*um
N= 1536 
Field=Begin(size,wavelength,N)

title = 'Circular Aperture: '
R= 1.5*mm; x_shift= 0*mm; y_shift = 0*mm
Field_z0= CircAperture(Field,R,x_shift,y_shift)
I_z0= Intensity(Field_z0)

#Propagation distance z= [0.045m; 4.5m; 450m] ==> F= [100, 1, 0.1]
z= 450*m
   
#Propagate the field
t0= time.time()
F1= Fresnel(Field_z0,z)
I1= Intensity(1,F1) 
print('Fresnel: Took ', round(time.time()-t0,4), 'sec')

Fresnel: Took  1.6018 sec

print('Airy disk first minima (in Fraunhofer regime) should be at',round((1.22*z*wavelength)/(2*R)/mm, 2),'mm')
pltmethod(size, wavelength, N, R, I1, z, 'Fresnel', title)

Airy disk first minima (in Fraunhofer regime) should be at 91.5 mm
Fresnel number (F): 0.01

`Forward(Fin,z,sizenew,Nnew)`¶

Propagates the field using direct integration of the Fresnel-Kirchoff integral
Number of operations is proportional to $N^4$, where $N$ is the grid sampling
To reduce the overall time, size of the grid can be adjusted to match the cross section of field distribution
Forward also allows arbitrary sampling and size of square grid at the input and output planes

size= 2.0*mm # a small grid size large enough to contain the aperture below
wavelength= 0.5*um
N= 32
Field=Begin(size,wavelength,N)

title = 'Circular Aperture: '
R= 1.5*mm; x_shift= 0*mm; y_shift = 0*mm
Field_z0= CircAperture(Field,R,x_shift,y_shift)
I_z0= Intensity(Field_z0)

#Propagation distance z= [0.045m; 4.5m; 450m] ==> F= [100, 1, 0.1]
z= 450*m

newsize= 480.0*mm
newN= 256  

#Propagate the field
t0 = time.time()
F2= Forward(Field_z0,z,newsize,newN) 
I2= Intensity(1,F2)
print('Forward: Took ', round(time.time()-t0,4), 'sec')

Forward: Took  20.1499 sec

print('Airy disk first minima (in Fraunhofer regime) should be at',round((1.22*z*wavelength)/(2*R)/mm, 2),'mm')
pltmethod(newsize, wavelength, newN, R, I2, z, 'Forward', title)

Airy disk first minima (in Fraunhofer regime) should be at 91.5 mm
Fresnel number (F): 0.01

Propagation through a Square aperture:¶

size= 320.0*mm 
wavelength= 0.5*um
N= 1024
Field=Begin(size,wavelength,N)

title = 'Square Aperture: '
sx= 3*mm; sy= 3*mm;
x_shift= 0*mm; y_shift = 0*mm
Field_z0= RectAperture(Field,sx,sy,x_shift,y_shift)
I_z0= Intensity(Field_z0)

#Propagation distance z= [0.045m; 4.5m; 450m] ==> F= [100, 1, 0.1]
z= 450*m 

#Propagate the field
t0 = time.time()
F3=Forvard(Field_z0,z)
I3= Intensity(1,F3)
print('Forvard: Took ', round(time.time()-t0,4), 'sec')

t0 = time.time()
F4= Fresnel(Field_z0,z)
I4= Intensity(1,F4)
print('Fresnel: Took ', round(time.time()-t0,4), 'sec')

Forvard: Took  0.1639 sec
Fresnel: Took  0.6053 sec

print('First minima (in Fraunhofer regime) should be at',round((z*wavelength)/sx/mm, 2),'mm')
#pltmethod(size, wavelength, N, sx/2, I3, z, 'Forvard', title)
pltmethod(size, wavelength, N, sx/2, I4, z, 'Fresnel', title)

First minima (in Fraunhofer regime) should be at 75.0 mm
Fresnel number (F): 0.01

Propagation through a Gaussian aperture:¶

size= 320.0*mm 
wavelength= 0.5*um
N= 1024
Field= Begin(size,wavelength,N)

title = 'Gaussian Aperture: '
R= 1.5*mm; x_shift= 0*mm; y_shift = 0*mm; T= 0.8
Field_z0= GaussAperture(R,x_shift,y_shift,T,Field) 
I_z0= Intensity(Field_z0)

#Propagation distance z= [0.045m; 4.5m; 450m] ==> F= [100, 1, 0.1]
z= 450*m

#Propagate the field
t0 = time.time()
F6= Fresnel(Field_z0,z)
I6= Intensity(1,F6)
print('Fresnel: Took ', round(time.time()-t0,4), 'sec')

Fresnel: Took  0.5901 sec

pltmethod(size, wavelength, N, R, I6, z, 'Fresnel', title)

Fresnel number (F): 0.01

A Gaussian beam remains Gaussian at every point as it propagates through an optical system.

Diffraction from a Circular Screen: Poisson/Fresnel/Arago spot ¶

A fascinating story on diffraction that played out in the French Academy in 1818, in which a well-crafted counterargument ironically became a pillar of the opposiing view

size=10.0*mm
wavelength=0.5*um
N=1024
F=Begin(size,wavelength,N)
F=GaussBeam(F,1*mm)

title = 'Circular Screen: '
R=1*mm
F=CircScreen(F,R)

z=200*mm
Farago= Fresnel(F,z)
Iarago= Intensity(1,Farago)

pltmethod(size, wavelength, N, R, Iarago, z, 'Fresnel', title)

Fresnel number (F): 10.0

size=200.0*mm
wavelength=0.5*um
N=2000
F=Begin(size,wavelength,N)
F=GaussBeam(F,5*mm)

title = 'Circular Screen: '
R=50*mm
F=CircScreen(F,R)

z= 1600*m
Farago= Fresnel(F,z)
Iarago= Intensity(1,Farago)

pltmethod(size, wavelength, N, R, Iarago, z, 'Fresnel', title)

Fresnel number (F): 3.125

Note that Fresnel-Kirchhoff’s does not correctly match the measured intensity of Poisson spot.

The Rayleigh-Sommerfeld formulation is another scalar theory of diffraction and applies when $a >>\lambda$. It makes a slightly different assumption about the field and its derivative at the boundaries of the aperture. As such, it overcomes the mathematical inconsistencies encountered by Fresnel-Kirchoff formulation and predicts the correct intensity of the Poisson's spot.

The Rayleigh-Sommerfeld formulation also reproduces the diffracted field right behind the aperture ($z \gtrsim \lambda$), which Fresnel-Kirchoff fails to do. For this reason, the Rayleigh-Sommerfeld formulation is considered to be a "full" scalar solution.

Experimentally though, the Fresnel-Kirchoff formulation gives more accurate results when the observation plane is many wavelengths away form the diffracting aperture. Also, the Rayleigh–Sommerfeld theory is limited to a plane surface, which limits its utility since we usually deal with curved surfaces in optics.

Lastly, neither Fresnel-Kirchoff nor Rayleigh-Sommerfield take into account the coupled electric and magnetic vector fields. The vector nature cannot be ignored, for instance when $a \sim \lambda$ or $z < \lambda$

Using `Forvard`¶

Forvard is a numerical solution to Fresnel-Kirchoff diffraction integral

With Forvard, intensity at the edges of the grid must be negligible. Choose grid size >> $a$ i.e., a large grid with lots of zeros

Forvard also takes a negative $z$ values, in which case it performs “back propagation” i.e., it will reconstruct the initial field from the one diffracted

Using `Fresnel`¶

Fresnel is a numerical solution to Fresnel diffraction integral

Fresnel only takes positive $z$ values and doesn't require a large grid size

Use Fresnel whenever possible to model near and far-field diffraction as long as $z>>a$

When $z \sim a$, use Forvard instead of Fresnel

Using `Forward`¶

Forward is quite slow as computation time scales as $N^4$
Forward only takes positive $z$ values
Ideal when different grid sizes and grid elements are needed in the object and image planes

`Steps(Fin,dz,nstep,refr)`¶

Finite difference method: Propagate a distance nstep x dz in an absorptive medium with a complex refractive index stored as the square array refr

Extremely slow and results can be noisy if step sizes are large
Use when dealing with absorptive medium (i.e., refractive index with complex components, eg., waveguides)
Steps takes both positive and negative $z$ values

Propagation through a thin lens¶

size= 8*mm    
wavelength0= 500*nm
wavelength1= 800*nm
N= 512  

F0= Begin(size,wavelength0,N)
F0= GaussBeam(F0,w0= 2*mm) # TEM_0,0 Hermite Gauss beam

F1= Begin(size,wavelength1,N)
F1= GaussBeam(F1,w0= 2*mm) # TEM_0,0 Hermite Gauss beam

#Add a lens
f= 100*mm
F0= Lens(F0,f)
F1= Lens(F1,f)

F= F0

#Propagate to various z-positions past the lens using Fresnel
z1= 0.33*f; Fa= Fresnel(F,z1)
z2= 0.67*f; Fb= Fresnel(F,z2)
z3= f;      Fc= Fresnel(F,z3)
z4= 1.33*f; Fd= Fresnel(F,z4)
z5= 1.67*f; Fe= Fresnel(F,z5)

plt.imshow(Intensity(1,F),cmap='rainbow')
plt.title('TEM_0,0 Hermite Gauss beam')

fig= plt.figure(figsize=(18,4))
fig.suptitle('Propagation through a f=' + str(f/mm) + 'mm lens; for ' 
             + str(round(wavelength0/nm)) + 'nm', fontsize=18)
ax1= fig.add_subplot(151); ax1.set_title(str(z1/mm) + ' mm'); ax1.imshow(Intensity(1,Fa),cmap='rainbow')
ax2= fig.add_subplot(152); ax2.set_title(str(z2/mm) + ' mm'); ax2.imshow(Intensity(1,Fb),cmap='rainbow')
ax3= fig.add_subplot(153); ax3.set_title(str(z3/mm) + ' mm'); ax3.imshow(Intensity(1,Fc),cmap='rainbow')
ax4= fig.add_subplot(154); ax4.set_title(str(z4/mm) + ' mm'); ax4.imshow(Intensity(1,Fd),cmap='rainbow')
ax5= fig.add_subplot(155); ax5.set_title(str(z5/mm) + ' mm'); ax5.imshow(Intensity(1,Fe),cmap='rainbow')
plt.tight_layout(pad=5)
plt.show()

F= F1

#Propagate to various z-positions past the lens using Fresnel
z1= 0.33*f; Fa= Fresnel(F,z1)
z2= 0.67*f; Fb= Fresnel(F,z2)
z3= f;      Fc= Fresnel(F,z3)
z4= 1.33*f; Fd= Fresnel(F,z4)
z5= 1.67*f; Fe= Fresnel(F,z5)

fig= plt.figure(figsize=(18,4))
fig.suptitle('Propagation through a f=' + str(f/mm) + 'mm lens; for ' 
             + str(round(wavelength1/nm)) + 'nm', fontsize=18)
ax1= fig.add_subplot(151); ax1.set_title(str(z1/mm) + ' mm'); ax1.imshow(Intensity(1,Fa),cmap='rainbow')
ax2= fig.add_subplot(152); ax2.set_title(str(z2/mm) + ' mm'); ax2.imshow(Intensity(1,Fb),cmap='rainbow')
ax3= fig.add_subplot(153); ax3.set_title(str(z3/mm) + ' mm'); ax3.imshow(Intensity(1,Fc),cmap='rainbow')
ax4= fig.add_subplot(154); ax4.set_title(str(z4/mm) + ' mm'); ax4.imshow(Intensity(1,Fd),cmap='rainbow')
ax5= fig.add_subplot(155); ax5.set_title(str(z5/mm) + ' mm'); ax5.imshow(Intensity(1,Fe),cmap='rainbow')
plt.tight_layout(pad=5)
plt.show()

Phase recovery using Gerchberg Saxton iteration¶

"""
    Phase recovery from two measured intensity distributions using Gerchberg Saxton algorithm
    Adapted from PhaseRecovery.py, copyright: (c) 2017 by Fred van Goor, released under MIT license
    Modifications by Raghav K. Chhetri, September 2020
"""
wavelength= 500*nm
size= 10*mm
N= 1024
z=2*m; #propagation distance between the two fields

# Input field (near)
w0= 1.5*mm
F0= Begin(size,wavelength,N)
F0= GaussAperture(F0, w0)
I0= Intensity(0,F0)

# Output field (far)
n= 2 #number of spots
spots_dx= 1*mm
sizeforspot= spots_dx
w1= 250*um #Gaussian aperture 1/e width

spotgridN= int(N*sizeforspot/size)
Fspot= Begin(sizeforspot,wavelength,spotgridN)
Fspot= GaussAperture(Fspot,w1)

F1= Begin(size,wavelength,N)
F1= FieldArray2D(F1, Fspot, n, n, spots_dx, 3*spots_dx)
I1= Intensity(F1)

#Plots
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
plt.imshow(I0,cmap='jet'); plt.axis('off')
plt.title('Measured Intensity (near)')
plt.subplot(1,2,2)
plt.imshow(I1,cmap='jet'); plt.axis('off')
plt.title('Measured Intensity (far)')
plt.show()

#3d surface of field intensities
NewGridDimension = int(N/8)
F0_interp= Interpol(size,NewGridDimension,0,0,0,1,F0)
I0_interp= Intensity(F0_interp)
F1_interp= Interpol(size,NewGridDimension,0,0,0,1,F1)
I1_interp= Intensity(F1_interp)

X= range(NewGridDimension); Y= range(NewGridDimension)
X, Y= np.meshgrid(X,Y) 

fig= plt.figure(figsize=(18,7)) 
gs= gridspec.GridSpec(1, 2, width_ratios=[1,1]) 
ax1= plt.subplot(gs[0], projection='3d')
ax1.plot_surface(X,Y,I0_interp,rstride=1,cstride=1,cmap='rainbow',linewidth=0.0)
ax1.set_zlabel('Intensity [a.u.]'); ax1.set_title('Field intensity (near)')
ax2= plt.subplot(gs[1], projection='3d')
ax2.plot_surface(X,Y,I1_interp,rstride=1,cstride=1,cmap='rainbow',linewidth=0.0)
ax2.set_zlabel('Intensity [a.u.]'); ax2.set_title('Field intensity (far)')

plt.tight_layout()
plt.show()

Using Forvard for forward and back propagation

N_iterations=100 #number of iterations
N_new= N
size_new= size

#Define a field with uniform amplitude- (1) and phase (0) distribution ie = plane wave
F=Begin(size,wavelength,N)

#The iteration:
for k in range(1,N_iterations):
    #print(k)
    F=SubIntensity(I1,F) #Substitute the measured far field into the field
    F=Interpol(size_new,N_new,0,0,0,1,F);#interpolate to a new grid
    F=Forvard(-z,F) #Propagate back to the near field
    F=Interpol(size,N,0,0,0,1,F) #interpolate to the original grid
    F=SubIntensity(I0,F) #Substitute the measured near field into the field
    F=Forvard(z,F) #Propagate to the far field

#Recovered far- and near field and their phase- and intensity
#distributions; phases are unwrapped (i.e. multiples of pi removed)
F3= F #recovered far field
I3= Intensity(0,F3)
Phase_3= Phase(F3)
Phase_3= PhaseUnwrap(Phase_3)

F2= Forvard(-z,F) #recovered near field
I2= Intensity(0,F2)
Phase_2= Phase(F2)
Phase_2=PhaseUnwrap(Phase_2)

Phase_2= np.interp(Phase_2, (Phase_2.min(), Phase_2.max()), (0, 255))
Phase_3= np.interp(Phase_3, (Phase_3.min(), Phase_3.max()), (0, 255))

#Plot the recovered intensity- and phase distributions
plt.figure(figsize=(18,12))

plt.subplot(2,2,1)
plt.imshow(I2,cmap='jet'); plt.axis('off')
plt.title('Recovered Intensity (near)')

plt.subplot(2,2,2)
plt.imshow(I3,cmap='jet'); plt.axis('off')
plt.title('Recovered Intensity (far)')

plt.subplot(2,2,3)
plt.imshow(Phase_2,cmap='jet'); plt.axis('off')
plt.title('Recovered phase (near)')

plt.subplot(2,2,4)
plt.imshow(Phase_3,cmap='jet'); plt.axis('off')
plt.title('Recovered phase (far)')

plt.show()

#3d surface of field intensities
#NewGridDimension= int(N_new/8)
NewGridDimension= int(N/8)

F2_interp= Interpol(size,NewGridDimension,0,0,0,1,F2)
I2_interp= Intensity(F2_interp)
F3_interp= Interpol(size,NewGridDimension,0,0,0,1,F3)
I3_interp= Intensity(F3_interp)

X= range(NewGridDimension)
Y= range(NewGridDimension)
X, Y= np.meshgrid(X,Y) 

fig= plt.figure(figsize=(18,7)) 
gs= gridspec.GridSpec(1, 2, width_ratios=[1,1]) 
ax1= plt.subplot(gs[0], projection='3d')
ax1.plot_surface(X,Y,I2_interp,rstride=1,cstride=1,cmap='rainbow',linewidth=0.0)
ax1.set_zlabel('Intensity [a.u.]'); ax1.set_title('Recovered Field intensity (near)')
ax2= plt.subplot(gs[1], projection='3d')
ax2.plot_surface(X,Y,I3_interp,rstride=1,cstride=1,cmap='rainbow',linewidth=0.0)
ax2.set_zlabel('Intensity [a.u.]'); ax2.set_title('Recovered Field intensity (far)')

plt.tight_layout()
plt.show()

Many more examples in the LightPipes documentation¶

LighPipes Commands¶

Axicon(Fin, phi, n1=1.5, x_shift=0.0, y_shift=0.0)
BeamMix(Fin1, Fin2)
Begin(size, labda, N)
CircAperture(Fin, R, x_shift=0.0, y_shift=0.0)
CircScreen(Fin, R, x_shift=0.0, y_shift=0.0)
Convert(Fin)
CylindricalLens(Fin, f, x_shift=0.0, y_shift=0.0, angle=0.0)
FieldArray2D(Fin, Ffield, Nfieldsx, Nfieldsy, x_sep, y_sep)
Forvard(Fin, z)
Forward(Fin, z, sizenew, Nnew)
Fresnel(Fin, z)
Gain(Fin, Isat, alpha0, Lgain)
GaussAperture(Fin, w, x_shift=0.0, y_shift=0.0, T=1.0)
GaussBeam(Fin, w0, n=0, m=0, xshift=0, yshift=0, tx=0, ty=0, doughnut=False, LG=False)
GaussHermite(Fin, w0, m=0, n=0, A=1.0)
GaussLaguerre(Fin, w0, p=0, l=0, A=1.0)
GaussScreen(Fin, w, x_shift=0.0, y_shift=0.0, T=0.0)
IntAttenuator(Fin, att=0.5)
Intensity(Fin, flag=0)
Interpol(Fin, new_size, new_N, x_shift=0.0, y_shift=0.0, angle=0.0, magnif=1.0)
Lens(Fin, f, x_shift=0.0, y_shift=0.0)
LensForvard(Fin, f, z)
LensFresnel(Fin, f, z)
MultIntensity(Fin, Intens)
MultPhase(Fin, Phi)
Normal(Fin)
Phase(Fin, unwrap=False, units='rad', blank_eps=0)
PhaseSpiral(Fin, m=1)
PhaseUnwrap(Phi)
PipFFT(Fin, index=1)
PlaneWave(Fin, w, tx=0.0, ty=0.0, x_shift=0.0, y_shift=0.0)
PointSource(Fin, x=0.0, y=0.0)
Power(Fin)
RandomIntensity(Fin, seed=123, noise=1.0)
RandomPhase(Fin, seed=456, maxPhase=3.141592653589793)
RectAperture(Fin, sx, sy, x_shift=0.0, y_shift=0.0, angle=0.0)
RectScreen(Fin, sx, sy, x_shift=0.0, y_shift=0.0, angle=0.0)
RowOfFields(Fin, Ffield, Nfields, sep, y=0.0)
Steps(Fin, z, nstep=1, refr=1.0, save_ram=False, use_scipy=False)
Strehl(Fin)
SubIntensity(Fin, Intens)
SubPhase(Fin, Phi)
SuperGaussAperture(Fin, w, n=2.0, x_shift=0.0, y_shift=0.0, T=1.0)
Tilt(Fin, tx, ty)
Zernike(Fin, n, m, R, A=1.0, norm=True, units='opd')
ZernikeFilter(F, j_terms, R)
ZernikeFit(F, j_terms, R, norm=True, units='lam')
ZernikeName(Noll)
ZonePlate(Fin, N_zones, f=None, p=None, q=None, T=1.0, PassEvenZones=True)
noll_to_zern(j)

Explore further using the built-in help() function

Other interesting open-source Optics packages in Python¶

POPPY: Physical Optics Propagation in Python
Simulation of Fraunhofer and Fresnel diffraction and point spread function formation, particularly in the context of astronomical telescopes

BioBeam
Multiplexed wave-optical simulations of light-sheet microscopy

pyOpTools
Simulation of optical systems by 3D non-sequential ray tracing