clf;
m = 50;
range = 6;

w = m-range:0.02:m+range;

plot(w,N(m,m,w));
hold on;
plot(w,N(m,m-2,w),'.');
plot(w,N(m,m+2,w));
plot(w,N(m,m-1,w),'-.');
plot(w,N(m,m+1,w),'--');
title ('Noise functions for the m=50th channel')
axis ([44 56 -0.021 0.052])
text(47.3919-0.07,-0.0062,'m-2');
text(48.5836-0.13,-0.0176,'m-1');
text(50.8284-0.05,-0.0176,'m+1');
text(51.8538-0.05,-0.0062,'m+2');
text(50.5958-0.1,0.0314,'m: main peak');
xlabel ('channel number');
print noise_functions.ps
ss = input('next... ');

clf
plot(w,N(m,m,w))
hold on
plot(w,N(m,m-2,w)+N(m,m+2,w)-N(m,m-4,w)-N(m,m+4,w),'.');
plot(w,N(m,m-1,w)+N(m,m+1,w)-N(m,m-2,w)-N(m,m+2,w),'--')
axis ([44 56 -0.013 0.051])
title ('Combination of noise functions')
text(51.1-0.18, 0.0168, '(odd+even)')
text(45.3025-0.3, -0.0046,'(only even)')
xlabel ('channel number');
print noise_combinations.ps
ss = input('next... ');

clf;
mmax = 300;
mint = mmax*4;
d_phase = 2*pi/mint;
ddsin = sin(d_phase*(1:mint));
ddcos = cos(d_phase*(1:mint));
% ++--++--++-- series
twist = 0;
AR(1) = 1;
AR(2) = 1;
for m = 3:mmax,
 AR(m) = -AR(m-2);
end;
% imaginary part in alternating parity
scale = 1
parity = 1;
for m=1:mmax,
 AI(m) = parity*scale*AR(m);
 parity = -parity;
end;
order = (1:mmax).^0;
for ip = 1:mint,
 ipp = mod((1:mmax)*ip-1,mint)+1;
 dsin = ddsin(ipp).*order;
 fr(ip) = dsin*AR';
end;
plot(1:mint,fr);
title (['perturbation force: ' num2str(mmax) ' channels, no "spin-wave" encoding']);
xlabel('time step')
ylabel('force')
print force_no_twist.ps
ss = input('next... ');

clf;
% spin wave rotation 
twist = 20;
damax = 2*pi/twist;
phi(1) = 0;
for m = 1:mmax,
 BR(m) = AR(m)*cos(phi(m))-AI(m)*sin(phi(m));
 BI(m) = AR(m)*sin(phi(m))+AI(m)*cos(phi(m));
 if rand(1)>0.5 
 phi(m+1) = phi(m)+damax;
 else 
 phi(m+1) = phi(m)-damax;
 end;
end;
order = (1:mmax).^0;
for ip = 1:mint,
 ipp = mod((1:mmax)*ip-1,mint)+1;
 dsin = ddsin(ipp).*order;
 fr(ip) = dsin*BR';
end;
plot(1:mint,fr);
hold on;
title (['perturbation force: ' num2str(mmax) ' channels, 2*pi/' num2str(twist) ' "spin-wave" encoding']);
xlabel('time step')
ylabel('force')
print force_twist.ps
ss = input('next... ');


clf;
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 dos_64_gamma.out;
dat = dos_64_gamma;
freq = 1:dat(2);
 bar((freq-0.5)*dat(3)*UFREQ_IN_THZ,dat(4:3+dat(2)));
 title ('Direct diagonalization of a 64-particle SiC supercell with k=0');
xlabel ('THz');
hold on
axis ([0 35 0 28]);
print bar_gamma_64.ps
ss = input('next... ');

clf;
clear dat, freq;
load dos_64_discrete.out;
dat = dos_64_discrete;
freq = 1:dat(2);
 bar((freq-0.5)*dat(3)*UFREQ_IN_THZ,dat(4:3+dat(2)));
 title ('Direct diagonalization of a 64-particle SiC supercell under a special k-point');
xlabel ('THz');
hold on
axis ([0 35 0 25]);
print bar_discrete_64.ps
ss = input('next... ');

