Split step Fourier propagation - Beam propagation method

Question

Hello I'm having a difficulty using the split step Fourier method. Assuming I want to propagate a Gaussian in free space, I'm supposed to use:

A(x,z) = F^-1 [exp((i*k^2*z)/(2*k_0))* F[A(x,0)]]

where F is the Fourier and F^-1 is the inverse Fourier.

What I want to do is plot I(x,z=0), I(x,z=3) and the intensity distribution in the x-z plane.

I tried doing it numerically (plotting I(x,z=0), I(x,z=3)) using the following code:

lambda = 0.5*10^-6;
k0 = 2*pi/lambda;
w = 10*10^-6;
N=500;
a=0.4*10^-4;
dx=a/N;
x = -a/2:dx:a/2-dx;
Dk_x = 2*pi*N/a;
dk_x=2*pi/a;
k_x=-Dk_x/2:dk_x:Dk_x/2-dk_x;
N = (k_x.^2)/(2*k0);
z = 0:(5*10^-3)/length(N):5*10^-3;
z(end) = [];
% A0 = A(x,z=0)
A0 = exp(-x.^2/w^2);
I_0 = A0.*conj(A0);
% F_A0 is the fourier of A0
F_A0 = fft(A0);
% A3 = A(x,z=3)
A3 = ifft(exp(1i*N*3).*F_A0);
I_3 = A3.*conj(A3);
figure
plot(x,I_3,x,I_0)

However, I_3 is not what I expected to receive which is another Gaussian of smaller intensity.

Also, I'm unsure how I'd plot the intensity distribution in the x-z plane. They suggest using the imagesc function which I suppose I'd use it like:

imagesc(x,z, abs(ifft(exp(1i*N.*z).*F_A0).^2)

but the third argument, which is supposed to be a matrix, is a vector in what I've written..

Can anyone please help me with this?

Thank you in advance.

the propagation is along z, the gaussian you want to propagate is 1D, the split step needs a for loop where you propagate each time with a small dz... — bla, Nov 28 '17 at 21:36
I think I understand what you mean now. I shall try now and update. thank you! — FlyGuy, Nov 28 '17 at 21:55
I think you mixed up the spatial dimension of the gaussian x, and the propagation dimension z. z should be someything like `z=linspace(0,3,Nz)` — bla, Nov 28 '17 at 21:57
From what I understand now, I think the whole z component that I've written is quite random, and I_3 doesn't really represent the the intensity at the output face, since all I've done was take the initial gaussian, and just z=3 and plug them into the equation called A3. I didn't really propagate the beam. What I'm going to try and do now is iterate over z like you said and take the step before as the psi0 each time. I hope this'll work. I will update once I'm done. — FlyGuy, Nov 28 '17 at 22:07

bla · Accepted Answer · 2017-11-28T22:20:50.903

2

here's an example for a split step free Gaussian propagation:

N=2^9; % x grid points
L=100; % box length
dx=L/N; %position grid interval
x=(-L/2+1/N):dx:L/2; %define position grid (centered around origin)

dk=2*pi/L; %momentum grid interval
k=(-N/2+1:1:N/2).*dk; %define momentum grid

A_z=exp(-x.^2); % initial gaussian
dz=0.01; % propagation step
z=0:dz:3; % propagation vector

% do the propagation using split step
for n=1:numel(z)
    A_z=ifft(fftshift(  exp(1i*k.^2*dz).*fftshift(fft(A_z))  ));
    I_z(n,:)=abs(A_z).^2;
end

imagesc(x,z,I_z)
xlim([-20 20]);
xlabel('x'); ylabel('z')

edited Nov 28 '17 at 22:20

answered Nov 28 '17 at 22:13

bla

25,846
10
70
101

That's amazing. Thank you so much. In what I tried doing now I did almost the same thing but in my loop I happened to save only the last iteration of the intensity and so I was receiving a vector of I(x,z=3).. I tried doing something with your code and plugging in the numbers (wavelength, k0, etc) from my original code above and I'm just getting a green screen. What could be the problem here? – FlyGuy Nov 28 '17 at 22:55
you need to use consistent units (say microns) , I dont know what numbers you actually use and why, but check that carefully. I'd start by checking that what you input is indeed a non zero vector. From seeing your example I have a feeling (and forgive me for being blunt) that you dont understand what these parameters mean. – bla Nov 28 '17 at 23:04
I think I got it working now. Thank you so much for the help bla! – FlyGuy Nov 28 '17 at 23:26
glad to help, please "accept" the answer... (on the left side) – bla Nov 28 '17 at 23:30
Another question, if I may. I'm having trouble propagating the beam in a medium that has a refractive index that is both dependent of x and z - n(x,z)=1+f(x,z) where f(x,z) is 1 when |x|<5/(0.5+z/5000) and 0 otherwise. I know how to propagate it when the refractive index is dependent only of x - n(x). There's an extra phase before the ifft of exp(-i*k0/n0*dn*dz). However, I'm having trouble with this. I thought of running two loops for each x and z and fill up a matrix for n(x,z) but then that wouldn't help with it. Any idea? – FlyGuy Nov 29 '17 at 00:37
the loop is only on z, think that you can replace z with time. x is your spatial coordinate. and the "splitting" takes your f(x) to f(k) and back to f(x) using fft because the dispersion is expressed in k space. So for implementing a z-dependant operation in x coordinate, you need an intermediate step, where you split the operator to have the forms `exp(i k^2/2 dz/2) exp(i f(x,z) dz) exp(ik^2/2 dz/2)` – bla Nov 29 '17 at 00:52
I'm not sure I understand.. I'm not dealing with free space anymore. I now have an extra term for the index of refraction of the form n=n0+dn(x). The general solution taken to this is: A_z=exp(-i*dn*dz).*ifft(fftshift( exp(1i*k.^2*dz).*fftshift(fft(A_z)) )); However, dn is not dn(x) but dn(x,z). I'm not sure why what you suggested is how it needs to be done. How would I deal with exp(-i*dn*dz)? The same way (exp(-i*dn(x,z)*dz))? – FlyGuy Nov 29 '17 at 01:12
Also, how will I check this f(x,z) and where? I mean, over each iteration, it's going to have the same value since the condition would be something like abs(x) < 5/(0.5+dz/5000), right? – FlyGuy Nov 29 '17 at 01:17

Split step Fourier propagation - Beam propagation method

1 Answers1