import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch, hanning
plt.rcParams.update({'font.size': 22})

# ***** read u

data = np.genfromtxt("u_time_interior.dat",comments="%")

vi1_j1=data[:,0]   
vi2_j1=data[:,1]   
vi3_j1=data[:,2]   
vi4_j1=data[:,3]   

vi1_j2=data[:,4]   
vi2_j2=data[:,4]   
vi3_j2=data[:,4]   
vi4_j2=data[:,4]  


#     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
dt= 6.250E-04
t_tot=dt*len(vi1_j1)
t = np.linspace(0,t_tot,len(vi1_j1))


#%%%%%%%%%%%%%%% plotting section %%%%%%%%%%%%%%%%%%%%%%%%%%
# plot v
fig1 = plt.figure("Figure 1")
plt.clf() #clear the figure
plt.plot(t,vi1_j1)
plt.plot(t,vi3_j1,'r--')
plt.xlabel('t')
plt.ylabel('v')
plt.savefig('v_time_pans_python.ps.eps')


#%%%%%%%%%%%%%%%%% auto correlation %%%%%%%%%%%%%%
fig2 = plt.figure("Figure 2")
plt.clf() #clear the figure

a=0.954 
b=(1-a**2)**0.5

# a=exp(-dt/T)
T=-dt/np.log(a)


# take the auto corr over 'imax' time steps
imax=500
two_a=np.zeros(imax)
for i in range(0,imax-1):
   t1=t[i]
   two_a[i]=np.exp(-t1/T)

two=np.correlate(vi1_j1,vi3_j1,'full')

# find max
nmax=np.argmax(two)
# and its value
two_max=np.max(two)

# two_max is symmwetric. Pick the right half and normalize
two_sym_norm=two[nmax:]/two_max


# compute integral length scale
T_int=0
for ii in range(0,imax):
   if two_a[ii] < 1e-6:
      break
   T_int=T_int+two_a[ii]*dt

print('Prescribed time scale at interface= ',T)
print('Integral time scale at node(i3,j1)= ',T_int)


plt.plot(t[0:imax],two_a[0:imax],'b-')
plt.xlabel('t')
plt.ylabel('two')
plt.axis([0,0.03,-0.1,1])
plt.savefig('auto_corr_python.ps.eps')