lens Following Field Produced by Forward Function

[17aacw]

I'm trying to model a 1000 nm wavelength columnated beam of width 1 cm travelling through a 10 cm focal length lens to a mirror located at f and then reflected through the 10 cm lens. Despite its computational expense, I've found I get the best results returning to the lens using the Forward function. However, I've found that the phase of the field from the Forward function does not match the input beam at the lens. Further, I've found that the beam produced by the Forward function does not behave as expected when propagated back through the lens.

Are there any recommendations on how to fix this system or fix this problem?

Attached is an example of the system output with the beam going into the system in blue and the returning beam in orange. The code used to produce this graphic is also attached.

Thanks for your time and assitance

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from LightPipes import *

#%% Setting Constants
#Trial Peramiters
f=10*cm #focal length of lens

zl=80cm-f #lens location
w0 = 1cm #Waste of the initial beam
lm = f #distance to mirror from lens

#starting conditions
wavelength=1um
StartSize=w05
Nstart=128
Nstart2=int(Nstart/2)
z0 = np.pi*w0**2/wavelength

#zoomed in condition
w=2um
size= w5
N=64

#%% Function for interpuplation and forvard
def InterpForvard(F0,z0,zfinal,w0,size):
N=64
dz=z0*16
zlist = np.linspace(0,zfinal,round(zfinal/dz)+1)

if zfinal < dz: 
    zlist = [zlist]
    
for z in zlist[1:len(zlist)]: 
       wnew=w0*np.sqrt(1+z**2/z0**2)
       sizenew=wnew*5
       F=Forward(z,sizenew,256,F0)

return [F,Nnew]

#%% intiate the beam
#input Beam
F0=Begin(StartSize,wavelength,Nstart,dtype=np.complex64)
F1=GaussBeam(F0,w0)
I1= Intensity(F1)
P1 = Phase(F1)

#Beam going into lens
FF02 = F1._w0**2/(zl*wavelength)
F2=Forvard(F1,zl)
I2 = Intensity(F2)
P2= Phase(F2)

#%% Beam at mirror
#Define beam spot size of thin lens from equation 7.6.17
Lw0 = np.sqrt(F2._w0**2*((1-lm/f)**2+(lm/z0)**2))

z0l = np.pi*Lw0**2/wavelength

sizeL0 = Lw05
#intensity at focus
Fm=Begin(sizeL0,wavelength,N,dtype=np.complex64)
F3=GaussBeam(Lw0,Fm)
F3._field=F3._fieldw0/Lw0 #normalise of power to area. Check this value is constant
I3 = Intensity(F3)
P3=Phase(F3)

#%%reflected beam to lens
[F4,Nnew]=InterpForvard(F3,z0l,lm,Lw0,sizeL0)

I4 = Intensity(F4)
P4 = Phase(F4)

#%% Beam at start
F5=Lens(f,F4)
F5=Forvard(zl,F5)
I5 = Intensity(F5)
P5 = Phase(F5)