//Cubic spline curve fitting
import java.io.*;
import javax.swing.*;

class NA49
{ 
       
  public static double [] thomas(double f[],double e[],double g[],double r[])
  {
  // 3 band matrix system of equation solving algorithm
  int n=f.length;
  double x[]=new double[n];
  for(int k=1;k<n;k++)
     {e[k]=e[k]/f[k-1];
      f[k]=f[k]-e[k]*g[k-1];
     }
  for(int k=1;k<n;k++)
     {r[k]=r[k]-e[k]*r[k-1];
     }
  x[n-1]=r[n-1]/f[n-1];   
  for(int k=(n-2);k>=0;k--)
     {x[k]=(r[k]-g[k]*x[k+1])/f[k];}
  return x;      
  }
  
  
  public static double [][] cubic_spline(double xi[],double yi[],double c0,double cn)
  {
  int n=xi.length;
  double h[]=new double[n];
  double w[]=new double[n];
  double f[]=new double[n];
  double e[]=new double[n];
  double g[]=new double[n];
  double d[]=new double[n];
  double x[]=new double[n];
  double S[][]=new double[4][n];
  int k;	  
  for(k=0;k<(n-1);k++)
     {h[k]=xi[k+1]-xi[k];
      w[k]=(yi[k+1]-yi[k])/h[k];      
	 }
  d[0]=c0;
  d[n-1]=cn;
  for(k=1;k<(n-1);k++)
     {d[k]=6.0*(w[k]-w[k-1]);}
  f[0]=1.0;
  f[n-1]=1.0;
  g[0]=0.0;
  g[n-1]=0.0;
  e[0]=0.0;
  e[n-1]=0.0;
  for(k=1;k<(n-1);k++)
  {f[k]=2.0*(h[k]+h[k-1]);e[k]=h[k-1];g[k]=h[k];}
  S[2]=thomas(f,e,g,d);
  S[3]=xi;
  for(k=0;k<(n-1);k++)
    {S[0][k]=(6.*yi[k+1]-h[k]*h[k]*S[2][k+1])/(6.0*h[k]);
     S[1][k]=(6.*yi[k]-h[k]*h[k]*S[2][k])/(6.0*h[k]);
    } 
  return S;    	 	     	 
  }
  
  public static double funcSpline(double S[][],double x)
  {
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y=0;
  double hk=0;
  for(int k=0;k<(n-1);k++)
     {if(S[3][k]<=x && x<=S[3][k+1]) 
      {hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       break;
       }
     }
       if(y==0 && S[3][n-2]<=x )
       {    
	   int k=n-2;
	   hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       }
  return y;
  }     
   
