% //***********************************************************************************************************************
% // ACcalculations2.cpp
% //
% // Combines the features of breakup2.cpp and ACcalculations.cpp into a single program.
% // This program takes heat current data and calculates the
% // heat-current auto-correlation function which can be used to determine thermal conductivity.
% //
% // Written by Alan McGaughey
% //
% // initial development fall 2000
% //
% // last updated February 8, 2004
% //***********************************************************************************************************************

number = 1;        % Number of output HC files required.
directions = 1;    % Number of different directions in the crystal.

% Opening input and output files.
inHeatCurrent = fopen('heatcurrent.txt', 'r');
outACData = fopen('ACDataOut.txt', 'w');
outIntData = fopen('IntOut.txt', 'w');

%Declaring some variables.
step = 1;
file_length = fscanf(inHeatCurrent, '%d', 1);    % Number of data lines in the input file. 
length = 300;                                    % Number of points in AC function.
ensemble = (file_length-length)/step + 1;        % Number of points in ensemble average.
Q = zeros(file_length,3);                        % Heat current array.
AC = zeros(length,3);                            % Autocorrelation array.
input_length = file_length;                      % Number of line in input file, data in.
a = file_length;                                 % Number of lines per output file.
b = 0;                                           % Overlap - Number of output files required is (input_length-b)/(a-b) (effect of step not included)
step2 = 1;
delt = 0.002;
gap = 10;                                        % Must be equal to tramp in Main MATLAB file.
ind = 1;                                         % Number of independant data sets.
integral = zeros(length,directions);
data = zeros(4,1);
timeStep = zeros(file_length,1);

if (directions == 2) quartz = 1;
else quartz = 0; end;

% Outputting the number of line in the output file.
fprintf(outACData, '%d\n', length);
fprintf(outIntData, '%d\n', length);

% This program runs for only 1 input data file. No break ups.
T = 10; 
fprintf('Begin to input heatcurrent data file.\n');
 
% Read in heat current data.
for i = 1:file_length
  data = fscanf(inHeatCurrent, '%d %e %e %e %e', 5);
  timeStep(i,1) = data(1,1);
  Q(i,1) = data(2,1);
  if (quartz == 1) Q(i,2) = data(3,1);
  else Q(i,2) = data(3,1); end;
  Q(i,3) = data(4,1);
end; 
  
fprintf('Begin to calculate correlation.\n');
for k = 1:ensemble
  if (mod(k,100)==0) fprintf('ensemble: %d', k); end;
  for i = 1:length
    for j = 1:3 AC(i,j) = AC(i,j)+Q(k*step,j)*Q(i+k*step-1,j); end;
  end;
end;

% Average over similar directions and time origins
if (quartz == 1)
  for i = 1:length
    AC(i,1) = AC(i,2)/(number*ensemble*ind);
    AC(i,2) = (AC(i,1)+AC(i,3))/(2.*number*ensemble*ind);
  end;
else
  for i = 1:length
    AC(i,1) = (AC(i,1)+AC(i,2)+AC(i,3))/(3.*number*ensemble*ind);
  end;
end;

fprintf('Begin to calculate integration.\n');
% Do integration.
for i = 1:directions
  integral(1,i) = 0.;
  for j = 2:length integral(j,i) = integral(j-1,i) + delt*gap*AC(j,i); end;
end;

if (quartz == 1)
  for i = 1:length fprintf(outACData, '%e\t %e\t %e\t %e\n', AC(i,1), integral(i,1), AC(i,2), integral(i,2)); end;
else
  for i = 1:length
    fprintf(outACData, '%d\t %e\n', timeStep(i,1), AC(i,1));
    fprintf(outIntData, '%d\t %e\n', timeStep(i,1), integral(i,1));
  end;
end;

fprintf('Finished succesfully\n');

% Closing all files.
fclose(inHeatCurrent);
fclose(outACData);
fclose(outIntData);