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

ni=66
nj=49
nk=66

# read u1 & transform u1 to a 3D array
uvw = sio.loadmat('u1_hybrid.mat')
tt=uvw['u1_hybrid']
u3d= np.reshape(tt,(ni,nj,nk))

u3d=np.swapaxes(u3d,0,2) 

# x coordinate direction = index i, first index
# y coordinate direction = index j, second index
# z coordinate direction = index k, third index

# compute average u by averaging in x (i) and z (k) directions
umean=np.mean(u3d, axis=(0,2))


# face coordinates
yc = np.loadtxt("yc_hybrid.dat")

y=np.zeros(nj)

for j in range(1,nj-1):
   y[j]=0.5*(yc[j]+yc[j-1])

y[nj-1]=yc[nj-1]
yplus=y*8000

#************ U
fig1 = plt.figure("Figure 1")
plt.clf() #clear the figure

plt.semilogx(yplus,umean,'b-')
plt.ylabel("$y^+$")

B=5.2
kappa=0.4
# find y+=30
ulog=np.zeros(nj)
for j in range(1,nj-1):
   ulog[j]=np.log(yplus[j])/kappa+B
   if yplus[j] < 30 and yplus[j+1] > 30:
      jj=j

plt.plot(yplus[jj:nj-1],ulog[jj:nj-1],'+')

plt.axis([1, 8200, 0, 31])
plt.xlabel("$U^+$")
plt.savefig('u_log_python.eps')