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 = 14000
# Air
Cp2 = 1005

# 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

h = np.zeros(nx)
h = Cp*T

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])
    Cp1_c    = Cp1
    Cp2_c    = Cp2
    T_c     = 0.5*(T[i]+T[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* (Cp1_c*T_c*gradY1 + Cp2_c*T_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 constant with temperature')
plt.grid()
plt.savefig('CpiConstWithTemperature.png')