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

# ***** read u

utime = np.genfromtxt("utime_backstep.dat", dtype=None)

u1=utime[:,0] #x=3.3, y=0.037
u2=utime[:,1] #x=8.7, y=0.037
u3=utime[:,2] #x=18.7, y=0.037
u4=utime[:,3] #x=3.3, y=0.85
u5=utime[:,4] #x=8.7, y=0.85
u6=utime[:,5] #x=18.7, y=0.85
u7=utime[:,6] #x=-3.94, y=1.017
v1=utime[:,7] #x=3.3, y=0.85
v2=utime[:,8] #x=8.7, y=0.85
v3=utime[:,9] #x=18.7, y=0.85
t=utime[:,10]  # time

n=len(u1)

umean=np.mean(u3);

dt=8.177E-04;
t_tot=dt*len(u1)
t = np.linspace(0,t_tot,len(u1))

# subtract the mean
u3_fluct=u3-umean

# compute var=rms**2. This is the energy of the signal in time  (physical space)
u2=np.var(u3_fluct)

fig1 = plt.figure("Figure 3")
nblock = 1024
overlap = 128
win = hann(nblock, True)
fs=1

f, Pxx = welch(u3_fluct, fs, window=win, noverlap=overlap, nfft=nblock, return_onesided=True)

# This is the energy of the signal in time  (physical space).Should be eual to u2
energy_wave=np.sum(Pxx)/len(Pxx)

dt=t[1]-t[0]
f_phys=1./dt*np.linspace(0.,1.,len(Pxx))
plt.loglog(f_phys, Pxx, 'r-')
plt.xlabel("f [1/s]")
plt.ylabel("E(u) [m^2/s]")
plt.axis([1.e0,100,1e-4,10])
plt.title('spectrum of u')

# add line with -5/3 slope
xx=[1,10]
y=np.zeros(2)
ynoll=1
y[0]=ynoll
y[1]=y[0]*(xx[1]/xx[0])**(-5/3)
plt.loglog(xx,y,'k--',linewidth=4)

plt.savefig('u_spectum_python.eps')