import javax.swing.*;

public class O2
{
final static double R=8.314;// kJ/Kmol

public static double[][] leastsquarecoeffs()
{
  double [] T1={298,300,400,500,600}; 	 //K 
  double [] Cp1={29.376,29.385,30.106,31.091,32.089};// J/Kmol  
  double [] T2={600,700,800,900,1000,1100,1200}; 	  
  double [] Cp2={32.089,32.981,33.734,34.354,34.870,35.302,35.667};	  
  int  Cv1l=Cp1.length;
  int  Cv2l=Cp2.length;
  double [] Cv1 =new double[Cv1l];
  double [] Cv2 =new double[Cv2l]; 
   for(int i=0;i<Cv1l;i++)
   {Cv1[i]=Cp1[i]-R;
   }
   for(int i=0;i<Cv2l;i++)
   {Cv2[i]=Cp2[i]-R;
   }
   
   
   
double[] val1 =PolynomialLSQ(T1,Cv1,5);
double[] val2 =PolynomialLSQ(T2,Cv2,5);	
double[] val3 =PolynomialLSQ(T1,Cp1,5);	
double[] val4 =PolynomialLSQ(T2,Cp2,5);	
double[][] depo=new double[4][6]; 

   for(int i=0;i<6;i++)
   {
    depo[0][i]=val1[i];
    depo[1][i]=val2[i]; 	   
    depo[2][i]=val3[i]; 
    depo[3][i]=val4[i]; 
    }
return depo;
} 


// values for  u
public static double u_numeric_fonk1(double x)
{
	double[] v1 =new double[6];
	double a[][] =leastsquarecoeffs();
	for(int i=0;i<6;i++)
	{
	v1[i]=a[0][i];	
	}
	double f1=v1[0]+v1[1]*x+v1[2]*x*x+v1[3]*x*x*x+v1[4]*x*x*x*x+v1[5]*x*x*x*x*x;
	return f1;
}

public static double u_numeric_fonk2(double x)
{
	double[] v2 =new double[6];
	double a2[][] =leastsquarecoeffs();
	for(int i=0;i<6;i++)
	{
	v2[i]=a2[1][i];	
	}
	double f2=v2[0]+v2[1]*x+v2[2]*x*x+v2[3]*x*x*x+v2[4]*x*x*x*x+v2[5]*x*x*x*x*x;
	return f2;
}

public static double u_analytic_fonk1(double T1,double T2)
{
	double[] v1 =new double[6];
	double a1[][] =leastsquarecoeffs();
	for(int i=0;i<6;i++)
	{
	v1[i]=a1[0][i];	
	}
	double f2=(v1[0]*(T2-T1))+(v1[1]*(T2-T1)*(T2-T1)/2)+(v1[2]*(T2-T1)*(T2-T1)*(T2-T1)/3)+(v1[3]*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)/4)+(v1[4]*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)/5)+(v1[5]*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)/6);
    return f2;
}

public static double u_analytic_fonk2(double T3,double T4)
{
	double[] v2 =new double[6];
	double a1[][] =leastsquarecoeffs();
	for(int i=0;i<6;i++)
	{
	v2[i]=a1[1][i];	
	}
	System.out.println(Matrix.toString(v2));
	double f3=  (v2[0]*(T4-T3)) +  (v2[1]*(T4-T3)*(T4-T3)/2) + (v2[2]*(T4-T3)*(T4-T3)*(T4-T3)/3) +  (v2[3]*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)/4)  +  (v2[4]*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)/5)  +   (v2[5]*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)/6) ;
	return f3;
}

public static double u_newton_cotes91(double a,double b, int n)
{       //Newton-cotes integral 9 points, n=8
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
	    double h=(b-a)/n;
        double h1=h/8.0;
        double sum=0;
        int p,i;
        double x[]=new double[9];
        double f[]=new double[9];
        for(p=0;p<n;p++)
        {  {for(i=0;i<9;i++) {x[i]=a+i*h1+p*h;f[i]=u_numeric_fonk1(x[i]);}}
           sum+=h/28350.0*(989.0*f[0]+5888.0*f[1]-928.0*f[2]+10496.0*f[3]-4540.0*f[4]+
           10496.0*f[5]-928.0*f[6]+5888.0*f[7]+989.0*f[8]); 
        }
        return sum;
}


public static double u_newton_cotes92(double a,double b, int n)
{       //Newton-cotes integral 9 points, n=8
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
	    double h=(b-a)/n;
        double h1=h/8.0;
        double sum=0;
        int p,i;
        double x[]=new double[9];
        double f[]=new double[9];
        for(p=0;p<n;p++)
        {  {for(i=0;i<9;i++) {x[i]=a+i*h1+p*h;f[i]=u_numeric_fonk2(x[i]);}}
           sum+=h/28350.0*(989.0*f[0]+5888.0*f[1]-928.0*f[2]+10496.0*f[3]-4540.0*f[4]+
           10496.0*f[5]-928.0*f[6]+5888.0*f[7]+989.0*f[8]); 
        }
        return sum;
}

