function main
[v,h,s,cp,u,cv,w] = wxpt1(50,500)
function [v,h,s,cp,u,cv,w] = wxpt1(p,t)
e=[0 0 0 0 0 0 0 0 1 1 1 1 1 1 2 2 2 2 2 3 3 3 4 4 4 5 8 8 21 23 29 30 31 32];
n=[-2 -1 0 1 2 3 4 5 -9 -7 -1 0 1 3 -3 0 1 3 17 -4 0 6 -5 -2 10 -8 -11 -6 -29 -31 -38 -39 -40 -41];
m=[0.14632971213167 -0.84548187169114 -0.37563603672040*10 0.33855169168385*10 -0.95791963387872 0.15772038513228 -0.16616417199501*10^-1 0.81214629983568*10^-3 0.28319080123804*10^-3 -0.60706301565874*10^-3 -0.18990068218419*10^-1 -0.32529748770505*10^-1 -0.21841717175414*10^-1 -0.52838357969930*10^-4 -0.47184321073267*10^-3 -0.30001780793026*10^-3 0.47661393906987*10^-4 -0.44141845330846*10^-5 -0.72694996297594*10^-15 -0.31679644845054*10^-4 -0.28270797985312*10^-5 -0.85205128120103*10^-9 -0.22425281908000*10^-5 -0.65171222895601*10^-6 -0.14341729937924*10^-12 -0.40516996860117*10^-6 -0.12734301741641*10^-8 -0.17424871230634*10^-9 -0.68762131295531*10^-18 0.14478307828521*10^-19 0.26335781662795*10^-22 -0.11947622640071*10^-22 0.18228094581404*10^-23 -0.93537087292458*10^-25];
c=0.461526;
p0=16.53;
t0=1386.0;
pi=p/p0;
tao=t0/t;
Epi=7.1-pi;
Etao=tao-1.222;
r=0;
rpi=0;
rtao=0;
rpi2=0;
rtao2=0;
rpitao=0;
for i=1:1:34
dizeng=e(i)*(Epi.^n(i))*(Etao.^m(i));
r=r+dizeng;
rpi=rpi-dizeng*n(i)/Epi;
rtao=rtao+dizeng*m(i)/Etao;
rpi2=rpi2+dizeng*n(i)*(n(i)-1)/Epi.^2;
rtao2=rtao2+dizeng*m(i)*(m(i)-1)/Etao.^2;
rpitao=rpitao-dizeng*m(i)*n(i)/Etao/Epi;
end
g=r*c*t;%吉布斯自由能
v=rpi*pi*(c*1000)*t/p;%比容
h=rtao*tao*c*t;%比焓
s=(h-g)/t;%比熵
cp=-tao*tao*rtao2*c;%定压比热
u = c*t*(tao*rtao-pi*rpi);%内能
cv=cp+c*(rpi-tao*rpitao).^2/rpi2;%定容比热
w=((c*t*rpi.^2)/((rpi-tao*rpitao).^2/(tao.^2*rtao2)-rpi2)).^0.5;%声速