#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