clf;
load dos_64_si.out;
load dos_64_c.out;
plot (dos_64_si(:,1),(dos_64_si(:,2)+dos_64_c(:,2))/2*sum(dat(4:3+dat(2)))*(dos_64_si(200,1)-dos_64_si(199,1)),'-.');
hold on
plot (dos_64_si(:,1),(dos_64_si(:,2)+dos_64_c(:,2))/2*sum(dat(4:3+dat(2)))*(dos_64_si(200,1)-dos_64_si(199,1)),'o');
title ('Using our method: base frequency = 0.11 THz')
xlabel ('THz');
axis ([0 35 0 25]);
print our_64.ps
ss = input('next... ');

clf;
load dos_64_15_si.out;
load dos_64_15_c.out;
plot (dos_64_15_si(:,1),(dos_64_15_si(:,2)+dos_64_15_c(:,2))/2*sum(dat(4:3+dat(2)))*(dos_64_15_si(8,1)-dos_64_15_si(7,1))/5,'*');
title ('Using our method: base frequency = 1.5 THz')
xlabel ('THz');
axis ([0 35 -0.5 25]);
print our_64_1.5.ps
ss = input('next... ');


clf;
load dos_si.out;
load dos_c.out;
freq = 1:dos_si(2);
Si = sum(dos_si(4:3+dos_si(2)))
C = sum(dos_c(4:3+dos_c(2)))
plot((freq-0.5)*dos_si(3)*UFREQ_IN_THZ,(dos_si(4:3+dos_si(2))/Si/dos_si(3)/UFREQ_IN_THZ+dos_c(4:3+dos_c(2))/C/dos_c(3)/UFREQ_IN_THZ)/2.);
title ('Total DOS of perfect SiC crystal from exact phonon dispersion calculation');
xlabel ('THz');
axis([ 0 35 -0. 0.3]);
print total.ps;
ss = input('next... ');

clf;
load dos_si.out;
freq = 1:dos_si(2);
Si = sum(dos_si(4:3+dos_si(2)))
plot((freq-0.5)*dos_si(3)*UFREQ_IN_THZ,dos_si(4:3+dos_si(2))/Si/dos_si(3)/UFREQ_IN_THZ);
title ('LDOS of Si in perfect crystal from exact phonon dispersion calculation');
xlabel ('THz');
axis([ 0 35 -0. 0.175]);
print exact_per_si.ps;
ss = input('next... ');

clf;
load dos_si_per.out;
plot(dos_si_per(:,1),dos_si_per(:,2),'o');
hold on;
plot(dos_si_per(:,1),dos_si_per(:,2),'-.');
axis([ 0 35 -0. 0.175]);
title ('Our method: 4096 particles, base frequency = 0.125 THz, 25 supercell k-points');
xlabel ('THz');
print our_per_si.ps;
ss = input('next... ');

clf;
load dos_c.out;
freq = 1:dos_c(2);
C = sum(dos_c(4:3+dos_c(2)))
plot((freq-0.5)*dos_c(3)*UFREQ_IN_THZ,dos_c(4:3+dos_c(2))/C/dos_c(3)/UFREQ_IN_THZ);
title ('LDOS of C in perfect crystal from exact phonon dispersion calculation');
axis([ 0 35 -0. 0.6]);
xlabel ('THz');
print exact_per_c.ps;
ss = input('next... ');

clf;
load dos_c_per.out;
plot(dos_c_per(:,1),dos_c_per(:,2),'o');
hold on;
plot(dos_c_per(:,1),dos_c_per(:,2),'-.');
axis([ 0 35 -0. 0.6]);
title ('Our method: 4096 particles, base frequency = 0.1 THz, 30 supercell k-points');
xlabel ('THz');
print our_per_c.ps;
ss = input('next... ');

clf;
load dos_si_def.out;
plot(dos_si_def(:,1),dos_si_def(:,2));
axis([ 0 40 -0. 0.175]);
title ('Change of Si LDOS in (111) direction after putting in an antisite-pair defect');
xlabel ('THz');
print our_def_si.ps;
ss = input('next... ');

clf;
load dos_c_def.out;
plot(dos_c_def(:,1),dos_c_def(:,2));
axis([ 0 40 -0. 0.6]);
title ('Change of C LDOS in (111) direction after putting in an antisite-pair defect');
xlabel ('THz');
print our_def_c.ps;
ss = input('next... ');


  