function cp = heat_capacity(x)
A = co2_sw(x);
cp = (-A.g(1)-x(1)*A.H(1,1)+x(3)*A.H(3,1))/x(3);
