写个UDF，调试通过，计算不行，大家给看看

jshln

很简单一个小程序，就是通过状态方程计算气体密度，编译时候可以通过的
但是就是一计算就出错误了，看了半天也找补出来问题，各位给看看，谢谢！
这个算法在VC6.0里面是正确的，计算结果和试验结果偏差不大的，但是在fluent环境下就是不好用，提前谢谢了。

#include "udf.h"
#define Rg 8.314
DEFINE_PROPERTY(cell_density, cell, thread)
{
double P;
double  T;
double Pc,Tc;
double Mw;
double omiga;
double A,B;
double density;
double m,at,ac;
double a,b;  
double h,z,z1;
int i;
P=C_P(cell,thread);
T=C_T(cell, thread);
Pc=72.8*101325;Tc=31+273.15;Mw=44;omiga=0.225;
m=0.37464+1.54226*omiga-0.26992*pow(omiga,2);
at=pow(1+m*(1-pow(T/Tc,0.5)),2);
ac=0.45724*Rg*Rg*Tc*Tc/Pc;
b=0.0778*Rg*Tc/Pc;
a=ac*at;
A=a*P/pow(Rg*T,2);
B=b*P/(Rg*T);
z=1,z1=0.5;
for(i=1;fabs(z-z1)>0.001;i++)
{
z1=z;
h=B/z;
z=1/(1-h)-A/B*(h/(1+2*h-h*h));
}
z=P/(z*Rg*T);
density=1000*Mw*z*1e-6;
return density;
}