% file autocorr.m
%
% arguments: (data,skip,tmax)
% data: column of time history of MC estimates
% skip: # of initial data to be disregarded
% tmax: autocorrelation constructed up to tmax
%
% return: [av,err,tauint]
% av:     average
% err:    error (corrected for autocorrelation)
% tauint: integrated autocorr. time
%
% summation window for tauint is asked interactively
%
function [av,err,tauint] = autocorr(data,skip,tmax)

if size(data,2) ~= 1
  disp('data not column, return');
  return;
end;

n=size(data,1);
d=data(skip+1:n);
av=mean(d);
d=d-av;
n=n-skip;
for t=0:tmax,
  gamma(t+1)=(d(1:n-t)'*d(1+t:n))/(n-t);
end
var=gamma(1);
gamma=gamma/var;
plot(0:tmax,gamma,'x');

window=input(' How far sum for tau_int?  ');
tauint=sum(gamma(2:window+1))+0.5;
err=sqrt(var*2*tauint/n);
