//=============================================================//
//     Numerical Analysis Package in Java                      //
//     Numerical analysis calculation class Numeric            //
//                                                             //
//     Author : Dr. Turhan Coban                               //
//     N106 WSEC, University of Nebraska- Lincoln              //
//     phone : 1-402-472 3793                                  //
//     email : MCOBAN2@unl.edu                                 //
//             MCOBAN2@netscape.net                            //
//=============================================================//
import java.io.*;

abstract class f_x
{
  //single function single independent variable
  // example f=x*x
  abstract double func(double x);
}

abstract class f_xi
{
  // multifunction multi independent variable
  // a single value is returned indiced to equation_ref
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[0]-x[1]
  // func(x,1) returns the value of f[1]
  // func(x,0) returns the value of f[0]

  abstract double func(double x[],int equation_ref);
}

abstract class fi_xi
{
  // multifunction multi independent variable
  // vector of dependent variables are returned
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[0]-x[1]
  // func(x) returns the value of f[0] and f[1]
  // as a two dimensional vector

  abstract double[] func(double x[]);
}

abstract class fij_xi
{
  // multifunction multi independent variable
  // for n independent variable n*n matrix of
  // dependent variables are returned
  // example f[0][0]=x[0]+Math.sin(x[1])   f[0][1]=x[0]-x[1]
  //         f[1][0]=x[0]*x[0]-x[1]        f[1][1]=Math.exp(x[0]+x[1]*x[1]
  // func(x) returns the value of f[0][0], f[0][1]
  //                              f[1][0], f[1][1]
  // as a two dimensional matrix

  abstract double[][] func(double x[]);
}


