function ctr
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Properties of the Crystal Truncation Rod 

figure
set(gcf,'papertype','a4','paperunits','centimeters','units','centimeters',...
        'position',[0.1 -8 21 26],'paperposition',[0.1 0.1 21 26]);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot rod from flat surface without (beta=0) and with (beta=0.2) absorption
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.2 0.55 0.6 0.35])

ell1=[0.01:0.001:0.99];      % beta=0, l range chosen to avoid Bragg peak at l=1
F_CTR=1./(1-exp(i*2*pi*ell1));
h1=line(ell1,F_CTR.*conj(F_CTR),'color','b','linewidth',1,'linestyle','-')
ell2=[1.01:0.001:1.99];
F_CTR=1./(1-exp(i*2*pi*ell2));
line(ell2,F_CTR.*conj(F_CTR),'color','b','linewidth',1,'linestyle','-')

ell=[0.01:0.001:1.99];       % beta=0.2, l range now includes Bragg peak at l=1
beta=0.2;
F_CTR=1./(1-exp(i*2*pi*ell)*exp(-beta));
h2=line(ell,F_CTR.*conj(F_CTR),'color','b','linewidth',1,'linestyle','-.')

set(gca,'Fontsize',16,'FontName','Times')
[h,obj]=legend([h1 h2],'\beta=0','\beta=0.2')
set(gca,'Fontsize',18,'FontName','Times')
set(gca,'FontName','Times','FontSize',18)
xlabel('{\it l} (r.l.u.)'); ylabel('|{\it F }^{CTR}|^2'); box on
axis([0.0 2.0 0.1 1000])
set(gca,'Ytick',[0.1 1 10 100 1000],'Yscale','Log')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot rod from flat surface + overlayer at different relative diplacements, z0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.2 0.12 0.6 0.35])

ell1=[0.01:0.001:0.99];         %l range chosen to avoid Bragg peak at l=1
F_CTR=1./(1-exp(i*2*pi*ell1));
line(ell1,F_CTR.*conj(F_CTR),'color','b','linewidth',1,'linestyle','-')
ell2=[1.01:0.001:1.99];
F_CTR=1./(1-exp(i*2*pi*ell2));
h1=line(ell2,F_CTR.*conj(F_CTR),'color','b','linewidth',1,'linestyle','-')

z0=0.05;                        % relative displacement of overlayer, z0=0.05
F_CTR=1./(1-exp(i*2*pi*ell1));
F_T=F_CTR+exp(-i*2*pi*(1+z0)*ell1);
line(ell1,F_T.*conj(F_T),'color','b','linewidth',1,'linestyle','--')
F_CTR=1./(1-exp(i*2*pi*ell2));
F_T=F_CTR+exp(-i*2*pi*(1+z0)*ell2);
h2=line(ell2,F_T.*conj(F_T),'color','b','linewidth',1,'linestyle','--')

z0=-0.05;                       % relative displacement of overlayer, z0=-0.05
F_CTR=1./(1-exp(i*2*pi*ell1));
F_T=F_CTR+exp(-i*2*pi*(1+z0)*ell1);
line(ell1,F_T.*conj(F_T),'color','b','linewidth',1,'linestyle','-.')
F_CTR=1./(1-exp(i*2*pi*ell2));
F_T=F_CTR+exp(-i*2*pi*(1+z0)*ell2);
h3=line(ell2,F_T.*conj(F_T),'color','b','linewidth',1,'linestyle','-.')

set(gca,'FontName','Times','FontSize',18)
xlabel('{\it l} (r.l.u.)'); ylabel('|{\it F }^{CTR}|^2'); box on
set(gca,'Fontsize',14,'FontName','Times')
[h,obj]=legend([h1 h2 h3],'{\it z}_0=0','{\it z}_0=0.05','{\it z}_0=-0.05');
set(gca,'Fontsize',18,'FontName','Times')
axis([0. 2.0 0.1 1000])
set(gca,'Ytick',[0.1 1 10 100 1000],'Yscale','log')