%LightPipes Simulation with LightPipes for Matlab
%February 2014. F.A. van Goor.
%longslit.m
%Diffraction at long slit example.

clear all;

m=1;
nm=1e-9*m;
mm=1e-3*m;
cm=1e-2*m;

lambda=1000*nm; %wavelength
size=50*mm;
N=300;
a=1.5*mm; %width of the slit

Field=LPBegin(size,lambda,N); %make plane wave
Field=LPRectAperture(a,size,0,0,0,Field);
I0=LPIntensity(2,Field);
%
% To calculate the intensity distribution in the focal plane of a lens we apply a trick to
% enhance the accuracy of the calculation. Because the extension of the field in the focal
% plane is very small compared to the grid dimension the number of grid points in the
% focus to describe the field is very small. Fortunately LightPipes offers a method to
% overcome this problem using 'spherical coordinates', implemented in the
% LPLensForvard and LPLensFresnel commands. When one of these comands is called
% it 'bends' the coordinate system in such a way that it follows the divergent or
% convergent wave front and propagates the field to a distance z in the transformed
% coordinates. The resulting field fits in a reduced (converging beam) or in a increased
% (diverging beam) grid size but with the same number of grid points. The new grid size
% can be extracted from the resulting field.
%
%In what follows we calculate the diffraction to the focus of a lens with
%focal length:
f=1*m;
%It is the combination of a weak lens,
f1=15*m;
%followed by a strong geometrical coordinate transform using spherical coordinates
%with:
f2=(1/f-1/f1)^-1;
%
Field=LPLens(f1,0,0,Field);
Field=LPLensFresnel(f2,f,Field);
%The new grid size is:
size_new=Field(1,N+1);
% To plot the intensity distribution in the focal plane of the lens we use a saturable gain
% sheet to suppress the very intens central lobe. The effect is that the side-lobes will be
% amplified more than the central part of the field.
Field=LPGain(40,10,1,Field);
I=LPIntensity(2,Field);

%Plot the results:
h=figure(1);
set(h,'Position',[680 558 800 420])
subplot(1,2,1);
imagesc(I0); %plot intensity behind slit
axis off equal;colormap(gray);
title('Intensity dstribution just after the slit');
subplot(1,2,2);
imagesc(I); %plot intensity at focus
Str=sprintf('Intensity distribution in the focal plane of a 1m lens.\nNote that the grid size is reduced to: %4.3f mm',size_new/mm);
title(Str);
axis off equal;colormap(gray);