public class Numeric
{
//This is a library of numerical methods
//
public static double derivative(f_x f_deriv,double x)
{
// This method calculates derivatives of a simple function
// accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_x
double h0=0.0256;
int i,m;
int n=7;
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
T[i][0]=( f_deriv.func(x + h[i]) - f_deriv.func( x - h[i]))/(2.0*h[i]);
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/
          (h[i]*h[i]- h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}

public static double derivative(f_xi f_deriv,int equation_ref,
double x[],int x_ref)
{
// This method calculates derivative of a function selected from a set of
// functions. Accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_xi
// df(equation_ref)/dx(x_ref)
double h0=0.0256;
int i,m;
int n=7;
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
T[i][0]=( f_deriv.func(x1,equation_ref) -
          f_deriv.func(x2,equation_ref))/(2.0*h[i]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] -
  h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i] - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}

public static double derivative(fi_xi f_deriv,
int equation_ref,double x[],int x_ref)
{
// This method calculates derivative of a function selected from a set of
// functions with respect to x[x_ref].
// Accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_xi
double h0=0.0256;
int i,m;
int n=7;
double f1[];
f1=new double[x.length];
double f2[];
f2=new double[x.length];
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
f1=f_deriv.func(x1);
f2=f_deriv.func(x2);
T[i][0]=( f1[equation_ref] - f2[equation_ref])/(2.0*h[i]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i]
  - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}

public static double integral(f_x f_xnt,double a,double b)
{
//integral of a function by using gauss-legendre quadrature
//
 double s[],w[];
 int i;
 s=new double[30];
 w=new double[30];
 s[ 0] = .15532579626752470000E-02;
 s[ 1] = .81659383601264120000E-02;
 s[ 2] = .19989067515846230000E-01;
 s[ 3] = .36899976285362850000E-01;
 s[ 4] = .58719732103973630000E-01;
 s[ 5] = .85217118808615820000E-01;
 s[ 6] = .11611128394758690000E+00;
 s[ 7] = .15107475260334210000E+00;
 s[ 8] = .18973690850537860000E+00;
 s[ 9] = .23168792592899010000E+00;
 s[10] = .27648311523095540000E+00;
 s[11] = .32364763723456090000E+00;
 s[12] = .37268153691605510000E+00;
 s[13] = .42306504319570830000E+00;
 s[14] = .47426407872234120000E+00;
 s[15] = .52573592127765890000E+00;
 s[16] = .57693495680429170000E+00;
 s[17] = .62731846308394490000E+00;
 s[18] = .67635236276543910000E+00;
 s[19] = .72351688476904450000E+00;
 s[20] = .76831207407100990000E+00;
 s[21] = .81026309149462140000E+00;
 s[22] = .84892524739665800000E+00;
 s[23] = .88388871605241310000E+00;
 s[24] = .91478288119138420000E+00;
 s[25] = .94128026789602640000E+00;
 s[26] = .96310002371463720000E+00;
 s[27] = .98001093248415370000E+00;
 s[28] = .99183406163987350000E+00;
 s[29] = .99844674203732480000E+00;
w[ 0] = .39840962480827790000E-02;
w[ 1] = .92332341555455000000E-02;
w[ 2] = .14392353941661670000E-01;
w[ 3] = .19399596284813530000E-01;
w[ 4] = .24201336415292590000E-01;
w[ 5] = .28746578108808720000E-01;
w[ 6] = .32987114941090080000E-01;
w[ 7] = .36877987368852570000E-01;
w[ 8] = .40377947614710090000E-01;
w[ 9] = .43449893600541500000E-01;
w[10] = .46061261118893050000E-01;
w[11] = .48184368587322120000E-01;
w[12] = .49796710293397640000E-01;
w[13] = .50881194874202750000E-01;
w[14] = .51426326446779420000E-01;
w[15] = .51426326446779420000E-01;
w[16] = .50881194874202750000E-01;
w[17] = .49796710293397640000E-01;
w[18] = .48184368587322120000E-01;
w[19] = .46061261118893050000E-01;
w[20] = .43449893600541500000E-01;
w[21] = .40377947614710090000E-01;
w[22] = .36877987368852570000E-01;
w[23] = .32987114941090080000E-01;
w[24] = .28746578108808720000E-01;
w[25] = .24201336415292590000E-01;
w[26] = .19399596284813530000E-01;
w[27] = .14392353941661670000E-01;
w[28] = .92332341555455000000E-02;
w[29] = .39840962480827790000E-02;
int n=30;
double z=0;
double x,y;
for(i=0;i<n;i++)
{
x=(b+a)/2.0+(b-a)/2.0*s[i];
y=f_xnt.func(x);
z+=(b-a)/2*w[i]*y;
}
for(i=0;i<n;i++)
{
x=(b+a)/2.0+(b-a)/2.0*(-s[i]);
y=f_xnt.func(x);
z+=(b-a)/2.0*w[i]*y;
}
return z;
}

public static void trap(f_x ff,double a,double h, int j,int m,double R[][])
{
    //  successive trapezoidal integration rule
    //  this program will be utilised in romberg integration
        double sum=0;
        int p;
        for(p=1;p<=m;p++)
          { sum+=ff.func(a+h*(2*p-1)); }
        R[j][0]=R[j-1][0]/2.0+h*sum;
}

public static double integral_romberg(f_x ff,double a,double b)
{
//romberg integration
int n=8;//increase n to increase accuracy
double R[][];
R=new double[n+1][n+1];
int m=1;
double h=b-a;
double close=1;
double tol=1e-40;
int j=0,k=0;
double ret=0;
R[0][0]=h/2.0*(ff.func(a)+ff.func(b));
do
{
  j++;
  h/=2.0;
  Numeric.trap(ff,a,h,j,m,R);
  m*=2;
  for(k=1;k<=j;k++)
  {
    R[j][k]=R[j][k-1]+(R[j][k-1]-R[j-1][k-1])/(Math.pow(4,k)-1.0);
    close=Math.abs(R[j-1][j-1]-R[j][j]);
    ret=R[j][k];
  }
}while(close>tol && j<n);
return ret;
}

//==================================================
//Finding Non-linear Roots of Functions

public static double newton(double x,f_x y,f_x dy)
{

// Newton - Raphson method  for single equation
// with single variable
// required function and its derivative

int nmax=500;
double tolerance=1.0e-15;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
dfx=dy.func(x);
x-=fx/dfx;
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}


public static double newton(double x,f_x y)
{
// Newton - Raphson method  for single equation
// required function only derivative is calculated
// numerically by using method derivative

int nmax=500;
double tolerance=1.0e-15;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
dfx=Numeric.derivative(y,x);
x-=fx/dfx;
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}


public static double newton_2nd_derivative(double x,f_x y,f_x dy)
{
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
//function and first order derivative is required
int nmax=500;
double dx=1e-15;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
dfx=dy.func(x);
dfxm=dy.func(x-dx);
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}

public static double newton_2nd_derivative(double x,f_x y)
{
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
// function is required, first order derivative calculated
// numerically by derivative method.
int nmax=500;
double dx=1e-3;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
dfx=Numeric.derivative(y,x);
dfxm=Numeric.derivative(y,(x-dx));
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}

public static double secant_2nd_derivative(double x,f_x y)
{
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
// function should be supplied
int nmax=500;
double dx=1.0e-3;
double dx1=2.0e-3;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
fxm1=y.func(x-dx1);
dfx=(fx-fxm)/dx;
dfxm=(fx-fxm1)/dx1;
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}


public static double bisection(double xl,double xu,f_x y)
{
//bisection method to find roots of a function y.func(x)
//function should be supplied
// defination of variables :
// xl : lower guess
// xu : upper guess
// xr : root estimate
// es : stopping criterion
// ea :approximate error
// maxit : maximum iterations
// iter  : number of iteration
double test;
double xr=0;
double es,ea;
double fxl,fxr;
int maxit=500,iter=0;
es=0.000001;
ea=1.1*es;
while((ea>es)&&(iter<maxit))
{
xr=(xl+xu)/2.0;
iter++;
if((xl+xu)!=0)
  { ea=Math.abs((xu-xl)/(xl+xu))*100;}
fxl= y.func(xl);
fxr= y.func(xr);
test= fxl*fxr;
if(test==0.0)       ea=0;
else if(test<0.0)  xu=xr;
else
  {
  xl=xr;
  }
} //end of while
return xr;
}

public static double bisectionsearch(double xl,double xu,f_x y)
{
//bisection search method to find roots of a function y.func(x)
//function should be supplied
// defination of variables :
// xl : lower guess
// xu : upper guess
// xr : root estimate
// es : stopping criterion
// ea :approximate error
// maxit : maximum iterations
// iter  : number of iteration
int maxit=100,iter=0;
double es=0.000001;
double dx=Math.abs((xu-xl)/maxit);
double dxp=dx;
double xr=xl;
double fx=y.func(xr);
double fxr=fx;
for(int i=0;i<maxit;i++)
{
xr+=dxp;
if(xr>xu) break;
fx=y.func(xr);
  if(fx*fxr<0)
  {
    if(dxp<es) return xr;
    else
    {
    xr-=dxp;
    dxp/=2.0;
    }
  }
}
return xr;
}

public static double[] newton( double x[],fi_xi y,fij_xi dy)
{
//solution of nonlinear system of equations
//by using newton raphson method.
//ti :weigth factor
//x independent value vector:
//y dependent function vector
//dy derivative of dependent function matrix
double ti=1.0;
int i;
int nmax=500;
double tolerance=1.0e-30;
int n=x.length;
double b[];
b=new double[n];
for(i=0;i<n;i++)
{
b[i]=1.0;
}
i=0;
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
    b=Matrix.multiply(Matrix.divide(y.func(x),dy.func(x)),-ti);
    x=Matrix.add(x,b);
}
if(i >= nmax) System.out.println("warning maximum number"+
              " of iterations results may not be valid");
return x;
}

