import javax.swing.*;
public class peng_robinson
{
//T derece K
//P bar
//v m^3/kmol

double R; // J/kmol K
double M; // kg/kmol
double Tc;  //derece K
double Pc;  //bar
double omega;
double fomega;
double u,w,b; 
double Tref;
double vref;
double uref;
public peng_robinson()
{R=8314.5; // J/kmol K
 M=44.01; // kg/kmol
 Tc=304.1;  //derece K
 Pc=73.8;  //bar
 omega=0.239;
 fomega=0.37464+1.54226*omega-0.26992*omega*omega;
 u=2;
 w=-1;
 b=0.0778*R*Tc/(Pc*1.0e5); 
 Tref=273.15; //derece K
 vref=22.4140732; //m^3/kmol
 uref=0.0; //J/kmol
}

public double a(double T)
{ double xx=(1+fomega*(1-Math.sqrt(T/Tc)));
double aa=0.45724*R*R*Tc*Tc/(Pc*1e5)*xx*xx;
return aa;
}

public double dadt(double T)
{ double xx=fomega*(1-Math.sqrt(T/Tc));
double aa=0.45724*R*R*Tc*Tc/(Pc*1e5)*(1+xx)*xx;
return aa;
}

public double P(double T,double v)
{double aa=R*T/(v-b)-a(T)/(v*v+u*b*v+w*b*b);
 return aa*1e-5;
}

public double func_v(double T,double v)
{double aa=(a(T)-dadt(T)*T)/(v*v+u*b*v+w*b*b);
return aa;
}

public double Cp(double T)
{double c[]={19.8,7.344e-2,-5.602e-5,1.715e-8};
double aa= c[0]+c[1]*T+c[2]*T*T+c[3]*T*T*T;
return aa;
}
public double Cv(double T)
{
double aa=Cp(T)-R*1e-3;
return aa;
}

public double integCv(double T)
{ //KJ/kmol
double c[]={(19.8-R*1e-3),7.344e-2,-5.602e-5,1.715e-8};
double aa= c[0]*(T-Tref)+c[1]*(T*T-Tref*Tref)/2.0+c[2]*(T*T*T-Tref*Tref*Tref)/3.0+c[3]*(T*T*T*T-Tref*Tref*Tref*Tref)/4.0;
return aa;
}

public static double[][] gauss_legendre_coefficients(double x1,double x2,int n)
{
	// calculates legendre gauss-coefficients as coefficients of the integral
	// for n terms
	// x1,x2 integral limits
	// n : number of roots
	double EPS=3.0e-15;
	int m,j,i;
	double z1,z,xm,xl,pp,p3,p2,p1;
    double a[][]=new double[2][n];
	m=(n+1)/2;
	xm=0.5*(x2+x1);
	xl=0.5*(x2-x1);
	for (i=1;i<=m;i++)  {
		z=Math.cos(Math.PI*((i-0.25)/(n+0.5)));
		do {
			p1=1.0;
			p2=0.0;
			for (j=1;j<=n;j++) {
				p3=p2;
				p2=p1;
				p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			pp=n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z=z1-p1/pp;
		} while (Math.abs(z-z1) > EPS);
		a[0][i-1]=xm-xl*z;
		a[0][n-i]=xm+xl*z;
		a[1][i-1]=2.0*xl/((1.0-z*z)*pp*pp);
		a[1][n-i]=a[1][i-1];
			}
return a;
}


public double integral(double T,double v)
{
int n=20;
double a[][]=new double[2][n];
a=gauss_legendre_coefficients(vref,v,n);
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*func_v(T,a[0][i]);}
return z;
}

public double u(double T,double v)
{
	double x1=integCv(T);
	double x2=integral(T,v);
	System.out.println("x1="+x1+"x2="+x2);
return uref+x1+x2*1e-3;

}

public double h(double T,double v)
{
return u(T,v)+P(T,v)*v;
}

public double funcP(double T,double Pi,double v)
{ //Pi bar
  return P(T,v)-Pi;
}
public double dfuncP(double T,double Pi,double v)
{ // derivative formula
  double h=0.02;
  double hh=1.0/h; 
  return (-funcP(T,Pi,(v+2.0*h))+8.0*funcP(T,Pi,(v+h))-8.0*funcP(T,Pi,(v-h))+funcP(T,Pi,(v-2.0*h)))/12.0*hh;
}
// entropi denklemi

public double dpdt(double T,double v)
{
double dadtt=dadt(T);
return 	(R/(v-b))-(dadtt/(v*v+u*b*v+w*b*b));
}

public double dpdtintegral(double T , double v)
{
int n=20;
double a[][]=new double[2][n];
a=gauss_legendre_coefficients(v,vref,n);
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*dpdt(T,a[0][i]);}
return z;
}

public double CvT(double T)
{
double dd=(Cp(T)-R*1e-3)/T;
return dd;	
}

public double integralCvT (double T)
{
int n=20;
double a[][]=new double[2][n];
a=gauss_legendre_coefficients(T,Tref,n);
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*CvT(a[0][i]);}
return z;
}

public double s (double T,double v)
{
  return sref+integralCvT(T)+dpdtintegral(T,v);
}

public  double v(double T,double Pi)
{ // Newton-Raphson
double v=R*T/(Pi*1e5);
int nmax=100;
double tolerance=1.0e-10;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
  fx= funcP(T,Pi,v);
  dfx=dfuncP(T,Pi,v);
  v-=fx/dfx;
  System.out.println("i="+i+"v="+v+"fx="+fx+"dfx = "+dfx);
  if(Math.abs(fx)<tolerance) { return v;}
}
JOptionPane.showMessageDialog(null,"Warning : Maximum number of iteration is exceeded \n"+
   " çözüm geçerli olmıyabilir","MAXİMUM ITERATION NUMBER WARNING",JOptionPane.WARNING_MESSAGE);
return v;
}


}