function cuoctr
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: CTR of O on Cu (110) and compares with data 
% Calls to: ff
% Data: Feidenhans'l et al., Phys. Rev. B., vol. 41, page 5420 (1990)

%----- Cu real and reciprocal lattice parameters

ac=3.615; ar=2*pi/(ac/sqrt(2)); br=2*pi/ac; cr=2*pi/(ac/sqrt(2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%(1,1) rod
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure; axes('position',[0.55 0.15 0.35 0.8]);

h=1; k=1; l=0.05:0.001:1.0; Q=sqrt(h^2*ar^2+k^2*br^2+l.^2*cr^2);

%----- Cu form factor

a=[13.338 7.1676 5.6158 1.6735]; b=[3.5828 0.2470 11.3966 64.82]; c=[1.1910];
f_Cu=ff(a,b,c,Q);

%------ Cu Debye-Waller factor for bulk (B) and surface (S)

DW_Cu_B=exp(-0.55*(Q/4/pi).^2); DW_Cu_S=exp(-1.70*(Q/4/pi).^2);

%----- O form factor

a=[3.0485 2.2868 1.5463 0.8670]; b=[13.2771 5.7011 0.3239 32.9089]; c=[0.2508];
f_O=ff(a,b,c,Q);

%----- Bulk CTR

Phi=pi*(h+k+l); F_CTR=f_Cu.*DW_Cu_B./(1-exp(i*Phi));
semilogy(l,4*abs(F_CTR).^2,'-.'); 
axis([0 1 6 100000 ])
set(gca,'FontName','Times','FontSize',16,'Ytick',[1e1 1e2 1e3 1e4 1e5])
text(0.1,50000,'(b) (1,1) rod','FontName','Times','FontSize',16)

%----- Add 1/2 a monolayer of Cu (no relaxation) 

F_S=0.5*f_Cu.*DW_Cu_S*exp(i*pi*(h+k)).*exp(-i*2*pi*0.5*l); 
F_T=F_CTR+F_S;
line(l,4*abs(F_T).^2,'linestyle','--');

%----- Add 1/2 a monolayer of Cu (relaxed to z0) plus O layer (relaxed to ´-z1)

z0=0.1445; z1=z0-0.133;
F_S=0.5*exp(i*pi*(h+k))*(f_Cu.*exp(-i*2*pi*(0.5+z0)*l).*DW_Cu_S...
                        +f_O.*exp(i*pi*k).*exp(-i*2*pi*(0.5+z1)*l));
F_T=F_CTR+F_S;
line(l,4*abs(F_T).^2,'linestyle','-');

%----- Add 1/2 a monolayer of Cu (relaxed to z0) plus O layer (relaxed to ´+z1)

z0=0.1445; z1=z0+0.133;
F_S=0.5*exp(i*pi*(h+k))*(f_Cu.*exp(-i*2*pi*(0.5+z0)*l).*DW_Cu_S...
                        +f_O.*exp(i*pi*k).*exp(-i*2*pi*(0.5+z1)*l));
F_T=F_CTR+F_S;
line(l,4*abs(F_T).^2,'linestyle',':');

data=[0.0787 4.5555;0.1517 3.8388;0.2247 3.5032;0.2978 3.1584;0.3539 3.0586;...
      0.3708 2.8954;0.4382 2.6051;0.5169 2.3692;0.5787 2.1878;0.6517 1.9700;...
      0.7360 1.9156];
line(data(:,1),10.^data(:,2),'marker','o','linestyle','none','markerfacecolor','w')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%(1,0) rod
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.10 0.15 0.35 0.8]); 

h=1; k=0; l=0.0:0.001:0.95; Q=sqrt(h^2*ar^2+k^2*br^2+l.^2*cr^2);

%----- Cu form factor

a=[13.338 7.1676 5.6158 1.6735]; b=[3.5828 0.2470 11.3966 64.82]; c=[1.1910];
f_Cu=ff(a,b,c,Q);

%----- Cu Debye-Waller factor for bulk (B) and surface (S)

DW_Cu_B=exp(-0.55*(Q/4/pi).^2); DW_Cu_S=exp(-1.70*(Q/4/pi).^2);

%----- O form factor

a=[3.0485 2.2868 1.5463 0.8670]; b=[13.2771 5.7011 0.3239 32.9089]; c=[0.2508];
f_O=ff(a,b,c,Q);

%----- Bulk unit cell SF

Phi=pi*(h+k+l); F_CTR=f_Cu.*DW_Cu_B./(1-exp(i*Phi));

semilogy(l,4*F_CTR.*conj(F_CTR),'-.'); 
axis([0 1 6 100000 ])
set(gca,'FontName','Times','FontSize',16,'Ytick',[1e1 1e2 1e3 1e4 1e5])
text(0.1,50000,'(a) (1,0) rod','FontName','Times','FontSize',16)
ylabel('Intensity (electron units)','Fontsize',18)
text(1.15,2,'l (r.l.u.)','FontName','Times','FontSize',18,...
           'HorizontalAlignment','Center')

%----- Add 1/2 a monolayer of Cu (not relaxed)

F_S=0.5*exp(i*pi*(h+k))*f_Cu.*DW_Cu_S.*exp(-i*2*pi*(0.5+z0)*l);
F_T=F_CTR+F_S.*DW_Cu_B;
line(l,4*abs(F_T).^2,'linestyle','--');

%----- Add 1/2 a monolayer of Cu (relaxed to z0) plus O layer (relaxed to ´-z1)

z0=0.1145; z1=z0-0.133;
F_S=0.5*exp(i*pi*(h+k))*(f_Cu.*exp(-i*2*pi*(0.5+z0)*l).*DW_Cu_S...
                        +f_O.*exp(i*pi*k).*exp(-i*2*pi*(0.5+z1)*l));
F_T=F_CTR+F_S;
line(l,4*abs(F_T).^2,'linestyle','-');

%----- Add 1/2 a monolayer of Cu (relaxed to z0) plus O layer (relaxed to ´+z1)

z0=0.1145; z1=z0+0.133;
F_S=0.5*exp(i*pi*(h+k))*(f_Cu.*exp(-i*2*pi*(0.5+z0)*l).*DW_Cu_S...
                        +f_O.*exp(i*pi*k).*exp(-i*2*pi*(0.5+z1)*l));
F_T=F_CTR+F_S;
line(l,4*abs(F_T).^2,'linestyle',':');

%----- Plot data

data=[0.0226 1.8793;0.0960 1.9973;0.1695 2.2785;0.2429 2.4781;0.3220 2.7139;...
      0.3955 2.8228;0.4689 2.9861;0.5480 3.1222;0.6158 3.2310;0.6949 3.4124;...
      0.7740    3.6574];
line(data(:,1),10.^data(:,2),'marker','o','linestyle','none','markerfacecolor','w')
