# Start writing code here...
import matplotlib.pyplot as plt
import scipy.optimize as opt
import numpy as np
# User specified constants
T0 = 298.15 # Standard Temperature
R = 8.314 # Ideal gas constant
T = 873 # Ending Temperature
P1 = 1.03215
# Free energies of formation, Appendix 4, Sandler, 2nd Edition
delGrxn = 3.8e3
# Calculation of equilibrium constant at T0
KaT0 = np.exp(-delGrxn/(R*T0))
delHstd = -49e3
#deltaCp = np.array([26.143, 10.415e-2, -9.209e-5, 27.791e-9])
#const = deltaCp[0] * T0 + (deltaCp[1] / 2)*T0**2 + (deltaCp[2] / 3)*T0**3 + (deltaCp[3] / 4) *T0**4
KaT = KaT0 * np.exp(-(delHstd/ R) * ((1 / T)-(1 / T0)))
def fun(x):
return (x**2)/(1-x**2) - KaT
out= opt.root(fun, 0.1)
x = float(abs(out.x)) # out.x is an array and float() simplifies output below
print(x)
KaT
0.0006928451563432711