// values for entropy

public static double s_numeric_fonk1(double x)
{
	double[] v1 =new double[6];
	double a[][] =leastsquarecoeffs();
	for(int i=0;i<6;i++)
	{
	v1[i]=a[2][i];	
	}
	double f1=(v1[0]/x) + (v1[1]) + (v1[2]*x) + (v1[3]*x*x) + (v1[4]*x*x*x) + (v1[5]*x*x*x*x);
	return f1;
}

public static double s_numeric_fonk2(double x)
{
	double[] v2 =new double[6];
	double a2[][] =leastsquarecoeffs();
	for(int i=0;i<6;i++)
	{
	v2[i]=a2[3][i];	
	}
	double f2=(v2[0]/x)+(v2[1])+(v2[2]*x)+(v2[3]*x*x)+(v2[4]*x*x*x)+(v2[5]*x*x*x*x);
	return f2;
}

public static double s_newton_cotes91(double a,double b, int n)
{       //Newton-cotes integral 9 points, n=8
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
	    double h=(b-a)/n;
        double h1=h/8.0;
        double sum=0;
        int p,i;
        double x[]=new double[9];
        double f[]=new double[9];
        for(p=0;p<n;p++)
        {  {for(i=0;i<9;i++) {x[i]=a+i*h1+p*h;f[i]=s_numeric_fonk1(x[i]);}}
           sum+=h/28350.0*(989.0*f[0]+5888.0*f[1]-928.0*f[2]+10496.0*f[3]-4540.0*f[4]+
           10496.0*f[5]-928.0*f[6]+5888.0*f[7]+989.0*f[8]); 
        }
        return sum;
}

public static double s_newton_cotes92(double a,double b, int n)
{       //Newton-cotes integral 9 points, n=8
        // ff integrating function
        // a,b integral limits
        // n : number of sub division h=(b-a)/n;
	    double h=(b-a)/n;
        double h1=h/8.0;
        double sum=0;
        int p,i;
        double x[]=new double[9];
        double f[]=new double[9];
        for(p=0;p<n;p++)
        {  {for(i=0;i<9;i++) {x[i]=a+i*h1+p*h;f[i]=s_numeric_fonk2(x[i]);}}
           sum+=h/28350.0*(989.0*f[0]+5888.0*f[1]-928.0*f[2]+10496.0*f[3]-4540.0*f[4]+
           10496.0*f[5]-928.0*f[6]+5888.0*f[7]+989.0*f[8]); 
        }
        return sum;
}

public static double s_analytic_fonk1(double T1,double T2)
{
	double[] v1 =new double[6];
	double a1[][] =leastsquarecoeffs();
	for(int i=0;i<6;i++)
	{
	v1[i]=a1[2][i];	
	}
	double f3=(v1[0]*Math.log(T2/T1))+  (v1[1]*(T2-T1))+(v1[2]*(T2-T1)*(T2-T1)/2)  +  (v1[3]*(T2-T1)*(T2-T1)*(T2-T1)/3)  + (v1[4]*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)/4)  +(v1[5]*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)*(T2-T1)/5);
    return f3;
}

public static double s_analytic_fonk2(double T3,double T4)
{
	double[] v1 =new double[6];
	double a1[][] =leastsquarecoeffs();
	for(int i=0;i<6;i++)
	{
	v1[i]=a1[3][i];	
	}
	double f3=  (v1[0]*Math.log(T4/T3))  +  (v1[1]*(T4-T3))+(v1[2]*(T4-T3)*(T4-T3)/2)  +  (v1[3]*(T4-T3)*(T4-T3)*(T4-T3)/3)  +   (v1[4]*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)/4)  +  (v1[5]*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)*(T4-T3)/5);
    return f3;
}

public static double pressureside(double p1,double p2)
{
return R*Math.log(p2/p1);	
}

public static double[] gausswithpartialpivot(double a[][],double b[])
  { //Gauss elimination with partial pivoting
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss elimination
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //partial pivoting
	p=k;
	buyuk=Math.abs(a[k][k]);
	for(ii=k+1;ii<n;ii++)
	{  dummy=Math.abs(a[ii][k]);
	   if(dummy > buyuk) {buyuk=dummy;p=ii;}
	}
	if(p!=k)
	{ for(jj=k;jj<n;jj++)
	  { dummy=a[p][jj];
	    a[p][jj]=a[k][jj];
	    a[k][jj]=dummy;
	  }
	  dummy=b[p];
	  b[p]=b[k];
	  b[k]=dummy;
	}
	//  
    for(i=k+1;i<n;i++)
    {  carpan=a[i][k]/a[k][k];
       a[i][k]=0;
       for(j=k+1;j<n;j++)
       {   a[i][j]-=carpan*a[k][j]; }
           b[i]   =b[i]   -carpan*b[k];
	}
  }
  //backward substitution
  x[n-1]=b[n-1]/a[n-1][n-1];
  for(i=n-2;i>=0;i--)
  {
    toplam=0;
    for(j=i+1;j<n;j++)
    {  toplam+=a[i][j]*x[j];}
    x[i]=(b[i]-toplam)/a[i][i];
  }
  return x;
}	
	