public static double[] newton( double x[],fi_xi y)
{
//solution of nonlinear system of equations
//by using newton raphson method.
//this function does not require derivatives
//it is utilised method derivative to calculate derivatives
double ti=1.0;
int i,ii,jj;
int nmax=500;
double tolerance=1.0e-15;
int n=x.length;
double b[];
b=new double[n];
double dy[][];
dy=new double[n][n];
i=0;
for(i=0;i<n;i++)
{
b[i]=1.0;
}
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
  for(ii=0;ii<n;ii++)
  {
    for(jj=0;jj<n;jj++)
    {
        dy[ii][jj]=Numeric.derivative(y,ii,x,jj);
    }
  }
    b=Matrix.multiply(Matrix.divide(y.func(x),dy),-ti);
    x=Matrix.add(x,b);
}
if(i >= nmax) System.out.println("warning maximum number of iterations"+
                                 " results may not be valid"          );
return x;
}

//========= least square curve fitting methods ==========

//-------- polynomial least square curve fitting --------

public static double[] poly_least_square
(double xi[],double yi[],int n)
//polynomial least square curve fitting
//variables xi,yi vector of data points
//degree of curve fitting
{
//xi vector of independent variable
//yi vector of dependent variable
//n : degree of curve fitting
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=Matrix.divide(B,A);
//X=B/A;

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 f_least_square(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 error(double x[],double y[],double e[])
{
//calculates absolute square root error of a least square approach
double n=x.length;
 int k;
 double total=0;
 for(k=0;k<n;k++)
 {
 total+=(y[k]-f_least_square(e,x[k]))*(y[k]-f_least_square(e,x[k]));
 }
total=Math.sqrt(total);
return total;
}

//-------End least square methods --------------------
//-------Orthogonal polynomial least square curve fitting ---

public static double[][] orthogonal_polynomial_least_square
(double xi[],double fi[],int m)
//orthogonal polynomial least square curve fitting
//this method do not require any matrix solution
//variables xi,fi vector of data points wi weight functions
//m degree of curve fitting
{
int i,j,k;
int n=xi.length;
int mp2=n+2;
int mp1=n+1;
//matrix<double> p(m+2,n,0);
double p[][];
p=new double[mp2][n];
//vector<double> gamma(m+1,0);
double gamma[];
gamma=new double[mp1];
//vector<double> beta(m+1,0);
double beta[];
beta=new double[mp1];
//vector<double> omega(m+1,0);
double omega[];
omega=new double[mp1];
//vector<double> alpha(m+1,0);
double alpha[];
alpha=new double[mp1];
//vector<double> b(m+1,0);
double b[];
b=new double[mp1];
//vector<double> wi(n,1.);
double wi[];
wi=new double[n];
//matrix<double> a(3,m+1);
double a[][];
a=new double[3][mp1];
for(i=0;i<n;i++)
  {
   p[1][i]=1.0;
   p[0][i]=0.0;
  }
gamma[0]=0;
for(i=0;i<n;i++)
  {
  gamma[0]+=wi[i];
  }
beta[0]=0.0;
for(j=0;j<m+1;j++)
  {
  omega[j]=0;
  for(i=0;i<n;i++)
    {
    omega[j]+=wi[i]*fi[i]*p[j+1][i];
    }
  b[j]=omega[j]/gamma[j];
if( j != m)
    {
    alpha[j+1]=0;
    for(i=0;i<n;i++)
      {
      alpha[j+1]+=wi[i]*xi[i]*p[j+1][i]*p[j+1][i]/gamma[j];
      }
    for(i=0;i<n;i++)
      {
      p[j+2][i]=(xi[i]-alpha[j+1])*p[j+1][i]-beta[j]*p[j][i];
      }
      gamma[j+1]=0;
    for(i=0;i<n;i++)
      {
      gamma[j+1]+=wi[i]*p[j+2][i]*p[j+2][i];
      }
    beta[j+1]=gamma[j+1]/gamma[j];
    }
}//end of j
for(j=0;j<m+1;j++)
  {
  a[0][j]=b[j];
  a[1][j]=alpha[j];
  a[2][j]=beta[j];
  }
return a;
}

public static double f_orhogonal_polynomial(double a[][],double x)
{

double yy=0;
int k;
int m=a[0].length-1;
int mp2=m+2;
double q[];
q=new double[mp2];
//vector<double> q(m+2,0.0);
for(k=m-1;k>=0;k--)
  {
  q[k]=a[0][k]+(x-a[1][k+1])*q[k+1]-a[2][k+1]*q[k+2];
  yy=q[k];
  }
return yy;
}

//minimization of a function
public static double amotry(double p[][], double y[], double psum[], int ndim,
	f_xi ff, int ihi, double fac)
{
	int j;
	double fac1,fac2,ytry;
	double ptry[]=new double[ndim];
	fac1=(1.0-fac)/ndim;
	fac2=fac1-fac;
	for (j=1;j<=ndim;j++) ptry[j-1]=psum[j-1]*fac1-p[ihi-1][j-1]*fac2;
	ytry=ff.func(ptry,0);
	if (ytry < y[ihi-1]) {
		y[ihi-1]=ytry;
		for (j=1;j<=ndim;j++) {
			psum[j-1] += ptry[j-1]-p[ihi-1][j-1];
			p[ihi-1][j-1]=ptry[j-1];
		}
	}
	return ytry;
}


public static void amoeba(double p[][], double y[], int ndim, double ftol,
	f_xi ff, int nfunk)
{
	int i,ihi,ilo,inhi,j,mpts=ndim+1;
	int NMAX=5000;
	double rtol,sum,swap,ysave,ytry;
	double psum[]=new double[ndim];
	nfunk=0;
	for (j=1;j<=ndim;j++)
	{
	  for (sum=0.0,i=1;i<=mpts;i++) sum += p[i-1][j-1];
	      psum[j-1]=sum;
	}
	for (;;) {
		ilo=1;
		//ihi = y[1]>y[2] ? (inhi=2,1) : (inhi=1,2);
		if(y[0]>y[1])
		{
		  inhi=2;
                  ihi=1;
		}
		else
		{
		  inhi=1;
                  ihi=2;
		}
		for (i=1;i<=mpts;i++) {
			if (y[i-1] <= y[ilo-1]) ilo=i;
			if (y[i-1] > y[ihi-1]) {
				inhi=ihi;
				ihi=i;
			} else if (y[i-1] > y[inhi-1] && i != ihi) inhi=i;
		}
		rtol=2.0*Math.abs(y[ihi-1]-y[ilo-1])
             /(Math.abs(y[ihi-1])+Math.abs(y[ilo-1]));
		if (rtol < ftol) {
			swap=y[0];
			y[0]=y[ilo-1];
                        y[ilo-1]=swap;
			for (i=1;i<=ndim;i++)
			{
			swap=p[0][i-1];
			p[0][i-1]=p[ilo-1][i-1];
                        p[ilo-1][i-1]=swap;
                        }
			break;
		}
		if (nfunk >= NMAX) System.out.println("NMAX exceeded");
		nfunk += 2;
		ytry=amotry(p,y,psum,ndim,ff,ihi,-1.0);
		if (ytry <= y[ilo-1])
			ytry=amotry(p,y,psum,ndim,ff,ihi,2.0);
		else if (ytry >= y[inhi-1]) {
			ysave=y[ihi-1];
			ytry=amotry(p,y,psum,ndim,ff,ihi,0.5);
			if (ytry >= ysave) {
				for (i=1;i<=mpts;i++) {
					if (i != ilo) {
						for (j=1;j<=ndim;j++)
							psum[j-1]=0.5*(p[i-1][j-1]+p[ilo-1][j-1]);
                                                        p[i-1][j-1]=psum[j-1];
						y[i]=ff.func(psum,0);
					}
				}
				nfunk += ndim;
				 for (j=1;j<=ndim;j++)
	                         {
	                         for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];
	                              psum[j]=sum;
                           	 }

			}
		} else --nfunk;
	}
}

// differential equation solutions
public static double[] RK4(f_xi fp,double x0,double xn,double f0,int N)
{
//4th order Runge Kutta Method
//fp :given derivative function df/dx(f,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
double fi[]=new double[N+1];
double xi[]=new double[2];
double k[];
k=new double[4];
int i;
double x;
fi[0]=f0;
for(x=x0,i=0;x<xn;x+=h,i++)
  {
  xi[0]=x;
  xi[1]=fi[i];
   k[0]=h*fp.func(xi,0);
  xi[0]=x+h/2.0;
  xi[1]=fi[i]+k[0]/2;
   k[1]=h*fp.func(xi,0);
  xi[1]=fi[i]+k[1]/2.0;
   k[2]=h*fp.func(xi,0);
  xi[0]=x+h;
  xi[1]=fi[i]+k[2];
   k[3]=h*fp.func(xi,0);
  fi[i+1]=fi[i]+k[0]/6.0+k[1]/3.0+k[2]/3.0+k[3]/6.0;
  }
return fi;
}

public static double[][] RK4(f_xi fp,double x0,double xn,double f0[],int N)
{
//4th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[4];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[j]=fi[j-1][i]+k[1]/2.0;
     k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]+k[2];
     k[3]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/6.0+k[1]/3.0+k[2]/3.0+k[3]/6.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {
  a[j][i]=fi[j-1][i];
  }
}

