求解湿球温度matlab小程序

%已知干球温度、相对湿度，求湿球温度
%该程序仅适用于常压、0-90℃范围内，通过迭代求得不是非常精确的湿球温度解
P=input('输入环境压力Pa：');%对应环境压力，Pa
Td=input('输入干球温度K：');%干球温度，K
RH=input('输入相对湿度（0-1）：');%相对湿度，0-1
Psd=10^(10.11564-1687.54/(Td+230.23-273.15));%干球温度对应的饱和蒸汽压力，Pa
d=0.622*RH*Psd/(P-RH*Psd);%干球温度条件下含湿量，kg/kg(dry air)
h=1.01*Td+(2500+1.84*Td)*d;%干球温度对应焓值，KJ/kg
Tw=298.15;%假定一个湿球温度，K
hw=h;%湿球温度焓值=干球温度焓值
i=1;
Twk=Tw;
j=1;
while i>10^-7
Tw=Twk;
dw=(hw-1.01*Tw)/(2500+1.84*Tw);%湿球温度条件下含湿量，kg/kg(dry air)
Psw=dw*P/(dw+0.622);%湿球温度对应的饱和蒸汽压力，Pa
Twk=42.92+1687.54/(5.11564-log10(Psw*10^(-5)));%迭代出的新湿球温度，K
i=abs(Twk-Tw)/Tw;
j=j+1;
end
fprintf('湿球温度是 %g K\n',Tw);
fprintf('干球含湿量是 %g kg/kg(dry air)\n',d);
fprintf('湿球含湿量是 %g kg/kg(dry air)\n',dw);