import java.io.*;
import java.util.*;
import javax.swing.*;


abstract class f_xr
{
  abstract double func(double x,int equation_ref);
}

class fh extends f_xr
{
	double func(double t,int i)
	{	
      double T=t+273.15;
	  double xx=0.0;
      if(i==0) xx=T; 
      else if(i==1) xx=1.0e-3/2.0*T*T;
      else if(i==2) xx=-1.0e5/T; 
      else if(i==3) xx=1.0e-6/3.0*T*T*T; 
      return xx;
	}      
}    

class fcp extends f_xr
{
	double func(double t,int i)
	{ // Cp(T) = xa[i]+xb[i]*1e-3*T+xc[i]*1.0e5/T^2+xd[i]*1e-6*T^2	
      double T=t;
	  double xx=0.0;
      if(i==0) xx=1.0; 
      else if(i==1) xx=1.0e-3*T;
      else if(i==2) xx=1.0e5/(T*T); 
      else if(i==3) xx=1.0e-6*T*T; 
      return xx;
	}      
} 

class Gas_GEN_Data
{
//this class reads and store and proses ideal gas data 
public double M;        // kg/kmol
public double h0;       // KJ/kmol at temperature a[0][0]=Tref
public double s0;       // KJ/kmolK at temperature a[0][0]=Tref
public double hf;       // KJ/kmol at 298 K
public double R=8.3145; // KJ/kgmolK
public double a[][];
public String GasName;
public String openName;
public String ekbilgi;
int fitoption;
f_xr f1;
f_xr f2;

public double cpl[][];  // ideal gas heat capacitiy
public double visl[][]; // saturated liquid viscosity
public double kl[][]; // saturated liquid thermal conductivity

//General partial continious least square method
   public Gas_GEN_Data(int fitoptioni,String GasNamei,String openNamei,double ai[][],int ngroup,double Mi, double h0i,double s0i,double hfi)
   {
   //fitoption : 0   Cp
   //other           h	   
   fitoption=fitoptioni;
   GasName=GasNamei;
   a=ai;
   ekbilgi=" ";
   openName=openNamei;
   int n=a.length;
   int i,j;
   if(fitoption==0) {f1=new fcp(); } 
   else             {f1=new fh(); }
   f2=new fcp(); 
   cpl=EKKgeneral(f1,a[0],a[2],ngroup);
   visl=EKKgeneral(f2,a[0],a[3],ngroup);   
   kl=EKKgeneral(f2,a[0],a[4],ngroup);
   M=Mi;
   h0=h0i;
   s0=s0i;      
   hf=hfi;    
   }
   	  
  public static double[] pivotlugauss(double a[][],double b[])
  { //gauss eliminiation with partial pivoting
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gausselimination
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //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;
	}
	//solving of gauss elimination 
    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];
	}
  }
  //back 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[][] EKKgeneral(f_xr f, double xi[],double yi[],int ngroup)
{
	//ngroup : number of data in each group
	//negri  : degree of curve fitting
	//groupnumber : number of group
	int negri=3;
	int n=negri;
	double y[]=new double[xi.length];
	for(int i=0;i<xi.length;i++)
	   {y[i]=yi[i];}
if(ngroup<negri) ngroup=negri;
int Nxi=xi.length; //total number of data
int groupnumber=Nxi/ngroup;
int ngroup1=ngroup+Nxi%ngroup;
//
int i,j,k;
int np1=4;
int np5=np1+2;
double A[][];
A=new double[np1][np1];
double B[];
B=new double[np1];
double X[][];
X=new double[groupnumber][np5];
System.out.println("group number = "+ groupnumber);
double Y[]=new double[np1];
int l=0;
double max=0;
//
for(l=0;l<groupnumber-1;l++)
{
for(i=0;i<np1;i++)
  {  B[i]=0;
	 for(j=0;j<np1;j++)
     {
	 A[i][j]=0.0;    
     for(k=l*ngroup;k<=(l+1)*ngroup;k++) {A[i][j]+=f.func(xi[k],i)*f.func(xi[k],j);}
     }
     for(k=l*ngroup;k<=(l+1)*ngroup;k++) B[i]+= f.func(xi[k],i)*y[k];
  } 
Y=pivotlugauss(A,B);//pivotlugauss(A,B);

int ii;
for(ii=0;ii<Y.length;ii++) {X[l][ii]=Y[ii];}
X[l][ii]=xi[l*ngroup];
X[l][ii+1]=xi[(l+1)*ngroup];
for(i=0;i<n+1;i++)
 if(Math.abs(X[l][i]) > max) max = Math.abs(X[l][i]);
for(i=0;i<n+1;i++)
 if((Math.abs(X[l][i]/max) > 0) && (Math.abs(X[l][i]/max) < 1.0e-100))X[l][i]=0;
}
l=groupnumber-1;
int k1=l*ngroup;
int k2=(l+1)*ngroup+Nxi%ngroup;
for(i=0;i<np1;i++)
  {  B[i]=0;
	 for(j=0;j<np1;j++)
     {
	 A[i][j]=0.0;   
     for(k=k1;k<k2;k++) {A[i][j]+=f.func(xi[k],i)*f.func(xi[k],j);}
     }
     for(k=k1;k<k2;k++) {B[i]+= f.func(xi[k],i)*y[k];}
  } 
Y=pivotlugauss(A,B);
int ii;
for(ii=0;ii<Y.length;ii++) {X[l][ii]=Y[ii];}
X[l][ii]=xi[k1];
X[l][ii+1]=xi[k2-1];

return X;
}

