import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# Domain size 
x_start, x_end = 0,1

# number of nodes
nx  = 101

# values of x
x = np.linspace(x_start, x_end, nx)

# Temperature
T = np.linspace(275,2000,nx)

# Specific heats 
# Hydrogen
Cp1 = np.zeros(nx)
Cp1 = 1000*(-1.64044e-18*pow(T,6) +1.22802e-14*pow(T,5) -3.70152e-11*pow(T,4) +5.6503e-8*pow(T,3) -4.4297e-5*pow(T,2) +0.0172878*T +11.8523)
# Air
Cp2 = np.zeros(nx)
Cp2 = 1000*(7.36435e-22*pow(T,6) -8.89851e-17*pow(T,5) +5.73158e-13*pow(T,4) -1.44058e-9*pow(T,3) +1.64424e-6*pow(T,2) -0.000626846*T +1.07935)

# Mass fractions
Y1 = np.linspace(0,1,nx)
Y2 = 1-Y1

# Thermal conductivity
kappa1 = 0.4
kappa2 = 5e-2

kappa = np.zeros(nx)
kappa = kappa1*Y1+kappa2*Y2

Cp = np.zeros(nx)
Cp = Cp1*Y1+Cp2*Y2

h1 = np.zeros(nx)
for i in range(0,nx):
    Te      = T[i]
    upLim   = 1000*(-1.64044e-18*pow(Te,7)/7 +1.22802e-14*pow(Te,6)/6 -3.70152e-11*pow(Te,5)/5 +5.6503e-8*pow(Te,4)/4 -4.4297e-5*pow(Te,3)/3 +0.0172878*pow(Te,2)/2 +11.8523*Te)
    Te      = T[0]
    lowLim  = 1000*(-1.64044e-18*pow(Te,7)/7 +1.22802e-14*pow(Te,6)/6 -3.70152e-11*pow(Te,5)/5 +5.6503e-8*pow(Te,4)/4 -4.4297e-5*pow(Te,3)/3 +0.0172878*pow(Te,2)/2 +11.8523*Te)
    h1[i]   = Cp1[0]*T[0]+upLim-lowLim

h2 = np.zeros(nx)
for i in range(0,nx):
    Te      = T[i]
    upLim   = 1000*(7.36435e-22*pow(Te,7)/7 -8.89851e-17*pow(Te,6)/6 +5.73158e-13*pow(Te,5)/5 -1.44058e-9*pow(Te,4)/4 +1.64424e-6*pow(Te,3)/3 -0.000626846*pow(Te,2)/2 +1.07935*Te)
    Te      = T[0]
    lowLim  = 1000*(7.36435e-22*pow(Te,7)/7 -8.89851e-17*pow(Te,6)/6 +5.73158e-13*pow(Te,5)/5 -1.44058e-9*pow(Te,4)/4 +1.64424e-6*pow(Te,3)/3 -0.000626846*pow(Te,2)/2 +1.07935*Te)
    h2[i]   = Cp2[0]*T[0]+upLim-lowLim

h = np.zeros(nx)
h = h1*Y1+h2*Y2

alpha = np.zeros(nx)
alpha = kappa/Cp

# Values at cell centres
Exact       = np.zeros(nx-1)
OpenFoam    = np.zeros(nx-1)
x_c         = np.zeros(nx-1)

for i in range(0,nx-1):
    kappa_c = 0.5*(kappa[i]+kappa[i+1])
    alpha_c = 0.5*(alpha[i]+alpha[i+1])
    Cp_c    = 0.5*(Cp[i]+Cp[i+1])
    h1_c    = 0.5*(h1[i]+h1[i+1])
    h2_c    = 0.5*(h2[i]+h2[i+1])
    
    gradT   = (T[i+1]-T[i])/(x[i+1]-x[i])
    gradY1   = (Y1[i+1]-Y1[i])/(x[i+1]-x[i])
    gradY2   = (Y2[i+1]-Y2[i])/(x[i+1]-x[i])
    gradH   = (h[i+1]-h[i])/(x[i+1]-x[i])

    x_c[i]      = x[i]+0.5*(x[i+1]-x[i])
    Exact[i]    = kappa_c*gradT + alpha_c* (h1_c*gradY1 + h2_c*gradY2)
    OpenFoam[i] = alpha_c*gradH 


plt.plot(x_c,Exact,'b-',label="Exact")
plt.plot(x_c,OpenFoam,'r--',label="OpenFOAM")
plt.legend()
plt.xlabel('X [m]')
plt.ylabel('Enthalpy from Heat conduction and Specie diffusion')
plt.title('Cpi variation with temperature')
plt.grid()
plt.savefig('CpiVaryWithTemperature.png')