close all
clear all
% ***** read u
load u_time_interior.dat
vi1_j1=u_time_interior(:,1);
vi2_j1=u_time_interior(:,2);
vi3_j1=u_time_interior(:,3);
vi4_j1=u_time_interior(:,4);

vi1_j2=u_time_interior(:,5);
vi2_j2=u_time_interior(:,6);
vi3_j2=u_time_interior(:,7);
vi4_j2=u_time_interior(:,8);

%     write(67,25)phi(17,10,k,v),phi(25,10,k,v),
%    .            phi(35,10,k,v),phi(45,10,k,v),
%    .            phi(17,30,k,v),phi(25,30,k,v),
%    .            phi(35,30,k,v),phi(45,30,k,v)



% time step=0.0049
dt= 6.250E-04;
n=length(vi1_j1);
% compute time array
t=dt:dt:n*dt;

%%%%%%%%%%%%%%%% plotting section %%%%%%%%%%%%%%%%%%%%%%%%%%
% plot v
figure(1)
plot(t,vi1_j1)
hold
plot(t,vi3_j1,'r--')
xlabel('t')
ylabel('v')
handle=gca;
set(handle,'fontsize',[20])
print v_time_pans.ps -deps


%%%%%%%%%%%%%%%%%% auto correlation %%%%%%%%%%%%%%
figure(2)

a=0.954; 
b=(1-a^2)^0.5;

% a=exp(-dt/T)
T=-dt/log(a);


% take the auto corr over 'imax' time steps
imax=500
for i=1:imax
   t1=t(i);
   two_a(i)=exp(-t1/T);
end


if exist('OCTAVE_VERSION');
  result = true;
else
  result = false
end
if result
% octave
   two_vvt=autocor(vi3_j1,imax);
else
% matlab
   two_vvt=autocorr(vi3_j1,imax);
end
% normalize auto corr
two_vvt=two_vvt/two_vvt(1);


% compute integral length scale
T_int=0;
for ii=1:imax
   if two_vvt(ii) < 1e-6
      break
   end
   T_int=T_int+two_vvt(ii)*dt;
end

disp(['Prescribed time scale at interface= ',num2str(T)])
disp([' '])
disp(['Integral time scale at node(i3,j1)= ',num2str(T_int)])


plot(t(1:imax),two_vvt(1:imax),'linewidth',2)
xlabel('t')
ylabel('two')
axis([0 0.03 -0.1 1])
handle=gca;
set(handle,'fontsize',[20])
print auto_corr.ps -deps