public static double funcEKKgeneral(f_xr f,double e[][],double x)
{
// this function calculates the value of
// least square curve fitting function
int n1=e.length;
int n2=e[0].length;
int n3=n2-2;
System.out.println("n1="+n1+"n2="+n2+"n3="+n3);
double ff=0;
double xlow,xhigh;
for(int j=0;j<n1;j++)
{
xlow=e[j][n3];
xhigh=e[j][n3+1];
System.out.println("xlow="+xlow+"xhigh="+xhigh);
if((x>=xlow) && (x<=xhigh))
{	
ff=0;
for(int i=0;i<n3;i++)
  {ff+=e[j][i]*f.func(x,i);System.out.println("ff="+ff+"i="+i);}
break;
}
}   
return ff;
}

public static double[][] xyfit(f_xr f,double x[],double yi[],double e[][])
{
//calculates absolute square root error of a least square approach
double n=x.length;
 int k;
 double total=0;
 double yy=0;
 double z[][]=new double[2][x.length];
 for(k=0;k<n;k++)
 {   z[0][k]=x[k];
	 z[1][k]=funcEKKgeneral(f,e,x[k]);
 }
return z;
}


public static double[][] dxyfit(f_xr f,double x[],double yi[],double e[][])
{
//calculates absolute square root error of a least square approach
double n=x.length;
 int k;
 double total=0;
 double yy=0;
 double z[][]=new double[2][x.length];
 for(k=0;k<n;k++)
 {   z[0][k]=x[k];
	 z[1][k]=yi[k]-funcEKKgeneral(f,e,x[k]);
 }
return z;
}

public static double error(f_xr f,double x[],double yi[],double e[][])
{
//calculates absolute square root error of a least square approach
double n=x.length;
 int k;
 double total=0;
 double yy=0;
 for(k=0;k<n;k++)
 {
	 yy=yi[k]-funcEKKgeneral(f,e,x[k]);
 total+=yy*yy;
 }
total=Math.sqrt(total/(n-1));
return total;
}
 public static double[][]  errorEKKPlot(f_xr f,double xi[],double yi[],int ngroup)
 {	
	//ngroup : Number of data in each group
	int n=xi.length;
	int l;
	double z[][];
    double E[][];//=EKKgeneral(xi,yi,polinomkatnumber,ngroup);	
    double dx=0;
    int k=0;
    int i; 
	E=EKKgeneral(f,xi,yi,ngroup);      
    z=dxyfit(f,xi,yi,E);
    Plot pp=new Plot(z[0],z[1]);
    pp.setPlabel("Cp=a0+a1*1e-3*T+a2*1e5/T^2+a3*1e-5*T^2  Tl<=T<=Th");
    pp.setXlabel("T degree K");
    pp.setYlabel("Error Cp KJ/kmolK");    
    pp.plot();
    return z;
 }
 
 public static double[][]  funcEKKPlot(f_xr f,double xi[],double yi[],int ngroup)
 {	
	int n=xi.length;
	int l;
	double z[][];
    double E[][];
    double dx=0;
    int k=0;
    int i; 
	E=EKKgeneral(f,xi,yi,ngroup);      
    z=xyfit(f,xi,yi,E);
    Plot pp=new Plot(z[0],z[1]);
       pp.setPlabel("Cp=a0+a1*1e-3*T+a2*1e5/T^2+a3*1e-5*T^2  Tl<=T<=Th");
       pp.setXlabel("T temperature degree K");
       pp.setYlabel("Cp specific heat KJ/kmolK");      
     pp.addData(xi,yi,1);  
    pp.plot();
    return z;
 }
 
 public static double[][]  funcEKKgeneral(f_xr f,double xi[],double yi[],int ngroup,int subintervalnumber)
 {
	int n=xi.length;
	int nn=(n-1)*(subintervalnumber+1)+1;
	double z[][]=new double[2][nn];
    double E[][]=EKKgeneral(f,xi,yi,ngroup);	
    Text.print(E,"EfuncEKK");
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {z[0][k]=xi[i];z[1][k]=funcEKKgeneral(f,E,z[0][k]);
       for(int j=0;j<subintervalnumber;j++) 
       {dx=(xi[i+1]-xi[i])/((double)subintervalnumber+1.0);
       k++;
       z[0][k]=z[0][k-1]+dx;z[1][k]=funcEKKgeneral(f,E,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=funcEKKgeneral(f,E,z[0][k]);
    return z;
 } 

   public double Cp_l(double t)
   { return funcEKKgeneral(f2,cpl,t);}     

   public double viscosity_l(double t)
   {return funcEKKgeneral(f2,visl,t);}

   public double k_l(double t)
   {return funcEKKgeneral(f2,kl,t);}
 
   
  public static void main(String args[]) throws IOException
  {
  String s1=JOptionPane.showInputDialog("name of the input file : ");	  
  double c[][]=Text.readDoubleT(s1);
  int n=c[0].length;
  Text.print(Text.T(c));
  double d[][]=EKKgeneral(new fcp(), c[0],c[1],n);
  funcEKKPlot(new fcp(), c[0],c[1],n);
  Text.print(d);
  errorEKKPlot(new fcp(), c[0],c[1],n);
  System.out.println(Matrix.toString(d));
  }

}