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)