N2, O2 and H2 density

o2h2n2density.m

clear;
clc;

#O2 and H2 density
O2Tcrit = 154.7;
H2Tcrit = 33.2;
N2crit = 126.2;
temp1 = linspace(0,O2Tcrit,75);
temp2 = linspace(O2Tcrit,500,100);
temp3 = linspace(0,H2Tcrit,15);
temp4 = linspace(H2Tcrit,500,200);
temp5 = linspace(0,N2crit,50);
temp6 = linspace(N2crit,500,120);

O2Dens1 = (1434.338998096669 + 30827.66466562366 ./ (temp1 - 186.3966881148979));
O2Dens2 = (44.24461555143480 + 7784.502442355128 ./ (temp2 - 136.05498694800465));

H2Dens1 = (136.4894046680936 + 3242.617524782929 ./ (temp3 - 67.46034096292647));
H2Dens2 = (0.741833633973125 + 642.4040759445162 ./ (temp4 - 17.5701803944558));

N2Dens1 = (734.3287921946625 + 9878.83146453045 ./ (temp5 - 146.65628914669438));
N2Dens2 = (17.85307222754917 + 4418.262547908501 ./ (temp6 - 107.28230096889227));

O2DensOld = 0.56 * [temp1,temp2].^2 - 134.0 .*[temp1,temp2] + 6900.0 + 1141.0;
H2DensOld = 0.03333 * [temp3,temp4].^2 - 4.3333 .*[temp3,temp4] + 73.3333 + 70.0;
N2DensOld = 807.0.*ones(1,length([[temp5],[temp6]]));


figure(1)
axis([0 500 2 1400]);

grid on;
grid minor;
hold on;
##
semilogy([temp1,temp2],[O2Dens1,O2Dens2],'bo-');
semilogy([temp1,temp2],O2DensOld,'ro-');
##
semilogy([temp3,temp4],[H2Dens1,H2Dens2],'b+-');
semilogy([temp3,temp4],H2DensOld,'r+-');
##
semilogy([temp5,temp6],[N2Dens1,N2Dens2],'bx-');
semilogy([temp5,temp6],N2DensOld,'rx-');
hold off;

xlabel("Temperature [K]");
ylabel("Density [g/l]");
legend("Hume O_2","Legacy SPSDK O_2","Hume H_2","Legacy SPSDK H_2","Hume N_2","Legacy SPSDK N_2");
title("Liquid Density Model Comparison");

N2Dens1(50)-N2Dens2(1)
N2Dens2(120)