Propagation of Fresnel diffraction pattern as a function of propagation distance
Matlab code
% Program following the evolution of Fresnel diffraction for a given input
% field 'win' defined by an array of complex numbers A exp(i phi).
% Uses the function 'prop_free_nf.m' from the CXS group, cSAXS
% beamline, Paul Scherrer Institute. See
% https://www.psi.ch/en/sls/csaxs/software, and in this, the
% file path cSAXS_matlab_base_package/+utils
clear
close all
% _________________________________________________________________________
% Prompt user to input parameters
prompt = 'Square (S) or circular (C) object [S]: ';
shapeType = input(prompt,'s');
if isempty(shapeType)
shapeType = 'S';
elseif (shapeType == 'S') || (shapeType == 's')
shapeType = 'S';
else
shapeType = 'C';
end
prompt = 'Linear length of detector in number of pixels [1024]: ';
linNoPix = input(prompt);
if isempty(linNoPix)
linNoPix = 1024;
end
win = ones(linNoPix);
prompt = 'Diameter of circular object or length of square object in pixels [100]: ';
width = input(prompt);
if isempty(width)
width = 100;
end
halfWidth = width/2; % If circle chosen, = radius
prompt = 'Transmission of object (0 < T < 1)[1]: ';
A = input(prompt);
if isempty(A)
A = 1;
end
prompt = 'Phase shift induced by object in radians [pi]: ';
objPhase = input(prompt);
if isempty(objPhase)
objPhase = pi;
end
prompt = 'Wavelength of radiation in Angstroms [2]: ';
lambda = input(prompt);
if isempty(lambda)
lambda = 2;
end
lambda = lambda*1e-7; % Convert from Angstroms to mm
prompt = 'Pixel size in mm [0.01]: ';
pixelSize = input(prompt);
if isempty(pixelSize)
pixelSize = 0.01;
end
% Set up field describing the diffracting object
if (shapeType == 'S') || (shapeType == 's')
% Range of central part of field in which the object lies
phaseRange = linNoPix/2 - halfWidth:linNoPix/2 + halfWidth;
% Central square diffracting object has its own phase and amplitude
win(phaseRange,phaseRange) = A*exp(1i*objPhase);
else
for jj = 1:linNoPix
for kk = 1:linNoPix
if ((jj-linNoPix/2)^2 + (kk-linNoPix/2)^2 <= halfWidth^2)
win(jj,kk) = A*exp(1i*objPhase);
end
end
end
end
% End of prompting user for inputs
% _________________________________________________________________________
vid = VideoWriter('fresnelPropagation.mp4','MPEG-4');
vid.FrameRate = 30; % Default 30
vid.Quality = 100; % Default 75
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]); % White background
set(gca,'linewidth',7);
% _________________________________________________________________________
% Set up strings for printing
str1 = 'Fresnel # = ';
str3 = 'z = ';
str5 = ' m';
str6 = 'Transparency = ';
str7 = num2str(A*100,'% .0f');
str8 = '%';
strTransp = [str6,str7,str8];
if (shapeType == 'S') || (shapeType == 's')
str9 = 'Square, size = ';
str10 = num2str(pixelSize*halfWidth*2,'% .2f');
else
str9 = 'circle, radius = ';
str10 = num2str(pixelSize*halfWidth,'% .2f');
end
str11 = ' mm';
strSqSize = [str9,str10,str11];
str12 = 'Phase = ';
str13 = num2str(objPhase/pi,'% .2f');
str14 = ' pi';
strPhase = [str12,str13,str14];
energy = 12.3984/(lambda*10^7); % Photon energy in keV
str15 = 'h\nu = ';
str16 = num2str(energy,'% .3f');
str17 = ' keV';
strEnergy = [str15,str16,str17];
% End setting up strings for printing
% _________________________________________________________________________
myGray = customcolormap([0 0.75 1], {'#ffffff','#808080','#000000'});
loopRange = 2.5:0.005:5.5; % Exponent for z = 10^loopRange in mm
ydim = size(loopRange,2);
xdim = 1+linNoPix/2;
% Initialize plot of central horizontal slice through 2D pattern as a
% function of propagation distance
FresnelCrossSectionVz = zeros(xdim,ydim);
ydimIndex = 0;
for jj = loopRange % Determines range of propagation lengths to be sampled
ydimIndex = ydimIndex + 1;
z = 10^jj; % in mm
% _____________________________________________________________________
% Core line of program!! See https://www.psi.ch/en/sls/csaxs/software
% Thanks to the CXS group and the cSAXS beamline, Paul Scherrer
% Institute
wout = prop_free_nf(win,lambda,z,pixelSize);
% _____________________________________________________________________
imshow((abs(wout)).^2,[0 4])
axis equal
axis off
colormap(myGray)
% Print present Fresnel number
Fnum = (pixelSize*halfWidth)^2/(z*lambda);
str2 = num2str(Fnum,'% .2f');
strTot = [str1,str2];
annotation('textbox',[0.625 0.17 0.14 0.07],'String',strTot,'EdgeColor',...
'none','FontSize',11,'FitBoxToText','on','HorizontalAlignment',...
'center','verticalAlignment','bottom','Color','w');
% Print present propagation length z in meters
str4 = num2str(z/1000,'% .2f'); % String of z in meters
strTot = [str3,str4,str5];
annotation('textbox',[0.625 0.14 0.14 0.07],'String',strTot,'EdgeColor',...
'none','FontSize',11,'FitBoxToText','on','HorizontalAlignment',...
'center','verticalAlignment','bottom','Color','w');
% Print properties of diffracting square and wavelength
annotation('textbox',[0.625 0.11 0.14 0.07],'String',[strSqSize,', ',strTransp],'EdgeColor',...
'none','FontSize',11,'FitBoxToText','on','HorizontalAlignment',...
'center','verticalAlignment','bottom','Color','w');
annotation('textbox',[0.625 0.08 0.14 0.07],'String',[strPhase,', ',strEnergy],'EdgeColor',...
'none','FontSize',11,'FitBoxToText','on','HorizontalAlignment',...
'center','verticalAlignment','bottom','Color','w');
frame = getframe(gcf);
writeVideo(vid,frame);
FresnelCrossSectionVz(:,ydimIndex) = (abs(wout(linNoPix/4:3*linNoPix/4,linNoPix/2))).^2;
end
imshow(FresnelCrossSectionVz',[0 4])
colormap(myGray)
saveas(gcf,'FresnelCrossSectionVz.png');
axis square
axis off
close(vid)