# --------------- Import modules ---------------- #
import math
import numpy as np
import matplotlib.pyplot as plt
#plt.rc('tex', usetex=True)
plt.rc('font', family='serif')
import os

cwd = os.getcwd()
# ------------- Material constants -------------- #
rho = 998; g = 9.81;
patm = 101325; dT = 2.5e-4;
pScale = 1e-3;

# ------------- Read surface data --------------- #
Q = np.loadtxt('postProcessing/flowRateInlet/0/surfaceFieldValue.dat')
time = Q[:,0]
Q = np.abs(Q[:,1]*1e3)

# ------------- Read probe data ----------------- #
pProbes = np.loadtxt('postProcessing/probes/0/p')*rho*pScale;

# --------- Read Numerical Wang data ------------ #
folder= cwd + '/validationData/';
wangData = np.genfromtxt(folder + 'wangEtAl_numericalData.csv',delimiter=',');
wangData[0,0]=0.00005; # error reading first value
wangData[:,0] = wangData[:,0] + 3 - wangData[0,0];
wangData[:,1] = wangData[:,1]*1e3;
wangData[:,2] = wangData[:,2]*rho*10*pScale;
#print(wangData)
expData = np.loadtxt(folder + 'petitEtAl_experimentalData.txt',comments='%');
expData[:,1] = expData[:,1]*pScale;

# ----------------- Plot Data ------------------- #

if not os.path.exists('Figures'):
    os.makedirs('Figures')

# plotting the points
plt.figure(1)
plt.plot(time, Q, 'k-', label = 'headLossPressure')
plt.plot(wangData[:,0], wangData[:,1],'r--', label = 'Wang et al. (2017)')
plt.xlabel('Time (s)')
plt.ylabel('Flow rate (l/s)')
plt.legend(loc="lower right") 
plt.xlim([0, 30]); plt.ylim([-20, 60])
plt.grid()
plt.savefig('Figures/Q.pdf')


plt.figure(2)
plt.plot(time, pProbes[:,1],'k-', label = 'headLossPressure')
plt.plot(wangData[:,0], wangData[:,2],'r--', label = 'Wang et al. (2017)')
plt.plot(expData[:,0], expData[:,1],'c-.', label = 'Petit et al. (2015)')
plt.xlabel('Time (s)')
plt.ylabel('Pressure (kPa)')
plt.legend(loc="upper right") 
plt.xlim([0, 30]);
plt.ylim([0, 50]);
plt.grid()
plt.savefig('Figures/p.pdf')


plt.show()


