!pip3 install --upgrade git+https://github.com/ptgodart/nasaPoly.git
%matplotlib inline

import nasaPoly

mgwhatever = nasaPoly.Species('Mg2SiO4(L)')
print(mgwhatever)

import nasaPoly # for thermodynamic functions
from matplotlib import pyplot as plt # for plotting
import numpy as np # for general math stuff

# Enter solutions below:
# Create class instances for each molecule:
H2 = nasaPoly.Species('H2')
N2 = nasaPoly.Species('N2')
CO2 = nasaPoly.Species('CO2')

# Define temperature range and other constants
T_range = np.linspace(200, 1000, 250) # K
R = 8.314 # J/mol-K

# Compute cv and cp for each molecule
cp_H2_range = [H2.cp_0(T) for T in T_range] # List comprehension
cv_H2_range = [H2.cp_0(T) - R for T in T_range]
cp_N2_range = [N2.cp_0(T) for T in T_range]
cv_N2_range = [N2.cp_0(T) - R for T in T_range]
cp_CO2_range = [CO2.cp_0(T) for T in T_range]
cv_CO2_range = [CO2.cp_0(T) - R for T in T_range]

plt.plot(T_range, cp_H2_range, '-', color='tab:blue', label='$c_p^{H_2}$')
plt.plot(T_range, cv_H2_range, '--', color='tab:blue', label='$c_v^{H_2}$')
plt.plot(T_range, cp_N2_range, '-', color='tab:orange', label='$c_p^{N_2}$')
plt.plot(T_range, cv_N2_range, '--', color='tab:orange', label='$c_v^{N_2}$')
plt.plot(T_range, cp_CO2_range, '-', color='tab:green', label='$c_p^{CO_2}$')
plt.plot(T_range, cv_CO2_range, '--', color='tab:green', label='$c_v^{CO_2}$')
plt.xlabel('T [K]')
plt.ylabel('Specific Heat [J/mol-K]')
plt.legend()

plt.show()

from scipy.integrate import quad

# Enter solutions below:
# (hint: you can integrate functions using quad(func, lower_bound, upper_bound))

T_lower = [200, 800] # K
T_upper = [300, 900] # K

for T_l, T_u in zip(T_lower, T_upper):
  Delta_Q_H2_p = quad(H2.cp_0, T_l, T_u)[0]
  Delta_Q_H2_v = quad(lambda T: H2.cp_0(T) - R, T_l, T_u)[0]
  Delta_Q_N2_p = quad(N2.cp_0, T_l, T_u)[0]
  Delta_Q_N2_v = quad(lambda T: N2.cp_0(T) - R, T_l, T_u)[0]
  Delta_Q_CO2_p = quad(CO2.cp_0, T_l, T_u)[0]
  Delta_Q_CO2_v = quad(lambda T: CO2.cp_0(T) - R, T_l, T_u)[0]

  print(f'Temperature Range, {T_l}-{T_u} K:')
  print(f'H2: Delta_q at constant pressure = {Delta_Q_H2_p:0.3g} J/mol')
  print(f'H2: Delta_q at constant volume = {Delta_Q_H2_v:0.3g} J/mol')
  print(f'N2: Delta_q at constant pressure = {Delta_Q_N2_p:0.3g} J/mol')
  print(f'N2: Delta_q at constant volume = {Delta_Q_N2_v:0.3g} J/mol')
  print(f'CO2: Delta_q at constant pressure = {Delta_Q_CO2_p:0.3g} J/mol')
  print(f'CO2: Delta_q at constant volume = {Delta_Q_CO2_v:0.3g} J/mol')
  print('\n')

# Enter solutions below:
# *liquid* H2O:
H2O_l = nasaPoly.Species('H2O(L)')

T_H2O_l_range = np.linspace(5+273.15, 95+273.15, 250)
cp_H2O_l_range = [H2O_l.cp_0(T) for T in T_H2O_l_range]

plt.plot(T_H2O_l_range, cp_H2O_l_range, label='$c_p^{H_2O_{(l)}}$')
plt.xlabel('T [K]')
plt.ylabel('Specific Heat [J/mol-K]')
plt.legend()

plt.show()

# Enter solutions below:
# Same as part b, but just for constant pressure
T_lower_d = [5+273.15, 90+273.15] # K
T_upper_d = [10+273.15, 95+273.15] # K

for T_l, T_u in zip(T_lower_d, T_upper_d):
  Delta_Q_H2O_p = quad(H2O_l.cp_0, T_l, T_u)[0]

  print(f'{T_l - 273.15}-{T_u - 273.15} degC:')
  print(f'H2O: Delta_Q at constant pressure = {Delta_Q_H2O_p:0.4g} J/mol')
  print('\n')

# Enter solutions below:
# Create class instances for C and O2:
C = nasaPoly.Species('C(gr)') # for graphite
O2 = nasaPoly.Species('O2')

# Define delta h of reaction h_products(T_final) - h_reactants(T_initial)
# where T_final is the adiabatic flame temperature
def h_rxn_i(T_final, T_initial):
  return np.abs(CO2.h_0(T_final) - (C.h_0(T_initial) + O2.h_0(T_initial)))

# Let's plot to see what this looks like. Remember, to find the adiabatic flame
# tempereature we need to find where the h_reactants = h_products, or where
# h_rxn_i = 0:
T_i = 300 # K
T_sweep = np.linspace(200, 8000, 250)
plt.plot(T_sweep, [h_rxn_i(T, T_i) for T in T_sweep], label='h_rxn')
plt.plot([T_sweep[0], T_sweep[-1]], [0, 0], '--', color='black', label='h_rxn=0')
plt.xlabel('T [K]')
plt.ylabel('h_rxn [J/mol]')
plt.grid()
plt.legend()
plt.show()

from scipy.optimize import fminbound

T_final_i = fminbound(h_rxn_i, 200, 10000, args=(T_i,))
print(f'Adiabatic flame temperature in pure O2: {T_final_i:0.4g} K')

T_i = 300 # K
T_sweep = np.linspace(200, 8000, 250)
plt.plot(T_sweep, [h_rxn_i(T, T_i) for T in T_sweep], label='h_rxn')
plt.plot([T_sweep[0], T_sweep[-1]], [0, 0], '--', color='black', label='h_rxn=0')
plt.plot([T_final_i], [h_rxn_i(T_final_i, T_i)], 'o', color='tab:red')
plt.xlabel('T [K]')
plt.ylabel('h_rxn [J/mol]')
plt.grid()
plt.legend()
plt.show()

# Just as before:
def h_rxn_ii(T_final, T_initial):
  return np.abs((CO2.h_0(T_final) + 3.773*N2.h_0(T_final)) - (C.h_0(T_initial) + O2.h_0(T_initial) + 3.773*N2.h_0(T_initial)))

T_final_ii = fminbound(h_rxn_ii, 200, 5000, args=(T_i,))
print(f'Adiabatic flame temperature in air: {T_final_ii:0.4g} K')