%uses the values calculated by base

nbins=150; %number of output values for the plot

%w are the results from the base script, which determines the dispersion
%relation for the fcc crystal
wn=w./sqrt(2); %factor to have values like those in debye
%there is an error in the dispersionrelation calculation in the function
%that we took from an old student project

wmin = min(wn(:)); %determines the smallest value
wmax = max(wn(:)); %determines the greatest value

wneu = linspace(wmin,wmax,nbins); %generates a vector from wmin to wmax with 150 elements

%DOS=sum(delta(w-w(i)))
binned = hist(wn(:),wneu);
%bins the elements in wn(:) in wneu equally spaced containers and returns
%the number of elements in each container as a row vector

dos = 4./(wmax-wmin).*nbins./(a*10^(-10)).^3./kvaln.*binned./10^6;
%* 10^6......m^3 => cm^3
%* kvaln.....number of considered k-vectors
%* a.........scaling to volume [angström]
%* (wmax-wmin)/nbins...dx (width of containers)
%* 4.........number of atoms in the fcc crystal
%            8 corner atoms *1/8 + 6 area atoms *1/2 = 4

plot(wneu./10^(12),smooth(dos)./10^9); %plots the dos
xlabel('\omega / 10^{12} s^{-1}')

%calculation of the standard values
%because we determinend the dispersionrelation with 
%a=4.49, C=200, m = 0.344*10^-24;
%so we have convert this values to a=1, C=1, m=1*10^-24
%perhaps it would be easier/better to calculate this values again - but it takes a
%lot of time
einhw=zeros(nbins,2);
einhw(:,1) = wneu./10^(12)*sqrt(m*10^24./C);
einhw(:,2) = smooth(dos).*sqrt(C./m./10^24)*a^3;
%plot(einhw(:,1), einhw(:,2))
