% 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) 
