function DebyeWaller
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates:  The Debye-Waller factor for Aluminium
% Calls to: phiDebye

figure
set(gcf,'papertype','a4','paperunits','centimeters','units','centimeters',...
        'position',[0.1 -8 21 26],'paperposition',[0.1 0.1 21 26]);
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot of phi(x) vs x.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.30 0.70 0.45 0.225]);

x=0.01:0.02:8;
for il=1:length(x), phi(il)=phiDebye(x(il));   end

line(x,phi,'color','b','linewidth',1.5)

axis([0 8 0 1.1]); grid on
set(gca,'FontName','Times','FontSize',16,'box','on')
ylabel('phi(x)','position',[-1.5 0.55 0]); xlabel('Theta/T')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot of sqrt(u^2) vs Temperature for Al
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.30 0.40 0.45 0.225]);

Theta_Al=394;            % Debye temperature of Al 
A=27;                    % Atomic mass
nnd=4.04/sqrt(2);        % Nearest neighbour distance
T=x*394;
B_Al=11492.*T.*phi/A/Theta_Al/Theta_Al+2873/A/Theta_Al;
rms=sqrt(3/8/pi/pi.*B_Al);
iT=find(T<933); iTg=find(T>=933);
line(T(iT),sqrt(2)*rms(iT)/4.04,'color','g','linewidth',1.5)
line(T(iTg),sqrt(2)*rms(iTg)/4.04,'color','g','linewidth',1.5,'linestyle',':')

axis([0 1050 0 0.1]); grid on
set(gca,'FontName','Times','FontSize',16,'box','on')
xlabel('Temperature (K)','position',[500 -0.017 0]) 
ylabel('rms','position',[-200 0.05 0])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Temperature dependence at different Q's
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

axes('position',[0.30 0.10 0.45 0.225]);

I400=exp(-(8/4.04/4.04.*B_Al))./exp(-(8/4.04/4.04.*B_Al(1)));
I800=exp(-(32/4.04/4.04.*B_Al))./exp(-(32/4.04/4.04.*B_Al(1)));

line(T(iT),I400(iT),'color','r','LineWidth',1.5)
line(T(iT),I800(iT),'color','r','Linestyle','--','LineWidth',1.5)

axis([0 1050 0 1.1]); grid on
set(gca,'FontName','Times','FontSize',16,'box','on')
xlabel('Temperature (K)','position',[500 -0.2 0]); 
ylabel('Relative Intensity','position',[-200 0.55 0])
text(650,0.30,'(8,0,0)','FontName','Times','FontSize',16)
text(650,0.70,'(4,0,0)','FontName','Times','FontSize',16)

function phi=phiDebye(x)
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Evaluates the integral to calculate phi(x)
% Calls to: phiDebyeInt

phi=quad8('phiDebyeInt',0.000000001,x)./x;

function y=phiDebyeInt(xi)
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Defines the integrabd used to evaluate phi(x)
% Note: Must be placed in a separate file called phiDebyeInt.m

y=xi./(exp(xi)-1);