//=============================================================//
//     Numerical Analysis Package in Java                      //
//     Numerical analysis calculation class Numeric            //
//                                                             //
//     Author : Dr. Turhan Coban                               //
//     TUBİTAK Ulusal Metroloji Enstitüsü                      //
//     phone : 90-262-679 5013                                 //
//     email : turhan.coban@ume.tubitak.gov.tr                 //
//             turhan_coban@yahoo.com                          //
//=============================================================//
//
import java.io.*;
import java.text.*;
import java.util.Locale;

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

abstract class f_xj
{
  // single function  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[]);
}

abstract class f_xa
{
  // single function  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 t,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
// Minimisation (Optimisation) of a function
public static double[] grads(f_xj Fn,fi_xi Gn,double P0[])
{
   return grads(Fn,Gn,P0,100,1.0e-10,1.0e-10,0);
}

public static double[] grads(f_xj Fn,fi_xi Gn,double P0[],double max1,double delta,double epsilon,int print)
{
    //GRADS   Gradient search for a minimum of non-linear multivariable function
    // Inputs
    //   Fn        name of the vector function
    //   Gn        gradient for Fn
    //   P0        starting point
    //   max1      maximum number of iterations
    //   delta     convergence tolerance for the independent variables
    //   epsilon   convergence tolerance for the dependent variable
    //   show      if show==1 the iterations are displayed
    //
    //   P0        point V0 for the minimum
    //   y0        function value  Fn(V0)
    //   h         minimum step size
    //   err       error  estimate
    //---------------------------------------------------------------------------

    int n = P0.length;
    int maxj = 20;
    double big = 1e8;
    double h = 1.0;
    double len = Matrix.norm(P0);
    double y0 = Fn.func(P0);
    if (len > 1e4) {
      h = len / 1e4;
    }
    double err = 1;
    double cnt = 0;
    double cond = 0;
    double P[] = new double[n];
    double P1[] = new double[n];
    double P2[] = new double[n];
    double Pmin[] = new double[n];
    double Y, y1, y2, ymin, hmin;
    double S[]=new double[n];
    double h0,h1,h2;
    double e0,e1,e2,d;
    int i;
    for(i=0;i<n;i++)
    { P[i]=P0[i]; }
    Y= y0;
    while (cnt<max1 && cond!=5 && (h>delta || err>epsilon))
    {
      double z[] = Gn.func(P0);
      for (i = 0; i < n; i++) {
            S[i] = z[i];
          }
      P1 = Matrix.add(P0, Matrix.multiply(S, h));
      P2 = Matrix.add(P0, Matrix.multiply(S, (2.0 * h)));
      y1 = Fn.func(P1);
      y2 = Fn.func(P2);
      cond = 0;
      int j = 0;
      while (j < maxj && cond == 0) {
        len = Matrix.norm(P0);
        if (y0 < y1) {
          for (i = 0; i < n; i++) {
            P2[i] = P1[i];
          }
          y2 = y1;
          h = h / 2.0;
          P1 = Matrix.add(P0, Matrix.multiply(S, h));
          y1 = Fn.func(P1);
        }
        else {
          if (y2 < y1) {
            for (i = 0; i < n; i++) {
              P1[i] = P2[i];
            }
            y1 = y2;
            h = 2 * h;
            P2 = Matrix.add(P0, Matrix.multiply(S, (2.0 * h)));
            y2 = Fn.func(P2);
          }
          else cond = -1;
        }
        j++;
        if (h < delta) {
          cond = 1;
        }
        if (Math.abs(h) > big || len > big) {cond = 5;}
      }
      if (cond == 5) {
        for (i = 0; i < n; i++) {
          Pmin[i] = P1[i];
        }
        ymin = y1;
      }
      else {
        d = 4 * y1 - 2 * y0 - 2 * y2; // Start of a long block:
        if (d < 0) {hmin = h * (4 * y1 - 3 * y0 - y2) / d;}
        else {
          cond = 4;
          hmin = h / 3.0;
        }
        Pmin = Matrix.add(P0, Matrix.multiply(S, hmin));
        ymin = Fn.func(Pmin);
        h0 = Math.abs(hmin);
        h1 = Math.abs(hmin - h);
        h2 = Math.abs(hmin - 2.0 * h);
        if (h0 < h) {
          h = h0;
        }
        if (h1 < h) {
          h = h1;
        }
        if (h2 < h) {
          h = h2;
        }
        if (h == 0) {
          h = hmin;
        }
        if (h < delta) {
          cond = 1;
        }
        e0 = Math.abs(y0 - ymin);
        e1 = Math.abs(y1 - ymin);
        e2 = Math.abs(y2 - ymin);
        if (e0 != 0 && e0 < err) {
          err = e0;
        }
        if (e1 != 0 && e1 < err) {
          err = e1;
        }
        if (e2 != 0 && e2 < err) {
          err = e2;
        }
        if (e0 == 0 && e1 == 0 && e2 == 0) {
          err = 0;
        }
        if (err < epsilon) {
          cond = 2;
        }
        if (cond == 2 && h < delta) {
          cond = 3;
        }
      } //end of if(cond==5)
      cnt++;
      for (i = 0; i < n; i++) {
        P0[i] = Pmin[i];
      }
      y0 = ymin;
      if(print==1)
       {System.out.println("x min = "+Matrix.toString(P0)+"/nymin = "+y0+"/nerror estimate ="+err+"/nminimum step size "+h); }
    } //end of while
    return P0;
}