 public static double[][]  funcSpline(double xi[],double yi[],int numberofmidpoints)
 {
	//numberofmidpoints :  in x--o--o--x--o--o--x chain if x's are esxperimental points
	// numberofmidpoints is 2
	int n=xi.length;
	int nn=(n-1)*(numberofmidpoints+1)+1;
	double z[][]=new double[2][nn];
    double S[][]=cubic_spline(xi,yi,0,0);	
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {  z[0][k]=xi[i];z[1][k]=funcSpline(S,z[0][k]);k++;
          for(int j=0;j<numberofmidpoints;j++) 
             {dx=(xi[i+1]-xi[i])/((double)numberofmidpoints+1.0);
          z[0][k]=z[0][k-1]+dx;z[1][k]=funcSpline(S,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=funcSpline(S,z[0][k]);
    return z;
 } 
  // Türev formülleri
// =================       
  public static double dfSpline(double xi[],double yi[],double c0,double cn,double x)
  { //kübik şerit türev formülü 
	double S[][]=cubic_spline(xi,yi,c0,cn);	
    return dfSpline(S,x);
  }
  
  public static double dfSpline(double S[][],double x)
  {
    //kübik şerit türev formülü 
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y=0;
  double hk=0;
  for(int k=0;k<(n-1);k++)
     {if(S[3][k]<=x && x<=S[3][k+1]) 
      {hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]-S[1][k]+(xx1*xx1*S[2][k+1]-xx2*xx2*S[2][k])/(2.0*hk);
       break;
       }
     }
       if(y==0 && S[3][n-2]<=x )
       {    
	   int k=n-2;
	   hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       }
  return y;
  } 
       
 public static double  dfSpline(double xi[],double yi[],double x)
 {
    //doğal kübik şerit türev formülü 
	double S[][]=cubic_spline(xi,yi,0,0);	
    return dfSpline(S,x);
 } 
  
  public static double intSpline(double S[][],double a,double b)
  {
    //kübik şerit integral formülü 
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y1,y2;
  double hk=0;
  double toplam=0;
  y1=0;y2=0;
  for(int k=0;k<(n-1);k++)
     { hk=(S[3][k+1]-S[3][k]);
       //System.out.println("x1="+S[3][k]+"x2="+S[3][k+1]+"a="+a+"b="+b);
      if(a>S[3][k+1])
      {toplam=0;//System.out.println("şart 1 toplam = "+toplam);
      }//şart1
	  else if(S[3][k]<=a && a<=S[3][k+1] && b>S[3][k+1])//şart 2 
       {
       xx1=(a-S[3][k]);
       xx2=(S[3][k+1]-a);
       y1=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       xx1=hk;
       xx2=0;
       y2=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       toplam+=(y2-y1);
       //System.out.println("şart 2toplam = "+toplam+"y1="+y1+"y2="+y2);
       }
       else if(S[3][k]<=a && a<=S[3][k+1] && S[3][k]<=b && b<=S[3][k+1])//şart 3 
       {
       xx1=(a-S[3][k]);
       xx2=(S[3][k+1]-a);
       y1=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       xx1=(b-S[3][k]);
       xx2=(S[3][k+1]-b);
       y2=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       toplam+=(y2-y1);
       //System.out.println("şart 3 toplam = "+toplam+"y1="+y1+"y2="+y2);
       }
       else if(a<S[3][k] && b>=S[3][k+1]) //şart 4
       {
       xx1=0;
       xx2=hk;
       y1=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       xx1=hk;
       xx2=0;
       y2=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
       toplam+=(y2-y1);
       //System.out.println("şart 4 toplam = "+toplam+"y1="+y1+"y2="+y2);
       }
       else if(S[3][k]<=b && b<=S[3][k+1])//şart 5
       {xx1=0;
        xx2=hk;
        y1=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
        xx1=(b-S[3][k]);
        xx2=(S[3][k+1]-b);
        y2=S[0][k]*xx1*xx1/2.0-S[1][k]*xx2*xx2/2.0+(xx1*xx1*xx1*xx1*S[2][k+1]-xx2*xx2*xx2*xx2*S[2][k])/(24.0*hk);
        toplam+=(y2-y1);
        //System.out.println("şart 5 toplam = "+toplam+"y1="+y1+"y2="+y2);
        }
        else break;
     }
  return toplam;
  }
  
  public static double s0(double S[][],double a,double b)
  {
  //kübik şerit integral formülü integral(S(x)/x)
  int n=S[0].length;
  double xx1=0;  //x
  //ak S[0][k]
  //bk S[1][k]
  //ck S[2][k]
  //xk s[3][k]
  double xx2=0;  
  double y1,y2;
  double hk=0;
  double toplam=0;
  y1=0;y2=0;
  double xk1=0,xk=0,ak=0,bk=0,ck=0,ck1=0;
  for(int k=0;k<(n-1);k++)
     { hk=(S[3][k+1]-S[3][k]);
       xk=S[3][k];
       xk1=S[3][k+1];
       ak=S[0][k];
       bk=S[1][k];
       ck=S[2][k];
       ck1=S[2][k+1];
       //System.out.println("x1="+S[3][k]+"x2="+S[3][k+1]+"a="+a+"b="+b);
      if(a>xk1)
      {toplam=0;//System.out.println("şart 1 toplam = "+toplam);
      }//şart1
	  else if(xk<=a && a<=xk1 && b>xk1)//şart 2 
       {
       xx1=(xk1-a);
       xx2=(xk1/a);
       y1=(ak-bk)*xx1+(bk*xk1-ak*xk)*Math.log(xx2);
       xx1=xk1;
       xx2=a;
       y1+=ck1*((xx1*xx1*xx1-xx2*xx2*xx2)/3.0-
       1.5*(xx1*xx1-xx2*xx2)*xk+
       3.0*xk*xk*(xx1-xx2)-
       xk*xk*xk*Math.log(xx1/xx2))/6/hk;
       y1+=ck*(xk1*xk1*xk1*Math.log(xx1/xx2)-
       3.0*xk1*xk1*(xx1-xx2)+1.5*xk1*(xx1*xx1-xx2*xx2)-
       1.0/3.0*(xx1*xx1*xx1-xx2*xx2*xx2))/6/hk;
       toplam+=y1;
       //System.out.println("şart 2toplam = "+toplam+"y1="+y1+"y2="+y2);
       }
       else if(S[3][k]<=a && a<=S[3][k+1] && S[3][k]<=b && b<=S[3][k+1])//şart 3 
       {
       xx1=(b-a);
       xx2=(b/a);
        y1=(ak-bk)*xx1+(bk*xk1-ak*xk)*Math.log(xx2);
       xx1=b;
       xx2=a;
       y1+=ck1*((xx1*xx1*xx1-xx2*xx2*xx2)/3.0-
       1.5*(xx1*xx1-xx2*xx2)*xk+
       3.0*xk*xk*(xx1-xx2)-
       xk*xk*xk*Math.log(xx1/xx2))/6/hk;
       y1+=ck*(xk1*xk1*xk1*Math.log(xx1/xx2)-
       3.0*xk1*xk1*(xx1-xx2)+1.5*xk1*(xx1*xx1-xx2*xx2)-
       1.0/3.0*(xx1*xx1*xx1-xx2*xx2*xx2))/6/hk;
       toplam+=y1;
       }
       else if(a<S[3][k] && b>=S[3][k+1]) //şart 4
       {
       xx1=(xk1-xk);
       xx2=(xk1/xk);
        y1=(ak-bk)*xx1+(bk*xk1-ak*xk)*Math.log(xx2);
       xx1=xk1;
       xx2=xk;
       y1+=ck1*((xx1*xx1*xx1-xx2*xx2*xx2)/3.0-
       1.5*(xx1*xx1-xx2*xx2)*xk+
       3.0*xk*xk*(xx1-xx2)-
       xk*xk*xk*Math.log(xx1/xx2))/6/hk;
       y1+=ck*(xk1*xk1*xk1*Math.log(xx1/xx2)-
       3.0*xk1*xk1*(xx1-xx2)+1.5*xk1*(xx1*xx1-xx2*xx2)-
       1.0/3.0*(xx1*xx1*xx1-xx2*xx2*xx2))/6/hk;
       toplam+=y1;
       //System.out.println("şart 4 toplam = "+toplam+"y1="+y1+"y2="+y2);
       }
       else if(S[3][k]<=b && b<=S[3][k+1])//şart 5
       {       xx1=(b-xk);
       xx2=(b/xk);
       y1=(ak-bk)*xx1+(bk*xk1-ak*xk)*Math.log(xx2);
       xx1=b;
       xx2=xk;
       y1+=ck1*((xx1*xx1*xx1-xx2*xx2*xx2)/3.0-
       1.5*(xx1*xx1-xx2*xx2)*xk+
       3.0*xk*xk*(xx1-xx2)-
       xk*xk*xk*Math.log(xx1/xx2))/6/hk;
       y1+=ck*(xk1*xk1*xk1*Math.log(xx1/xx2)-
       3.0*xk1*xk1*(xx1-xx2)+1.5*xk1*(xx1*xx1-xx2*xx2)-
       1.0/3.0*(xx1*xx1*xx1-xx2*xx2*xx2))/6/hk;
       toplam+=y1;
        //System.out.println("şart 5 toplam = "+toplam+"y1="+y1+"y2="+y2);
        }
        else break;
     }
  return toplam;
  }
    
  public static double intSpline(double xi[],double yi[],double c0,double cn,double a,double b)
  {double S[][]=cubic_spline(xi,yi,c0,cn);	
   return intSpline(S,a,b);
  }
  
  //Optimisation
  public static double golden(double xi[],double yi[],double c0,double cn,double a,double b)
  {double S[][]=cubic_spline(xi,yi,c0,cn);	
   return golden(S,a,b);
  }
  
public static double golden(double S[][],double a,double b)
{
//  find the minimum of the function
//  note  maximum f(x) = minimum (-f(x))
    double epsilon=1.0e-10;
    double delta=1.0e-8;
    int print=0;
    double r1 = (Math.sqrt(5.0)-1.0)/2.0; // golden ratio
    double r2 = r1*r1;
    double h = b - a;
    double ya = -funcSpline(S,a);
    double yb = -funcSpline(S,b);
    double c = a + r2*h;
    double d = a + r1*h;
    double yc = -funcSpline(S,c);
    double yd = -funcSpline(S,d);
    int k = 1;
    double dp,dy,p,yp;
    while ((Math.abs(yb-ya)>epsilon) || (h>delta))
    {
      k++;
      if (yc<yd)
        {
        b = d;
        yb = yd;
        d = c;
        yd = yc;
        h = b - a;
        c = a + r2 * h;
        yc = -funcSpline(S,c);
        }
      else
       {
         a = c;
         ya = yc;
         c = d;
         yc = yd;
         h = b - a;
         d = a + r1 * h;
         yd = -funcSpline(S,d);
       }//end of if
    }//end of while
    dp = Math.abs(b-a);
    dy = Math.abs(yb-ya);
    p = a;
    yp = ya;
    if (yb<ya)
    {
      p = b;
      yp = yb;
    }
    if(print==1)
    {System.out.println("x min = "+p+"ymin = "+yp+"error of x ="+dp+"error of y"+dy); }
    return p;
}
    
  public static void main(String args[]) throws IOException
  {
  //double x[]={0.9,1.3,1.9,2.1,2.6,3.0,3.9,4.4,4.7,5.0,6.0,7.0,8.0,9.2,10.5,11.3,11.6,12.0,12.6,13.0,13.3};
  //double y[]={1.3,1.5,1.85,2.1,2.6,2.7,2.4,2.15,2.05,2.1,2.25,2.3,2.25,1.95,1.4,0.9,0.7,0.6,0.5,0.4,0.25};  
  String s1=JOptionPane.showInputDialog("file name : ");	  
  double c[][]=Text.readDoubleT(s1);
  double z3[][]=funcSpline(c[0],c[1],20);
  System.out.println("Cubic spline \n"+Matrix.toStringT(z3));
  Text.print(Text.T(z3),"Curve fitting with cubic spline interpolation polinomials");
  Plot pp=new Plot(c[0],c[1]);
  pp.setPlabel("Curve fitting with cubic spline interpolation polinomials");
  pp.setPlotType(0,20); 
  pp.addData(z3[0],z3[1]);
  pp.plot();
  }
}
