
%fn stands for the filename which has two columns.
%First column is KE from simulation
%Second column is corresponding total energies from simulation.

function energy_drift(fn)
A = textread(fn);
KE = A(:,3);
TE = A(:,4);

KE_S = size(KE);

ke_sum = 0;

for i=1:KE_S(1),
    ke_sum = ke_sum + KE(i);
end
ke_sum = ke_sum/KE_S(1);

f_te = TE(1)
te_size = size(TE);

NTE = TE - f_te;

indx = 0;
drift = [];
for j=1:te_size(1),
    drift(j) = NTE(j)*100/ke_sum;
end

plot(drift);
ke_sum
NTE(1000)
     TE(1000)
f_te
drift(1000)
