function kedge
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Photoelectric absorption cross-section of Kr - comparison
%             between hydrogen-like model of K shell contribution and the
%             self-consistent Dirac-Hartree-Fock theory
%             (C.T. Chantler, J. Phys. Chem. Ref. Data vol. 24, 71 (1995))
% Calls to: loaddata

[en,sigma,thom_comp]=loaddata;               % load theoretical values

%----- Plot theoretical photoelectron and (Thomson+Compton) cross-sections 

line(en,sigma,'linestyle','-.'); line(en,thom_comp,'linestyle','--')

%----- Find L shell photoelectric contribution below 14.3 keV

xl=find(en<14.3); xg=find(en>14.32);
xlog=log10(en(xl));
ylog=log10(sigma(xl));
[P,S]=polyfit(xlog,ylog,1);
ylogfit=polyval(P,log10(en));
y=10.^ylogfit;
line(en(xg),y(xg),'linestyle',':')           % Plot L contribution for E>14.32

%----- Plot theorectical absorption for K shell from Stoppe theory

en=en(xg); y=y(xg);
r0=2.82e-5;                                  % Thomson scattering length
ek=14.32;                                    % K edge of Kr energy in keV
lambda=12.398./en;
xi=sqrt(ek./(en-ek));
f=2*pi*sqrt(ek./en).*exp(-4*xi.*acot(xi))./(1-exp(-2*pi*xi));
sigmaa=256/3.*lambda.*(ek./en).^2.5.*f*r0*1e8;
line(en,sigmaa+y,'color','r','linewidth',2)
axis([5 40 50 100000])
set(gca,'FontName','Times','FontSize',18,'Xtick',[5 10 20 50],'Xscale','log','Yscale','log')
grid on; box on
xlabel('Photon energy [keV]');ylabel('Absorption cross-section [barn]')

text(6,150,'Thomson+Compton','FontName','Times','Fontsize',18)
text(10,20000,'L edges','FontName','Times','Fontsize',18,'horizontalalignment','center')
text(25,20000,'K + L edges','FontName','Times','Fontsize',18,'horizontalalignment','center')

function [en,sigma,thom_comp]=loaddata
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Photoelectric and Thomson+Compton cross-sections for Kr

en=[5.3 5.66 6.05 6.47 6.92 7.39 7.9 8.45 9.03 9.65 10.3 11 11.8 12.6 13.5 14 14.3...
   14.3 14.4 14.6 15.4 16.5 17.6 18.8 20.1 21.5 23 24.6 26.3 28.1 30 32.1 34.3 36.7 39.2];
sigma=[3.83e+004 3.20e+004 2.67e+004 2.23e+004 1.87e+004 1.57e+004 1.31e+004 1.09e+004...
       9.04e+003 7.43e+003 6.11e+003 5.03e+003 4.14e+003 3.42e+003 2.82e+003 2.51e+003...
       2.40e+003 2.37e+003 1.81e+004 1.74e+004 1.49e+004 1.25e+004 1.05e+004 8.87e+003...
       7.42e+003 6.20e+003 5.17e+003 4.32e+003 3.60e+003 3.00e+003 2.48e+003 2.05e+003...
       1.69e+003 1.40e+003 1.15e+003];
 thom_comp=[...
       4.72e+002 4.49e+002 4.27e+002 4.04e+002 3.82e+002 3.61e+002 3.40e+002 3.20e+002...
       3.00e+002 2.81e+002 2.63e+002 2.45e+002 2.29e+002 2.13e+002 1.98e+002 1.89e+002...
       1.86e+002 1.85e+002 1.84e+002 1.81e+002 1.71e+002 1.58e+002 1.47e+002 1.36e+002...
       1.26e+002 1.16e+002 1.07e+002 9.92e+001 9.18e+001 8.49e+001 7.85e+001 7.27e+001...
       6.73e+001 6.24e+001 5.80e+001]; 