public static double[] PolynomialLSQ(double xi[],double yi[],int n)
{ //Polynomial least square
int l=xi.length;
int i,j,k;
int np1=n+1;
double A[][];
A=new double[np1][np1];
double B[];
B=new double[np1];
double X[];
X=new double[np1];
for(i=0;i<n+1;i++)
  {  for(j=0;j<n+1;j++)
     {if(i==0 && j==0) A[i][j]=l;
     else  for(k=0;k<l;k++) A[i][j] += Math.pow(xi[k],(i+j));
     }
     for(k=0;k<l;k++) { if(i==0)  B[i]+= yi[k];
                       else       B[i] += Math.pow(xi[k],i)*yi[k];}
  } 
X=gausswithpartialpivot(A,B);
double max=0;
for(i=0;i<n+1;i++)
 if(Math.abs(X[i]) > max) max = Math.abs(X[i]);
for(i=0;i<n+1;i++)
 if((Math.abs(X[i]/max) > 0) && (Math.abs(X[i]/max) < 1.0e-100)) X[i]=0;
return X;
}
	
public static double funcPolynomialLSQ(double e[],double x)
{
// this function calculates the value of
// least square curve fitting function
int n=e.length;
double ff;
if(n!=0.0)
  { ff=e[n-1];
  for(int i=n-2;i>=0;i--)
   { ff=ff*x+e[i]; }
  }
else
  ff=0;
return ff;
}	
	
public static double[][]  funcPolynomialLSQ(double xi[],double yi[],int polinomkatsayisi,int aradegersayisi)
 {
	//aradegersayisi:  x--o--o--x--o--o--x zincirinde x deneysel noktalar ise
	// ara değer sayısı 2 dir
	 int n=xi.length;
	int nn=(n-1)*(aradegersayisi+1)+1;
	double z[][]=new double[2][nn];
    double E[]=PolynomialLSQ(xi,yi,polinomkatsayisi);	
    System.out.println("katsayılar :\n"+Matrix.toStringT(E));
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {z[0][k]=xi[i];z[1][k]=funcPolynomialLSQ(E,z[0][k]);k++;
       for(int j=0;j<aradegersayisi;j++) 
       {dx=(xi[i+1]-xi[i])/((double)aradegersayisi+1.0);
       z[0][k]=z[0][k-1]+dx;z[1][k]=funcPolynomialLSQ(E,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=funcPolynomialLSQ(E,z[0][k]);
    return z;
 } 
 
 
public static void main(String args[])
{
 
	double f1=u_newton_cotes91(298,600,100);
	double f2=u_newton_cotes92(600,1200,100);
	double f3=u_analytic_fonk1(298,600);
	double f4=u_analytic_fonk2(600,1200);
	double uref=0.0-R*298;     //uref=href-R*Tref!!! 
	double u_numeric=(uref+f1+f2)/1000;
	double u_analytic=(uref+f3+f4)/1000;
	/////////////////////////////////////////
	double sref=205.147; //J/Kmol referance entropy value for O2 
	double pref=1;//atm
	double p2=5;//atm
	double sint1=s_newton_cotes91(298,600,100);
	double sint2=s_newton_cotes92(600,1200,100);
	double a1=s_analytic_fonk1(298,600);
	double a2=s_analytic_fonk2(600,1200);
	double pside=pressureside(pref,p2);
	double sa=sref+a1+a2-pside;
	double sn=sref+sint1+sint2-pside;
	System.out.print(a2);
	String s1="";
	s1+=" O2 için İç Enerji değerleri u(T)=u(1200 K)  \n";
	s1+=" Analitik integral değeri  "+u_analytic+" kJ/kmol \n";
	s1+=" Sayısal integral değeri "+u_numeric+" kJ/kmol \n";
	s1+="\n";
	s1+=" O2 için Entropy değerleri s(T,p)=s(1200 K,5 atm) \n";
	s1+=" Analitik integral değeri  "+sa+" J/kmolK \n";
	s1+=" Sayısal integral değeri "+sn+" J/kmolK \n";
	JOptionPane.showMessageDialog(null,s1);
}

}