function darabs
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Darwin reflectivity curve of Si (111), including absorption

set(gcf,'papertype','a4','paperunits','centimeters','units','centimeters')
set(gcf,'position',[0.1 -8 21 26],'paperposition',[0.1 0.1 21 26])

axes('Position',[0.2 0.60 0.6 0.40])

% Case 1: lambda=1.5405 Angs
r0=2.82E-5;                               % Thompson scattering length in Angs
V=160.1966;                               % unitcell volume in Ang^3 e.g. 160.1966 for Si
d=3.13562;                                % d spacing for Si (111)
m=1;                                      % order of reflection, ie 1 for (111), 3 for (333)
F_hkl=abs(4-4*i)*(10.54+0.25-i*0.33);     % Complex structure factor for 111
F_0=8*(14+0.25-i*0.33);                   % Complex structure factor for 000
g=(2*d*d/m)*(r0/V)*F_hkl;
g0=g*(F_0/F_hkl);
[x,R]=darwin(g,g0,m); line(x,R,'color','b','linestyle','--')

axis([-2 2 0 1.1])
set(gca,'Ytick',[0 0.5 1],'FontSize',20,'FontName','Times')
xlabel('x','position',[0 -0.1 0]); ylabel('Intensity reflectivity','position',[-2.75 0.5 0])
box on; grid on

% Case 2: lambda=0.70926 Angs
F_0=8*(14+0.082-i*0.071);                 % Complex structure factor for 000
F_hkl=abs(4-4*i)*(10.54+0.082-i*0.071);   % Complex structure factor for 111
g=(2*d*d/m)*(r0/V)*F_hkl;
g0=g*(F_0/F_hkl);
[x,R]=darwin(g,g0,m); 
line(x,R,'color','b','linestyle','-');
text(-1.80,0.95,'(a)','FontName','Times','Fontsize',24)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot as a function of energy and angular variable in milli degrees

axes('Position',[0.2 0.12 0.6 0.40])

lambda=12.398/5.000;                      % 5 keV
theta=asin(m*lambda/2/d);
F_hkl=abs(4-4*i)*(10.54+0.38-i*0.8029);   % Complex structure factor for 111
F_0=8*(14+0.3807-i*0.8029);               % Complex structure factor for 000
g=(2*d*d/m)*(r0/V)*F_hkl;
g0=g*(F_0/F_hkl);
[x,R]=darwin(g,g0,m);
line(x*real(g/m/pi)*tan(theta)*180/pi*1e3,R,'color','b','linestyle','--')

lambda=12.398/10.000;                      % 10 keV
theta=asin(m*lambda/2/d);
F_hkl=abs(4-4*i)*(10.54+0.1943-i*0.2169);  % Complex structure factor for 111
F_0=8*(14+0.1943-i*0.2169);                % Complex structure factor for 000
g=(2*d*d/m)*(r0/V)*F_hkl;
g0=g*(F_0/F_hkl);
[x,R]=darwin(g,g0,m);
line(x*real(g/m/pi)*tan(theta)*180/pi*1e3,R,'color','b','linestyle','--')

lambda=12.398/50.000;                         % 50 keV
theta=asin(m*lambda/2/d);
F_hkl=abs(4-4*i)*(10.54+0.0027-i*0.0076);     % Complex structure factor for 111
F_0=8*(14+0.0027-i*0.0076);                   % Complex structure factor for 000
g=(2*d*d/m)*(r0/V)*F_hkl;
g0=g*(F_0/F_hkl);
[x,R]=darwin(g,g0,m);
line(x*real(g/m/pi)*tan(theta)*180/pi*1e3,R,'color','b','linestyle','-')

axis([-3 3 0 1.1])
set(gca,'Xtick',[-2 -1 0 1 2],'Ytick',[0 0.5 1],'FontSize',20,'FontName','Times')
xlabel('\omega (milli degrees)'); 
ylabel('Intensity reflectivity','position',[-4.15 0.55 0])
box on; grid on
text(-2.85,0.25,'5 keV','FontName','Times','Fontsize',18)
text(-2.10,0.15,'10 keV','FontName','Times','Fontsize',18)
text(-1.30,0.05,'50 keV','FontName','Times','Fontsize',18)
text(-2.70,0.95,'(b)','FontName','Times','Fontsize',24)

function [x,R]=darwin(g,g0,m);
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Darwin reflectivity R vs x  (absorption effects included)

x_m=[-5:0.01:-1]; zeta=real((g*x_m+g0)/m/pi); xc_m=m*pi*zeta/g-g0/g;
rc_m=xc_m+sqrt(xc_m.^2-1);
x_t=[-1:0.01:1]; zeta=real((g*x_t+g0)/m/pi); xc_t=m*pi*zeta/g-g0/g; 
rc_t=xc_t-i*sqrt(1-xc_t.^2);
x_p=[1:0.01:5]; zeta=real((g*x_p+g0)/m/pi); xc_p=m*pi*zeta/g-g0/g;
rc_p=xc_p-sqrt(xc_p.^2-1);
x=[x_m x_t x_p]; rc=[rc_m rc_t rc_p];
R=abs(rc).^2;
