function par_kin
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Parratt and kinematical reflectivities from a multilayer
%             Specific case of W/Si, 10 bilayers of [10 Angs W, 40 Angs Si]

r0=2.82e-5;                             % Thompson scattering length in Angs
Q=0.01:0.001:0.3;                       % Wavevector transfer in 1/Angs
lambda=1.54;	                         % wavelength in Angs
rhoA=4.678; muA=33.235e-6;              % density and absorption coefficient of W
rhoB=0.699; muB=1.399e-6;               % density and absorption coefficient of Si

bl=[rhoA*r0+i*muA rhoB*r0+i*muB];                 % bilayer scattering factor   
dbl=[10 40];                                      % bilayer d-spacings
ml=[bl bl bl bl bl bl bl bl bl bl 0.1e-20];       % multilayer scattering factor 
dml=[dbl dbl dbl dbl dbl dbl dbl dbl dbl dbl];    % multilayer d-spacings
sml=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];  % roughness at each interface

%------ Parratt reflectivity

R=parratt(Q,lambda,ml,dml,sml);

axes('position',[0.2 0.15 0.7 0.4]); line(Q,R)
axis([0 0.3 8e-6 2]); 
set(gca,'FontName','Times','FontSize',18,'box','on')
set(gca,'Ytick',[1e-4 1e-3 1e-2 1e-1 1],'Yscale','log')
text(0.15,5e-7,'Wavevector transfer Q (Å^{-1})',...
	'FontName','Times','FontSize',18,'horizontalalignment','center')
text(-0.05,1,'Reflectivity','FontName','Times',...
   'FontSize',18,'horizontalalignment','center','rotation',90)
text(0.20,0.7,'(b) Parratt','FontName','Times','FontSize',16)

%----- kinematical reflectivity

sld=bl;
sigma=0; N=10; Lambda=50; Gamma=0.2;

R=kinematicalR(Q,lambda,sld,sigma,N,Lambda,Gamma);

axes('position',[0.2 0.55 0.7 0.4]); line(Q,R)
set(gca,'FontName','Times','FontSize',18,'box','on')
axis([0 0.3 8e-6 2]);
set(gca,'Xticklabel','')
set(gca,'Ytick',[1e-4 1e-3 1e-2 1e-1 1],'Yscale','log')
text(0.20,0.7,'(a) Kinematical ','FontName','Times','FontSize',16)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [R]=kinematicalR(Q,lambda,sld,sigma,N,Lambda,Gamma)
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: kinematical reflectivity of a multilayer
% Inputs:  Q       wavevector transfer            1/Angs
%          lambda  wavelength of radiation        Angs
%          sld     scattering length density      1/Angs^2
%                  sld=[sldA+i*muA sldB+i*muB]     
%          sigma   rouhgness                      Angs
%          N       number of bilayers
%          Lambda  length of bilayer              Angs
%          Gamma   fraction of bilayer that is A
% Outputs: R       Intensity reflectivity

muA=imag(sld(1));
muB=imag(sld(2));

Dsld=real(sld(1))-real(sld(2));

zeta=Q/2/pi*Lambda;

beta=2*Lambda*Lambda*(muA*Gamma+muB*(1-Gamma))/lambda./zeta;

r_1=-2*i*Dsld*Lambda*Lambda*Gamma./zeta;
r_1=r_1.*sin(pi*Gamma*zeta)./(pi*Gamma*zeta);

r_N=r_1.*(1-exp(i*2*pi*zeta*N).*exp(-beta*N))./(1-exp(i*2*pi*zeta).*exp(-beta));

r_N=r_N.*exp(-((Q*sigma).^2/2));
R=r_N.*conj(r_N);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [RR]=parratt(Q,lambda,sld,d,sigma)
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Parratt reflectivity of a multilayer
% Inputs: Q      wavevector transfer                1/Angs
%         lambda wavelength of radiation            Angs
%         sld    scattering length density          1/Angs^2
%                sld=[sld1+i*mu1 sld2+i*mu2 ....]
%         d      thickness of layer                 Angs
%                d=[d1 d2 .....];
%         sigma  rouhgness                          Angs
% Outputs:R      Intensity reflectivity

k=2*pi/lambda;

%----- Calculate refractive index n of each layer

delta=lambda^2*real(sld)/(2*pi); beta=lambda/(4*pi)*imag(sld);
n=size(sld,2);
nu=1-delta+i*beta;

%----- Wavevector transfer in each layer

Q=reshape(Q,1,length(Q));
x=asin(Q/2/k);
for j=1:n
   Qp(j,:)=sqrt(Q.^2-8*k^2*delta(j)+i*8*k^2*beta(j));   
end
Qp=[Q;Qp];

%----- Reflection coefficients (no multiple scattering)

for j=1:n
   r(j,:)=((Qp(j,:)-Qp(j+1,:))./(Qp(j,:)+Qp(j+1,:))).*...
      exp(-0.5*(Qp(j,:).*Qp(j+1,:))*sigma(j)^2);
end

%----- Reflectivity from first layer

RR=r(1,:);
if n>1
   R(1,:)=(r(n-1,:)+r(n,:).*...
      exp(i*Qp(n,:)*d(n-1)))./(1+r(n-1,:).*r(n,:).*exp(i*Qp(n,:)*d(n-1)));
end

%----- Reflectivity from more layers

if n>2
  for j=2:n-1
     R(j,:)=(r(n-j,:)+R(j-1,:).*...
        exp(i*Qp(n-j+1,:)*d(n-j)))./(1+r(n-j,:).*R(j-1,:).*exp(i*Qp(n-j+1,:)*d(n-j)));
  end
end

%------ Intensity reflectivity

if n==1
  RR=r(1,:);
else
  RR=R(n-1,:);
end

RR=(abs(RR).^2)';