  public static double quadmin(f_x f,double a,double b,double epsilon)
  {return quadmin(f,a,b,epsilon,1.0e-10,0);}

  public static double quadmin(f_x f,double a,double b)
  {return quadmin(f,a,b,1.0e-10,1.0e-10,0);}

  public static double quadmin(f_x f,double a,double b,double epsilon,double delta,int print)
  {
      // Minimum of a single variable function by using polynomial quadratures
      double p0 = a;
      int maxj = 20;
      int maxk = 30;
      double big = 1e6;
      double err = 1;
      int k = 1;
      int cond = 0;
      double  h = 1;
      double f1,p1,p2,pmin,y0,y1,y2,ymin,d,hmin,dh,h0,h1,h2;
      double e0,e1,e2;
      double p,dp,dy,yp;
      if (Math.abs(p0)>1e4)
        { h = Math.abs(p0)/1e4;}
      while (k<maxk && err>epsilon && cond !=5 )
      {
        f1 = (f.func(p0 + 0.00001) - f.func(p0 - 0.00001)) / 0.00002;
        if (f1 > 0) {h = -Math.abs(h);}
        p1 = p0 + h;
        p2 = p0 + 2 * h;
        pmin = p0;
        y0 = f.func(p0);
        y1 = f.func(p1);
        y2 = f.func(p2);
        ymin = y0;
        cond = 0;
        int j = 0;
        while (j < maxj && Math.abs(h) > delta && cond == 0)
        {
          if (y0 <= y1)
          {
            p2 = p1;
            y2 = y1;
            h = h / 2;
            p1 = p0 + h;
            y1 = f.func(p1);
          }
          else
          {
            if (y2 < y1)
            {
              p1 = p2;
              y1 = y2;
              h = 2 * h;
              p2 = p0 + 2 * h;
              y2 = f.func(p2);
            }
            else
            {  cond = -1;}
          } //end of if(y0 <= y1)
          j = j + 1;
        if (Math.abs(h) > big || Math.abs(p0) > big) {cond = 5;}
        }//end of while(j<maxj
        if (cond == 5)
        {
              pmin = p1;
              ymin = f.func(p1);
        }
        else
        {
          d = 4 * y1 - 2 * y0 - 2 * y2;
          //Start of a long block:
          if (d < 0) {
            hmin = h * (4 * y1 - 3 * y0 - y2) / d;
          }
          else
          {
            hmin = h / 3;
            cond = 4;
          }
          pmin = p0 + hmin;
          ymin = f.func(pmin);
          dh = Math.abs(h);
          h0 = Math.abs(hmin);
          h1 = Math.abs(hmin - h);
          h2 = Math.abs(hmin - 2 * h);
          if (h0 < dh) {h = h0;}
          if (h1 < dh) {h = h1;}
          if (h2 < dh) {h = h2;}
          if (h == 0) {h = Math.abs(hmin);}
          if (h < delta) {cond = 1;}
          if (Math.abs(h) > big || Math.abs(pmin) > big) {cond = 5;}
          e0 = Math.abs(y0 - ymin);
          e1 = Math.abs(y1 - ymin);
          e2 = Math.abs(y2 - ymin);
          if (e0 != 0 && e0 < err) {err = e0;}
          if (e1 != 0 && e1 < err) {err = e1;}
          if (e2 != 0 && e2 < err)  {err = e2;}
          if (e0 != 0 && e1 == 0 && e2 == 0) {err = 0;}
          if (err < epsilon) {cond = 2;}
          p0 = pmin;
          k = k + 1;
      } // End of the long else block.
      if (cond == 2 && h < delta) {cond = 3;}
  }//end of the while
  p = p0;
  dp = h;
  yp = f.func(p);
  dy = err;
  if(print==1)
  {System.out.println("x min = "+p+"ymin = "+yp+"error of x ="+dp+"error of y"+dy); }
  return p;
  }


public static double golden(f_x f,double a,double b,double epsilon)
{return golden(f,a,b,epsilon,1.0e-10,0);}
public static double golden(f_x f,double a,double b)
{return golden(f,a,b,1.0e-10,1.0e-10,0);}
public static double golden(f_x f,double a,double b,double epsilon,double delta,int print)
{
// f abstract function of class f_x (f(x) type of single function single variable)
// fibonacchi golden ratio search
    double r1 = (Math.sqrt(5.0)-1.0)/2.0; // golden ratio
    double r2 = r1*r1;
    double h = b - a;
    double ya = f.func(a);
    double yb = f.func(b);
    double c = a + r2*h;
    double d = a + r1*h;
    double yc = f.func(c);
    double yd = f.func(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 = f.func(c);
       }
      else
       {
         a = c;
         ya = yc;
         c = d;
         yc = yd;
         h = b - a;
         d = a + r1 * h;
         yd = f.func(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 double[] nelder(f_xj fnelder,double x[][],int maxiteration,double tolerance,int printlist)
{
   // Nelder mead multidimensional simplex minimisation method
   //  Nelder & Mead 1965 Computer J, v.7, 308-313.

   // Parameter identification
   // fnelder : abstract multidimensional function f(x)
   // x : for n dimension n+1 simlex points
   // maxiteration : maximum number of iteration points
   // tolerance : convergence tolerance
         int NDIMS = x.length-1;
         int NPTS = x.length;
         int FUNC = NDIMS;
         int ncalls = 0;
            ////// set up the starting simplex //////////////////
            double p[][]=new double[NPTS][NPTS]; // [row][col] = [whichvx][coord,FUNC]
            double z[]=new double[NDIMS];
            double best = 1E99;
            //////////////// initialize the funcvals ////////////////
            for (int i=0; i<NPTS; i++)
            {
              for (int j=0; j<NDIMS; j++)
                  {p[i][j] = x[i][j];}
              p[i][NDIMS] = fnelder.func(p[i]);
            }
            int iter=0;
            for (iter=1; iter<maxiteration; iter++)
            {
                /////////// identify lo, nhi, hi points //////////////
                int  ilo=0, ihi=0, inhi = -1; // -1 means missing
                double flo = p[0][FUNC];
                double fhi = flo;
                double pavg,sterr;
                for (int i=1; i<NPTS; i++)
                {
                   if (p[i][FUNC] < flo)
                     {flo=p[i][FUNC]; ilo=i;}
                   if (p[i][FUNC] > fhi)
                     {fhi=p[i][FUNC]; ihi=i;}
                }
                double fnhi = flo;
                inhi = ilo;
                for (int i=0; i<NPTS; i++)
                  if ((i != ihi) && (p[i][FUNC] > fnhi))
                    {fnhi=p[i][FUNC]; inhi=i;}
                ////////// exit criterion //////////////
                if ((iter % 4*NDIMS) == 0)
                {
                    //standart error of yi  should be less than set criteria (tolerance)
                    //calculate average (including highest point)
                    pavg=0;
                    for(int i=0;i<NPTS;i++)
                      pavg+=p[i][FUNC];
                    pavg/=NPTS;
                    double tot=0;
                    if(printlist!=0)
                    {  System.out.print(toString(iter,6));
                       for (int j=0; j<=NDIMS; j++)
                          { System.out.print(toString(p[ilo][j],25,12));}
                       System.out.println("");
                    }
                    for(int i=0;i<NPTS;i++)
                    { tot=(p[i][FUNC]-pavg)*(p[i][FUNC]-pavg);}
                    sterr=Math.sqrt(tot/NPTS);
                    //if(sterr < tolerance)
                      { for (int j=0; j<NDIMS; j++)
                           { z[j]=p[ilo][j];}
                        //break;
                      }
                    best = p[ilo][FUNC];
                }

                ///// compute ave[] vector excluding highest vertex //////

                double ave[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  ave[j] = 0;
                for (int i=0; i<NPTS; i++)
                  if (i != ihi)
                    for (int j=0; j<NDIMS; j++)
                       ave[j] += p[i][j];
                for (int j=0; j<NDIMS; j++)
                   ave[j] /= (NPTS-1);


                ///////// try reflect ////////////////

                double r[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  r[j] = 2*ave[j] - p[ihi][j];
                double fr = fnelder.func(r);

                if ((flo <= fr) && (fr < fnhi))  // in zone: accept
                {
                    for (int j=0; j<NDIMS; j++)
                      p[ihi][j] = r[j];
                    p[ihi][FUNC] = fr;
                    continue;
                }

                if (fr < flo)  //// below zone; try expand, else accept
                {
                    double e[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      e[j] = 3*ave[j] - 2*p[ihi][j];
                    double fe = fnelder.func(e);
                    if (fe < fr)
                    {
                       for (int j=0; j<NDIMS; j++)
                         p[ihi][j] = e[j];
                       p[ihi][FUNC] = fe;
                       continue;
                    }
                    else
                    {
                       for (int j=0; j<NDIMS; j++)
                         p[ihi][j] = r[j];
                       p[ihi][FUNC] = fr;
                       continue;
                    }
                }

                ///////////// above midzone, try contractions:

                if (fr < fhi)  /// try outside contraction
                {
                    double c[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      c[j] = 1.5*ave[j] - 0.5*p[ihi][j];
                    double fc = fnelder.func(c);
                    if (fc <= fr)
                    {
                        for (int j=0; j<NDIMS; j++)
                          p[ihi][j] = c[j];
                        p[ihi][FUNC] = fc;
                        continue;
                    }
                    else   /////// contract
                    {
                        for (int i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (int j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                        continue;
                    }
                }

                if (fr >= fhi)   /// over the top; try inside contraction
                {
                    double cc[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      cc[j] = 0.5*ave[j] + 0.5*p[ihi][j];
                    double fcc = fnelder.func(cc);
                    if (fcc < fhi)
                    {
                        for (int j=0; j<NDIMS; j++)
                          p[ihi][j] = cc[j];
                        p[ihi][FUNC] = fcc;
                        continue;
                    }
                    else    ///////// contract
                    {
                        for (int i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (int j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                    }
                }
            }
        return z;
}
      public static double[] nelder(f_xj fnelder,double x[][],double tolerance)
      {return nelder(fnelder,x,300,tolerance,0);}

      public static double[] nelder(f_xj fnelder,double x[][])
      {return nelder(fnelder,x,300,1.0e-15,0);}

      public static double[] nelder(f_xj fnelder,double x[][],int maxiteration,double tolerance)
      {return nelder(fnelder,x,maxiteration,tolerance,0);}

//derivatives of a function

public static double derivative(f_x f_deriv,double x)
{
// df/dx
// This method calculates derivatives of a simple function f(x)
// 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)
{
// dfi/dxj
// 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 f(x)dx
//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)
{
//integral f(x)dx
//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)
{
// x in f(x)=0 given also df(x)/dx
// 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)
{
// x in f(x)=0 , df(x)/dx is calcuated by derivative
// 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)
{
// x in f(x)=0 given also df(x)/dx d2f(x)/d2x is calculated by method
// 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)
{
// x in f(x)=0 given, df(x)/dx is calculated by derivative method
// d2f(x)/d2x is calculated by method
// 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)
{
// x in f(x)=0
// 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)
{
// x in f(x)=0 bracket xl and xu is given
//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.001;
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)
{
// x in f(x)=0 bracket xl and xu is given
//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)
{
// xi solution matrix in fj(xi)=0 set of non-linear equation
// derivative matrix dfj/dxi is required
//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)
{
// xi solution matrix in fj(xi)=0 set of non-linear equation
// derivative matrix dfj/dxi is calculated by derivative method
//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)
// data will be fit into
// y=f(x)=a0+a1*x+a2*x2+a3*x^3+....+an*x^n
// polynomial least square curve fitting
// variables xi,yi vector of data points
// degree of curve fitting
{
// data input variables:
//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;
}


// 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;
}


static String toString(int n, int w)
// converts an int to a string with given width.
{
    String s = Integer.toString(n);
    while (s.length() < w)
       s = " " + s;
    return s;
}

static String toString(int left)
// converts an int to a string with a constant width.
{
    return toString(left,6);
}

public static String toString(double left, int w, int d)
// converts a double to a string with given width and decimals.
{
	NumberFormat df=NumberFormat.getInstance(Locale.US);
	df.setMaximumFractionDigits(d);
    df.setMinimumFractionDigits(d);
    df.setGroupingUsed(false);
    String s = df.format(left);
    while (s.length() < w)
      s = " " + s;
    if (s.length() > w)
    {
        s = "";
        for (int i=0; i<w; i++)
          s = s + "-";
    }
    return s;
}

public static String toString(double left)
{// converts a double to a string with a constant width and constant decimals.
  return toString(left,15,10);}

//End of Class Numeric
}