private double secondPartialDerivativeRidr(FittingFunction f,int i,int j,double[] args,double delt){
  double[] arg=new double[args.length];
  double[][] a=new double[NTAB][NTAB];
  double hh=delt;
  double errt;
  double ans=0.0;
  double fac;
  System.arraycopy(args,0,arg,0,args.length);
  double center=f.evaluate(arg);
  arg[i]+=delt;
  arg[j]+=delt;
  double ff1=f.evaluate(arg);
  arg[j]-=2 * delt;
  double ff2=f.evaluate(arg);
  arg[i]-=2 * delt;
  arg[j]+=2 * delt;
  double ff3=f.evaluate(arg);
  arg[j]-=2 * delt;
  double ff4=f.evaluate(arg);
  if (Double.isNaN(ff1)) {
    ff1=center;
  }
  if (Double.isNaN(ff2)) {
    ff2=center;
  }
  if (Double.isNaN(ff3)) {
    ff3=center;
  }
  if (Double.isNaN(ff4)) {
    ff4=center;
  }
  a[0][0]=(ff1 - ff2 - ff3 + ff4) / (4.0 * delt * delt);
  double err=BIG;
  for (int ii=1; ii < NTAB; ii++) {
    hh/=CON;
    System.arraycopy(args,0,arg,0,args.length);
    arg[i]+=hh;
    arg[j]+=hh;
    ff1=f.evaluate(arg);
    arg[j]-=2 * hh;
    ff2=f.evaluate(arg);
    arg[i]-=2 * hh;
    arg[j]+=2 * hh;
    ff3=f.evaluate(arg);
    arg[j]-=2 * hh;
    ff4=f.evaluate(arg);
    if (Double.isNaN(ff1)) {
      ff1=center;
    }
    if (Double.isNaN(ff2)) {
      ff2=center;
    }
    if (Double.isNaN(ff3)) {
      ff3=center;
    }
    if (Double.isNaN(ff4)) {
      ff4=center;
    }
    a[0][ii]=(ff1 - ff2 - ff3 + ff4) / (4.0 * hh * hh);
    fac=CON2;
    for (int jj=1; jj < ii; jj++) {
      a[jj][ii]=(a[jj - 1][ii] * fac - a[jj - 1][ii - 1]) / (fac - 1.0);
      fac=CON2 * fac;
      errt=Math.max(Math.abs(a[jj][ii] - a[jj - 1][ii]),Math.abs(a[jj][ii] - a[jj - 1][ii - 1]));
      if (errt < err) {
        err=errt;
        ans=a[jj][ii];
      }
    }
    if (Math.abs(a[ii][ii] - a[ii - 1][ii - 1]) >= SAFE * err) {
      break;
    }
  }
  return ans;
}

private double secondPartialDerivativeRidr ( FittingFunction f , int i , int j , double [ ] args , double delt ) { double [ ] arg = new double [ args . length ] ; double [ ] [ ] a = new double [ NTAB ] [ NTAB ] ; double hh = delt ; double errt ; double ans = 0.0 ; double fac ; System . arraycopy ( args , 0 , arg , 0 , args . length ) ; double center = f . evaluate ( arg ) ; arg [ i ] += delt ; arg [ j ] += delt ; double ff1 = f . evaluate ( arg ) ; arg [ j ] -= 2 * delt ; double ff2 = f . evaluate ( arg ) ; arg [ i ] -= 2 * delt ; arg [ j ] += 2 * delt ; double ff3 = f . evaluate ( arg ) ; arg [ j ] -= 2 * delt ; double ff4 = f . evaluate ( arg ) ; if ( Double . isNaN ( ff1 ) ) { ff1 = center ; } if ( Double . isNaN ( ff2 ) ) { ff2 = center ; } if ( Double . isNaN ( ff3 ) ) { ff3 = center ; } if ( Double . isNaN ( ff4 ) ) { ff4 = center ; } a [ 0 ] [ 0 ] = ( ff1 - ff2 - ff3 + ff4 ) / ( 4.0 * delt * delt ) ; double err = BIG ; for ( int ii = 1 ; ii < NTAB ; ii ++ ) { hh /= CON ; System . arraycopy ( args , 0 , arg , 0 , args . length ) ; arg [ i ] += hh ; arg [ j ] += hh ; ff1 = f . evaluate ( arg ) ; arg [ j ] -= 2 * hh ; ff2 = f . evaluate ( arg ) ; arg [ i ] -= 2 * hh ; arg [ j ] += 2 * hh ; ff3 = f . evaluate ( arg ) ; arg [ j ] -= 2 * hh ; ff4 = f . evaluate ( arg ) ; if ( Double . isNaN ( ff1 ) ) { ff1 = center ; } if ( Double . isNaN ( ff2 ) ) { ff2 = center ; } if ( Double . isNaN ( ff3 ) ) { ff3 = center ; } if ( Double . isNaN ( ff4 ) ) { ff4 = center ; } a [ 0 ] [ ii ] = ( ff1 - ff2 - ff3 + ff4 ) / ( 4.0 * hh * hh ) ; fac = CON2 ; for ( int jj = 1 ; jj < ii ; jj ++ ) { a [ jj ] [ ii ] = ( a [ jj - 1 ] [ ii ] * fac - a [ jj - 1 ] [ ii - 1 ] ) / ( fac - 1.0 ) ; fac = CON2 * fac ; errt = Math . max ( Math . abs ( a [ jj ] [ ii ] - a [ jj - 1 ] [ ii ] ) , Math . abs ( a [ jj ] [ ii ] - a [ jj - 1 ] [ ii - 1 ] ) ) ; if ( errt < err ) { err = errt ; ans = a [ jj ] [ ii ] ; } } if ( Math . abs ( a [ ii ] [ ii ] - a [ ii - 1 ] [ ii - 1 ] ) >= SAFE * err ) { break ; } } return ans ; }
this method implements ridder ' s algorithm for computing the second order partial derivatives . it is a translation of the c program in section 5 . 7 of numerical recipes in c . it is more robust than the above method in that it searches for a perferred value of delt . but based on our experience to date with sem fitting functions , the above method seems to be adequately accurate and faster that this one .
