import scipy.io as sio
import numpy as np
import sys
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 26})

# Channel flow
viscos=1/8000
data=np.loadtxt("y_u_k_eps_8000.dat")

j=data[:,0]
y=data[:,1]
u=data[:,2]
k=data[:,3]
eps=data[:,4]

# Q2/1 ##################################################################
# compute u_tau
dudy=np.gradient(u,y)
u_tau=(viscos*dudy[0])**0.5
print('\nu_tau=',u_tau)
yplus=y*u_tau/viscos

# Q2/2 ##################################################################
# 
dudy=np.gradient(u,y)

# the velocity graient is largest at the wall and smallest at the centerline
# check it
print('\ndudy largest at j=', np.argmax(dudy))
print('\ndudy smallest at j=', np.argmin(abs(dudy)))
#************
# velocity profile plot
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.25,bottom=0.20)
plt.plot(u,y,'b-')
plt.title('Velocity profile')
plt.axis([0,30,0,2])
plt.xlabel('$V_1$') 
plt.ylabel('$x_2$') 
plt.savefig('velprof.eps')

# Q2/3 ##################################################################

# wall b.c. for eps
eps_wall_bottom=2*viscos*k[0]/y[0]**2
ymax=2
eps_wall_top=2*viscos*k[0]/(ymax-y[-1])**2

# they agree with the data (i.e. with eps[0] and eps[-1]))

# Q2/4 ##################################################################

# damping function fmu, see Eq. F.3
fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.25,bottom=0.20)

ueps=(eps*viscos)**0.25
# compute the disytance to the closest wall
dist=np.minimum(y,y[-1]-y)
ystar=ueps*dist/viscos
rt=k**2/eps/viscos
fmu=((1.-np.exp(-ystar/14.))**2)*(1.+5./rt**0.75*np.exp(-(rt/200.)**2))
plt.plot(fmu,y,'b-')
plt.title('damping function')
plt.axis([0,1.1,0,0.05])
plt.xlabel('$f_\mu$') 
plt.ylabel('$x_2$') 
plt.savefig('fmu.eps')


# Q2/5 ##################################################################

fig1,ax1 = plt.subplots()
plt.subplots_adjust(left=0.2,bottom=0.25)
dkdy=np.gradient(k,y)
d2kdy2=np.gradient(dkdy,y)
# viscous diffusion term
diff_term=viscos*d2kdy2
plt.plot(diff_term,yplus,'b-')
plt.title('diffusion term')
plt.axis([-1000,1000,0,100])
plt.xlabel(r'$\frac{\partial^2 k}{\partial x_2^2}$') 
plt.ylabel('$x_2$') 
plt.savefig('diff-term.eps')