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

//prototip fonksiyon sınıfı
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)
	{	
      double T=t+273.15;
	  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
{
//Genel en küçük kareler metodu en kucuk kareler metodu	
	  
  public static double[] pivotlugauss(double a[][],double b[])
  { //kısmi pivotlu gauss eleme yöntemi
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss eleme
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //pivotlama
	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;
	}
	//gauss elemeyi çözme  
    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];
	}
  }
  //geriye doğru yerine koyma
  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[][] EKKgenel(f_xr f, double xi[],double yi[],int ngurup)
{
	//ngurup : bir guruptaki veri sayısı
	//negri  : eğrinin derecesi
	//gurupsayisi : gurupların sayısı
	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(ngurup<negri) ngurup=negri;
int Nxi=xi.length; //toplam veri sayısı
int gurupsayisi=Nxi/ngurup;
int ngurup1=ngurup+Nxi%ngurup;
//
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[gurupsayisi][np5];
System.out.println("gurup sayisi = "+ gurupsayisi);
double Y[]=new double[np1];
int l=0;
double max=0;
//
for(l=0;l<gurupsayisi-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*ngurup;k<=(l+1)*ngurup;k++) {A[i][j]+=f.func(xi[k],i)*f.func(xi[k],j);}
     }
     for(k=l*ngurup;k<=(l+1)*ngurup;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*ngurup];
X[l][ii+1]=xi[(l+1)*ngurup];
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=gurupsayisi-1;
int k1=l*ngurup;
int k2=(l+1)*ngurup+Nxi%ngurup;
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 funcEKKgenel(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-4;
//System.out.println("n1="+n1+"n2="+n2+"n3="+n3);
double ff=0;
double xlow,xhigh;
//for(int j=0;j<n1;j++)
//{for(int i=0;i<n3;i++){System.out.print(e[j][i]+" ");}
//System.out.println();
//}

for(int j=0;j<n1;j++)
{
xlow=e[j][n3];
xhigh=e[j][n3+1];
if((x>=xlow) && (x<=xhigh))
{	
ff=0;
for(int i=0;i<n3;i++)
  {ff+=e[j][i]*f.func(x,i);}
break;
}
}   
return ff;
}


public static double hata(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]-funcEKKgenel(f,e,x[k]);
 total+=yy*yy;
 }
total=Math.sqrt(total/(n-1));
return total;
}

 public static void  funcEKKPlot(f_xr f,double xi[],double yi[],int ngurup)
 {	//aradegersayisi:  x--o--o--x--o--o--x zincirinde x deneysel noktalar ise
	// ara değer sayısı 2 dir
	//ngurup bir guruptaki veri sayısı
	//polinom katsayısı : polinomun derecesi
	int n=xi.length;
	int l;
	double z[][]=new double[10][n];
    double E[][];//=EKKgenel(xi,yi,polinomkatsayisi,ngurup);	
    double dx=0;
    int k=0;
    int i; 
	E=EKKgenel(f,xi,yi,ngurup); 
	double hata=hata(f,xi,yi,E);
	System.out.println("hata="+hata);
	for(int ii=0;ii<E.length;ii++)
	{for(int jj=0;jj<E[0].length;jj++){System.out.print(E[ii][jj]+ " ");}
	System.out.println(" ");
	} 
	 
    for(i=0;i<n;i++)
       {z[0][i]=xi[i];          
        z[1][i]=(yi[i]-funcEKKgenel(f,E,z[0][i]))/yi[i]*100;
        //z[2][i]=(yi[i]-makPS1.Ps1(z[0][i]))/yi[i]*100;
        //z[3][i]=(yi[i]-makPS1.Ps2(z[0][i]))/yi[i]*100;
       }
     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 derece K");
       pp.setYlabel("hata %");   
    for(int j=2;j<7;j++)
    {
	ngurup/=2; 
	E=EKKgenel(f,xi,yi,ngurup);
    for(int ii=0;ii<E.length;ii++)
	{for(int jj=0;jj<E[0].length;jj++){System.out.print(E[ii][jj]+ " ");}
	System.out.println(" ");
	} 
	hata=hata(f,xi,yi,E);
	System.out.println("hata="+hata+"ngurup = "+ngurup);     
    for(i=0;i<n;i++)
       {z[0][i]=xi[i];          
       z[j][i]=(yi[i]-funcEKKgenel(f,E,z[0][i]))/yi[i]*100;
        //z[2][i]=(yi[i]-makPS.Ps1(z[0][i]))/yi[i]*100;
        //z[3][i]=(yi[i]-makPS.Ps2(z[0][i]))/yi[i]*100;
       }
     pp.addData(z[0],z[j],j);
       }
    //pp.addData(z[0],z[2]);
    //pp.addData(z[0],z[3]);  
    pp.plot();
 }
 
 public static double[][]  funcEKKgenel(f_xr f,double xi[],double yi[],int ngurup,int aradegersayisi)
 {
	// aradegersayisi:  x--o--o--x--o--o--x zincirinde x deneysel noktalar ise
	// ara değer sayısı 2 dir
	// ngurup bir guruptaki veri sayısı
	// polinom katsayısı : polinomun derecesi
	int n=xi.length;
	int nn=(n-1)*(aradegersayisi+1)+1;
	double z[][]=new double[2][nn];
    double E[][]=EKKgenel(f,xi,yi,ngurup);	
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {z[0][k]=xi[i];z[1][k]=funcEKKgenel(f,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]=funcEKKgenel(f,E,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=funcEKKgenel(f,E,z[0][k]);
    return z;
 } 

  public static void main(String args[]) throws IOException
  {
  //çok büyük veri setlerinde veriyi dosyadan okumamız gerekir
  String s1=JOptionPane.showInputDialog("dosya adı : ");	  
  double c[][]=Text.readDoubleT(s1);
  //R134a
  //double Tc=374.23;
  //double Pc=4059.3;
  //R290 Propane
  double Tc=369.83;
  double Pc=4247.1;
 fcp f=new fcp();
  funcEKKPlot(f,c[0],c[1],1);
  //double xi[],double yi[],double Tc,double Pc,int ngurup
  //double b[][]=EKKgenel(Tc,Pc,1408);  
  // System.out.println(b.length+"  "+b[0].length);
  //System.out.println("Genel EKK Denklem katsayıları : \n"+Matrix.toString(b));
  //System.out.println("y="+funcEKKgenel(b,50.0));
  //System.out.println("hata = "+hata(c[0],c[1],b));
  /*
  //funcEKKPlot(c[0],c[1],6,59);
  //System.out.println(Matrix.toStringT(funcEKKgenel(c[0],c[1],7,20,1)));
  double z[][]=funcEKKgenel(c[0],c[1],7,2);
  System.out.println("Genel E.K.K. çözüm seti \n"+Matrix.toStringT(z)); 
  Plot pp=new Plot(c[0],c[1]);
  pp.setPlabel("Genel En küçük kareler yöntemi ile R134a sıcaklık - doyma sıvı yoğunluğu eğrisi uydurma");
  pp.setXlabel("T derece C");
  pp.setYlabel("h entalpi"); 
  pp.addData(z[0],z[1]);
  pp.plot();
  */
  }

}