PR方程计算纯物质参数

#include
#include
#include

double tb=111.63tc=190.58 pc=4604 zc=0.228 Mn=16.000 w=0.011; //定义物质物性临界参数，注意w是偏心因子
double R=8.3144;
double Rg1;
Rg1=R/Mn;

//定义关于PR方程的全局变量

double zmin（double xdouble ydouble z）
{
	x=（x	x=（x	return x;
}  //比较三个数的最小值的函数
//*******************************************************
double zmax（double xdouble ydouble z）
{
	x=（x>y）?x:y;
	x=（x>z）?x:z;
	return x;
}   //比较三个数的最大值的函数
//*******************************************************
double zfg（double tdouble pint i）  //求解关于z的三次方程，使用PR方程，
{
	double zmax（double xdouble ydouble z）;
	double zmin（double xdouble ydouble z）;
	double aalphamb;         //定义PR状态方程中用到的常数。
	double ABprtr;               //定义A、B，在组建关于Z的三次方程的时候，其常数项和A、B有关。
	double z1z2z3;                //定义方程的三个解
	pr=p/pc; tr=t/tc;

	a=0.457235*（pow（（Rg*tc）2.0）/pc）;
	m=0.37464+1.54226*w-0.26992*w*w;
	alpha=pow（（1.0+m*（1-sqrt（tr）））2.0）;
	b=0.077796*Rg*tc/pc;
	A=a*alpha*p/（Rg*Rg*t*t）;
	B=b*p/（Rg*t）;

	double k1k2k3k4;   //定义关于z的三次方程的四个常数项。按照原始公式字母，k1=a;k2=b;k3=c;k4=d
	k1=1; k2=-（1-B）; k3=A-2*B-3*B*B; k4=-（A*B-B*B-B*B*B）;

	double c1c2c3c4;    //定义关于盛金公式的三个判别式和一个用delta求解的公式（c1=A;c2=B;c3=C;c4=delta）
	c1=k2*k2-3*k1*k3; c2=k2*k3-9*k1*k4; c3=k3*k3-3*k2*k4;
	c4=c2*c2-4*c1*c3;
	//以下对盛金公式求解不同根进行讨论。
	double zgzf;
/*    if（c4>0）
	{
		double Y1Y2;
		Y1=c1*k2+3*k1*（-c2+pow（（c4）0.5））/2;
    	Y2=c1*k2+3*k1*（-c2-pow（（c4）0.5））/2;
		z1=（-k2-（pow（Y11/3）+pow（Y21/3）））/（3*k1）;
		zg=z1;zf=0.01;
	}
	if（c4=0）
	{
		double K=c2/c1;
		z1=-k2/k1+K;
		z2=z3=-K/2;
		zf=zmin（z1z2z3）;  zg=zmax（z1z2z3）;
	}
	if（c4<0）
	{*/
	double T;
	T=（2*k2*c1-3*k1*c2）/（2*pow（c11.5））;//k[0]=a;k[1]=b;k[2]=c;k[3]=d;c[0]=A;c[1]=B;c[2]=C;c[3]=delta
	double th;
	th=acos（T）;
	z3=（-k2-2*sqrt（c1）*cos（th/3））/3;
	z2=（-k2+sqrt（c1）*（cos（th/3）+sqrt（3）*sin（th/3）））/3;
	z1=（-k2+sqrt（c1）*（cos（th/3）-sqrt（3）*sin（th/3）））/3;
	zf=zmin（z1z2z3）;  zg=zmax（z1z2z3）;//注意公式别写错
	if （i）return （zg）; else return （zf）;
//	}//这里的i是起到控制返回值的作用。
}

double df（double tdouble p）                //求逸度差的函数
{
	double aalphamb;         //定义PR状态方程中用到的常数。
	double ABprtr;               //定义A、B，在组建关于Z的三次方程的时候，其常数项和A、B有关。
	pr=p/pc; tr=t/tc;
	a=0.457235*（pow（（Rg*tc）2.0）/pc）;
    m=0.37464+1.54226*w-0.26992*w*w;
	alpha=pow（（1.0+m*（1-sqrt（tr）））2.0）;
	b=0.077796*Rg*tc/pc;
	A=a*alpha*p/（Rg*Rg*t*t）;
	B=b*p/（Rg*t）;
    double zfzg;
    zf=zfg（tp0）;  zg=zfg（tp1）;
	double lnzflnzg;

	lnzf=zf-1-log（zf-B）-A/（2*B*sqrt（2））*log（（zf+（1+sqrt（2））*B）/（zf-（sqrt（2）-1）*B））;

	lnzg=zg-1-log（zg-B）-A/（2*B*sqrt（2））*log（（zg+（1+sqrt（2））*B）/（zg-（sqrt（2）-1）*B））;

    double df;
    df=lnzf-lnzg;
    return （df）;
}

double dh（double tdouble pint i）            //求理想-实际焓差的函数
{
	double aalphamb;         //定义PR状态方程中用到的常数。
	double ABprtr;               //定义A、B，在组建关于Z的三次方程的时候，其常数项和A、B有关。
	pr=p/pc; tr=t/tc;
	a=0.457235*（pow（（Rg*tc）2.0）/pc