% Exact calculation of K for reaction N2 + 3 H2 = 2 NH3
hf298_h2 = 0; gf298_h2 = 0; sf298_h2 = 0;
hf298_n2 = 0; gf298_n2 = 0; sf298_n2 = 0;
hf298_nh3 = -46110; gf298_nh3 = -16410; sf298_nh3 = (hf298_nh3-gf298_nh3)/298.15;
% Cp-data from Reid: cp(T) = cpa + cpb*T + cpc*T^2 + cpd*T^3
cpa_h2 = 2.714e1; cpb_h2 = 9.274e-3; cpc_h2 = -1.381e-5; cpd_h2 = 7.645e-9;
cpa_n2 = 3.115e1; cpb_n2 =-1.357e-2; cpc_n2 =  2.680e-5; cpd_n2 =-1.168e-8;
cpa_nh3= 2.731e1; cpb_nh3= 2.383e-2; cpc_nh3=  1.707e-5; cpd_nh3=-1.185e-8;
% Integrate, see (A.39) in 1. utgave
T=700; T0=298.15; 
hfT_h2 = hf298_h2 + cpa_h2*(T-T0) + cpb_h2*(T^2-T0^2)/2 + cpc_h2*(T^3-T0^3)/3 + cpd_h2*(T^4-T0^4)/4;
hfT_n2 = hf298_n2 + cpa_n2*(T-T0) + cpb_n2*(T^2-T0^2)/2 + cpc_n2*(T^3-T0^3)/3 + cpd_n2*(T^4-T0^4)/4;
hfT_nh3= hf298_nh3+ cpa_nh3*(T-T0)+ cpb_nh3*(T^2-T0^2)/2+ cpc_nh3*(T^3-T0^3)/3+ cpd_nh3*(T^4-T0^4)/4;
sfT_h2 = sf298_h2 + cpa_h2*log(T/T0) + cpb_h2*(T-T0) + cpc_h2*(T^2-T0^2)/2 + cpd_h2*(T^3-T0^3)/3;
sfT_n2 = sf298_n2 + cpa_n2*log(T/T0) + cpb_n2*(T-T0) + cpc_n2*(T^2-T0^2)/2 + cpd_n2*(T^3-T0^3)/3;
sfT_nh3= sf298_nh3+ cpa_nh3*log(T/T0)+ cpb_nh3*(T-T0)+ cpc_nh3*(T^2-T0^2)/2+ cpd_nh3*(T^3-T0^3)/3;
% reaction N2 + 3 H2 = 2 NH3
dhrT = 2*hfT_nh3 - hfT_n2 - 3*hfT_h2
dsrT = 2*sfT_nh3 - sfT_n2 - 3*sfT_h2
dgrT = dhrT - T*dsrT
K = exp(-dgrT/(8.3145*T))
