function Hrx = HtRxnP(T,P)
% heat of reaction coefficients from spreadsheet fit
% basis: Hysys molar enthalpy predictions at 150 atm
% includes pressure effect using Pitzer correlation
R = 8.31446;  % kJ/kgmol/K
aa = -1.9314e5;
bb = 4.8403e5;
cc = -9.944e5;
dd = 8.8054e5;
ee = -2.9078e5;
Tk = T + 273.15; % K
Tk1 = Tk/1000;
Pk = P *101.325; % kPa
HrxT = aa + bb*Tk1 + cc*Tk1^2 + dd*Tk1^3 + ee*Tk1^4;
Tc = [ 405.6 126.2 33.2 ];
Pc = [ 112.5 33.5 12.8 ];
Omega = [ 0.25 0.04 0.0 ];
Tr = Tk ./ Tc;
Pck = Pc * 101.325;
B0 = 0.1445 - 0.33 ./ Tr - 0.1385 ./ Tr.^2 - 0.0121 ./ Tr.^3;
B1 = 0.073 + 0.46 ./ Tr - 0.5 ./ Tr.^2 - 0.097 ./ Tr.^3 - 0.0073 ./ Tr.^8;
dB0 = 0.33 ./ Tr.^2 + 0.277 ./ Tr.^3 + 0.0363 ./ Tr.^4;
dB1 = -0.46 ./ Tr.^2 + 1 ./ Tr.^3 + 0.291 ./ Tr.^4 + 0.0584 ./ Tr.^9;
H = R * Tk * ( 1-Pk) ./ Pck .* ( (dB0 - B0 ./ Tr) + Omega .* (dB1 - B1 ./ Tr));
Hrx = HrxT + [ 2 -1 -3 ]*H';


