FUNCTION INSERT /* Insert a new point adjacent to the current minimum. * the point may lie to the left or right of x[index]. * Return TRUE if the new Y is a new minimum. */ int insert(double X, double Y){ assert(npoints > 2); assert((x[index-1]< X) && (X < x[index+1])); assert((0 < index) && (index < npoints-1)); assert(X != x[index]); if(X > x[index]) return insert_right(X, Y); else return insert_left(X, Y); } BISECTION // bisect the interval between elements k-1 and k STATE bisect_left(double (*f)(double), int k){ assert((k > 0) && (k < npoints-1) && (npoints > 2)); double X = (x[k] + x[k-1])/2; double Y = f(X); insert_left(X, Y); return PARA_FIT; } // bisect the interval between elements k and k+1 STATE bisect_right(double (*f)(double), int k){ assert((k > 0) && (k < npoints-1) && (npoints > 2)); double X = (x[k] + x[k+1])/2; double Y = f(X); insert_right(X, Y); return PARA_FIT; } // bisect the largest interval adjacent to the current minimum STATE bisect(double (*f)(double)){ assert((npoints > 2) && (index > 0) && (index < npoints-1)); double dx0 = x[index] - x[index-1]; double dx1 = x[index+1] - x[index]; if(dx0 > dx1) return bisect_left(f, index); else return bisect_right(f, index); } THE INTERPOLATION ROUTINES /* do a parabolic fit of the three points centered at index * If the fit produces a new minimum, repeat. * Else do a bisection. */ STATE para_fit(double (*f)(double), double &dx){ #define MAX_PARA 1 assert(npoints >= 3); assert((index > 0) && (index < npoints-1)); static int para_count = MAX_PARA; double X, Y, Yguess; double dx0 = x[index] - x[index-1]; double dx1 = x[index+1] - x[index]; double span = dx0 + dx1; double m0 = (y[index] - y[index-1])/dx0; double m1 = (y[index+1] - y[index])/dx1; double m = (y[index+1] - y[index-1])/span; double c = (m1 - m0)/span; double oldy = y[index]; // compute the estimated minimum X = (x[index-1] + x[index+1] - m/c)/2.0; dx = X - x[index]; Yguess = y[index-1] - c*X*X; Y = f(X); // make sure results are in range // don't insert if not if((X < x[index-1]) || (X > x[index+1])) return BISECT; if(!insert(X, Y)) return BISECT; if(!(--para_count)){ para_count = MAX_PARA; return BISECT; } if((Yguess < oldy) && (Y > oldy)) return FLIP; else return PARA_FIT; } // Flip the step to the opposite side of index STATE flip(double (*f)(double), double dx){ double X = x[index] - 2*dx; double Y = f(X); insert(X, Y); return PARA_FIT; } // do a "wide" parabolic fit of the three points x0, x2, and x4 void para_wide(double (*f)(double)){ assert(npoints ==5); double X, Yguess, Y; double m0 = (y[2] - y[0])/(x[2] - x[0]); double m1 = (y[4] - y[2])/(x[4] - x[2]); double m = (y[4] - y[0])/(x[4] - x[0]); double c = (m1 - m0)/(x[4] - x[0]); X = (x[0] + x[4] - m/c)/2.0; Yguess = y[0] - c*X*X; Y = f(X); insert(X, Y); } // Fit straight lines through the outer two pairs of points, // and find the intersection. // return the results in X and Y void line_fit(double (*f)(double)){ double X, Yguess, Y; assert(npoints ==5); double m0 = (y[1] - y[0])/(x[1] - x[0]); double m3 = (y[4] - y[3])/(x[4] - x[3]); cout << m0 << ' ' << m3 << endl; X = (m3*x[3] + m0*x[0] - (y[3]-y[0]))/(m3-m0); Yguess = y[0] + m0*(X-x[0]); Y = f(X); insert(X, Y); }