%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the phonon dispersion curve for zinc-blend structure %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dispersion
new;
clf;
clear;
global dm;

UENERGY_IN_JOULE = 1.602192E-19;
ULENGTH_IN_METER = 1.E-10;
UMASS_IN_KG = 1.66053E-27;
UTIME_IN_SECOND = ULENGTH_IN_METER*sqrt(UMASS_IN_KG/UENERGY_IN_JOULE);
UFREQ_IN_HZ = 1./UTIME_IN_SECOND/2/pi;
UFREQ_IN_THZ = UFREQ_IN_HZ / 1E12;

load dm.out;
A = 4*pi/dm(1,2);
HA = A/2.;
QA = A/4.;

GAMMA = [0 0 0];
X = [0 0 HA];
G1 = [HA HA -HA];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The X on the left-hand side is (0, 0, 2*pi/a).           %
% The X on the right-hand side is its translation by a     %
% reciprocal vector (2*pi/a, 2*pi/a, -2*pi/a). So actaully %
% X-K-Gamma is in one straight line.                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X2 = X+G1;
K = [QA QA 0]*1.5;
L = [QA QA QA];


KD = 0:50+18+53+43;

f = phonon(GAMMA);
for I = 1:50+18+53+43,
 if I<=50  
  q = GAMMA + (X-GAMMA)*I/50.;
 elseif I<=50+18
  q = X2 + (K-X2)/18.*(I-50);
 elseif I<=50+18+53
  q = K + (GAMMA-K)/53.*(I-50-18);
 else
  q = GAMMA + (L-GAMMA)/43.*(I-50-18-53);
 end
f = [f phonon(q)];
end

hold on;
for branch = 1:6,
 plot (KD, f(branch,:)*UFREQ_IN_THZ);
end;

ylabel ('THz');
stitle ('Phonon Dispersion Curve for \beta -SiC');

v = axis;
axis ([0 50+18+53+43 v(3) v(4)]);

yline = v(3):(v(4)-v(3))/199:v(4);

plot(ones(200,1)*50,yline,'w');

plot(ones(200,1)*(50+18),yline,'w');

plot(ones(200,1)*(50+18+53),yline,'w');

plot(ones(200,1)*(50+18+53+43),yline,'w');

plot(1:50+18+53+43, ones(50+18+53+43,1)*v(4),'w');

sxlabel ('\Gamma                \Delta           X          K            \Sigma               \Gamma           \Lambda             L');

%%%%%%%%%%%%%%%%%
% FOR HARD COPY %
%%%%%%%%%%%%%%%%%
ss=input('Print on printer NUKE: ','s');
if isstr(ss)
 for branch = 1:6,
  plot (KD, f(branch,:)*UFREQ_IN_THZ,'w');
 end
 whitebg ('w');
 unix ('xdpr -Pnuke &');
end
end;
 

