Initial commit from Polytechnique Montreal

1 files changed, 139 insertions, 0 deletions

diff --git a/Utilib/src/zeroin.c b/Utilib/src/zeroin.c
new file mode 100644
index 0000000..dc250ff
--- /dev/null
+++ b/Utilib/src/zeroin.c
@@ -0,0 +1,139 @@
+/*
+freesteam - IAPWS-IF97 steam tables library
+Copyright (C) 2004-2009  John Pye
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*/
+
+#define FREESTEAM_BUILDING_LIB
+#include "zeroin.h"
+
+#include <math.h>
+#include <stdio.h>
+
+#ifndef DBL_EPSILON
+	#define DBL_EPSILON 2e-16
+#endif
+
+char zeroin_solve(ZeroInSubjectFunction *func, void *user_data, double lowerbound, double upperbound, double tol, double *solution, double *error){
+
+    double a, b, c;	/*  Abscissae, descr. see above. */
+    double fa;
+    double fb;
+    double fc;
+
+    a = lowerbound;
+    b = upperbound;
+    fa = (*func)(a,user_data);
+	    fb = (*func)(b,user_data);
+    c = a;
+    fc = fa;
+
+	if(fa == 0.){
+		*error = 0.; /* used by getError */
+		*solution = a;
+		return 0;
+	}
+
+    /*  Main iteration loop */
+
+    for (;;) {
+	    double prev_step = b - a;    /*  Distance from the last but one to the last approximation */
+	    double tol_act;              /*  Actual tolerance */
+	    double p;                    /*  Interpolation step is calculated in the form p/q; division */
+	    double q;                    /*  operations is delayed until the last moment */
+	    double new_step;             /*  Step at this iteration */
+
+	    if (fabs(fc) < fabs(fb)) {
+		    a = b;
+		    b = c;
+		    c = a;  /* Swap data for b to be the best approximation */
+		    fa = fb;
+		    fb = fc;
+		    fc = fa;
+	    }
+
+	    /* DBL_EPSILON is defined in math.h */
+	    tol_act = 2.0* DBL_EPSILON * fabs(b) + tol / 2.0;
+
+	    new_step = (c - b) / 2.0;
+
+	    if (fabs(new_step) <= tol_act || fb == 0.) {
+			*error = fb;
+			*solution = b;
+		    return 0;
+	    }
+	    /*  Decide if the interpolation can be tried */
+
+	    if (fabs(prev_step) >= tol_act   /* If prev_step was large enough and was in true direction, */
+	        && fabs(fa) > fabs(fb))      /* Interpolatiom may be tried */
+	    {
+		    register double t1, t2;
+		    double cb;
+
+		    cb = c - b;
+		    if (a == c) {
+			    /* If we have only two distinct points
+			       then only linear interpolation can be applied */
+			    t1 = fb / fa;
+			    p = cb * t1;
+			    q = 1.0 - t1;
+		    } else {
+			    /* Quadric inverse interpolation */
+
+			    q = fa / fc;
+			    t1 = fb / fc;
+			    t2 = fb / fa;
+			    p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0));
+			    q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
+		    }
+		    if (p > 0.) {
+			    /* p was calculated with the opposite sign; make p positive- */
+			    q = -q;	/* and assign possible minus to q */
+		    } else {
+			    p = -p;
+		    }
+
+		    if (p < (0.75 * cb * q - fabs(tol_act * q) / 2.0)
+		        && p < fabs(prev_step * q / 2.0)
+		       ) {
+			    /* If b+p/q falls in [b,c] and
+			       isn't too large it is accepted */
+			    new_step = p / q;
+		    }
+		    /* If p/q is too large then the bissection procedure can
+		       reduce [b,c] range to more extent */
+	    }
+
+	    if (fabs(new_step) < tol_act) {	/* Adjust the step to be not less */
+		    if (new_step > 0.)	/* than tolerance */
+			    new_step = tol_act;
+		    else
+			    new_step = -tol_act;
+	    }
+
+	    a = b;
+	    fa = fb;		/* Save the previous approx. */
+	    b += new_step;
+	    fb = (*func)(b,user_data);	/*  Do step to a new approxim. */
+
+	    if ((fb > 0. && fc > 0.)
+	        || (fb < 0. && fc < 0.)) {
+		    c = a;
+		    fc = fa;	/*  Adjust c for it to have a sign opposite to that of b */
+	    }
+    }
+    /* (((we never arrive here))) */
+}