#include
#include
#include

double?tb=111.63tc=190.58?pc=4604?zc=0.228?Mn=16.000?w=0.011;?//定義物質(zhì)物性臨界參數(shù)，注意w是偏心因子
double?R=8.3144;
double?Rg1;
Rg1=R/Mn;

//定義關(guān)于PR方程的全局變量

double?zmin（double?xdouble?ydouble?z）
{
	x=（x	x=（x	return?x;
}??//比較三個(gè)數(shù)的最小值的函數(shù)
//*******************************************************
double?zmax（double?xdouble?ydouble?z）
{
	x=（x>y）?x:y;
	x=（x>z）?x:z;
	return?x;
}???//比較三個(gè)數(shù)的最大值的函數(shù)
//*******************************************************
double?zfg（double?tdouble?pint?i）??//求解關(guān)于z的三次方程，使用PR方程，
{
	double?zmax（double?xdouble?ydouble?z）;
	double?zmin（double?xdouble?ydouble?z）;
	double?aalphamb;?????????//定義PR狀態(tài)方程中用到的常數(shù)。
	double?ABprtr;???????????????//定義A、B，在組建關(guān)于Z的三次方程的時(shí)候，其常數(shù)項(xiàng)和A、B有關(guān)。
	double?z1z2z3;????????????????//定義方程的三個(gè)解
	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;???//定義關(guān)于z的三次方程的四個(gè)常數(shù)項(xiàng)。按照原始公式字母，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;????//定義關(guān)于盛金公式的三個(gè)判別式和一個(gè)用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;
	//以下對(duì)盛金公式求解不同根進(jìn)行討論。
	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）;//注意公式別寫錯(cuò)
	if?（i）return?（zg）;?else?return?（zf）;
//	}//這里的i是起到控制返回值的作用。
}

double?df（double?tdouble?p）????????????????//求逸度差的函數(shù)
{
	double?aalphamb;?????????//定義PR狀態(tài)方程中用到的常數(shù)。
	double?ABprtr;???????????????//定義A、B，在組建關(guān)于Z的三次方程的時(shí)候，其常數(shù)項(xiàng)和A、B有關(guān)。
	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）????????????//求理想-實(shí)際焓差的函數(shù)
{
	double?aalphamb;?????????//定義PR狀態(tài)方程中用到的常數(shù)。
	double?ABprtr;???????????????//定義A、B，在組建關(guān)于Z的三次方程的時(shí)候，其常數(shù)項(xiàng)和A、B有關(guān)。
	pr=p/pc;?tr=t/tc;
	a=0.457235*（pow（（Rg*tc）2.0）/pc