return a;
}


public static double[][] RK6(f_xi fp,double x0,double xn,double f0[],int N)
{
//6th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[6];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/4.0+k[1]/4.0;
    k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]-k[1]+2.0*k[2];
    k[3]=h*fp.func(xi,j-1);
    xi[0]=x+2.0/3.0*h;
    xi[j]=fi[j-1][i]+7.0/27.0*k[0]+10.0/27.0*k[1]+1.0/27.0*k[3];
    k[4]=h*fp.func(xi,j-1);
    xi[0]=x+1.0/5.0*h;
    xi[j]=fi[j-1][i]+28.0/625.0*k[0]-1.0/5.0*k[1]+546.0/625.0*k[2]+54.0/625.0*k[3]-378/625.0*k[4];
    k[5]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/24.0+5.0*k[3]/48.0+27.0*k[4]/56.0+125.0*k[5]/336.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {
  a[j][i]=fi[j-1][i];
  }
}

return a;
}


public static double[][]
RKF45(f_xi fp,double x0,double xn,double f0,int N) throws IOException
{
// Runge Kutta - Fehlberg Method
// fp  : derivative function df/dx(f,x)  (defined as a class)
// Tol :Tolerance
//Hmax : maximum possible step size
//Hmin : minimum possible step size
//k    : Runge kutta coefficients
//Err  : error
//x[i] : input variable to fp = df/dx(f,x) = df/dx(x[i])
//j    : actual step size
//fi[][]:solution matrix
//fj[][]:solution matrix in exact size(j) (return function)
double h=(xn-x0)/N;
double fi[][]=new double[2][1000];
double Tol=2.0e-3;
double Hmin=h/64.0;
double Hmax=h*64.0;
double k[]=new double[6];
double Err;
double s;
double xi[]=new double[2];
int j=0;
fi[0][j]=x0;
fi[1][j]=f0;
while(fi[0][j] < xn )
  {
  if(fi[0][j]+h > xn) h=xn-fi[0][j];
  xi[0]=fi[0][j];
  xi[1]=fi[1][j];
   k[0]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+h/4.0;
  xi[1]=fi[1][j]+k[0]/4.0;
   k[1]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+3.0/8.0*h;
  xi[1]=fi[1][j]+3.0/32.0*k[0]+9.0/32.0*k[1];
   k[2]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+12.0/13.0*h;
  xi[1]=fi[1][j]+1932.0/2197.0*k[0]-7200.0/2197.0*k[1]+7296.0/2197.0*k[2];
   k[3]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+h;
  xi[1]=fi[1][j]+439.0/216.0*k[0]-8.0*k[1]+3680.0/513.0*k[2]-845.0/4104.0*k[3];
   k[4]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+0.5*h;
  xi[1]=fi[1][j]-8.0/27.0*k[0]+2.0*k[1]
  -3544/2565*k[2]+1859.0/4104.0*k[3]-11.0/40.0*k[4];
   k[5]=h*fp.func(xi,0);
   Err=Math.abs(1.0/360.0*k[0]-128/4275*k[2]
   -2197.0/75240*k[3]+1.0/5.0*k[4]+2.0/55.0*k[5]);
  if(Err<Tol || h<2*Hmin)
     {
     fi[0][j+1]=fi[0][j]+h;
     fi[1][j+1]=fi[1][j]+16.0/135.0*k[0]+6656.0/12825.0*k[2]+
     28561.0/56430.0*k[3]-9.0/50.0*k[4]+2.0/55.0*k[5];
     j++;
     }
  else
     {
     if(Err==0) s=0;
     else s=0.84*(Math.pow(Tol*h/Err,0.25));
     if(s<0.75 && h>(2*Hmin)  ) h/=2.0;
     if(s>1.5  && Hmax >(2*h) ) h*=2.0;
     }
  }
double fj[][]=new double[2][j+1];
for(int i=0;i<=j;i++)
  {
  fj[0][i]=fi[0][i];
  fj[1][i]=fi[1][i];
  }
return fj;
}

//End of Class Numeric
}


