star-wf

commit 7479559cea358a9f29145a31b078a2ffec8a5ef3
parent 09198e97681e4bf3ba62851a7f2e2252f6af643e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 11 Mar 2024 16:26:22 +0100

Detect unexpected numerical errors in quadratic extrapolation

Check that a solution exists (i.e. that the tabulation is not too
coarse) and that the extrapolated x lies within the tabulated interval
used to caculate it.

Diffstat:
M
src/swf_tabulation.c
 | 
16
++++++++++++----

1 file changed, 12 insertions(+), 4 deletions(-)
diff --git a/src/swf_tabulation.c b/src/swf_tabulation.c
@@ -60,15 +60,23 @@ quadratic_extrapolation
    const struct item* a,
    const struct item* b)
 {
-  double x0, x1,  y0, y1;
+  double x;
+  double d0, d1, x0, x1, y0, y1;
   ASSERT(a && b);
 
-  y0 = (-a->dfx + sqrt(a->dfx*a->dfx - 2*a->d2fx*(a->fx - fx))) / (a->d2fx);
-  y1 = (-b->dfx + sqrt(b->dfx*b->dfx - 2*b->d2fx*(b->fx - fx))) / (b->d2fx);
+  d0 = a->dfx*a->dfx - 2*a->d2fx*(a->fx - fx);
+  d1 = b->dfx*b->dfx - 2*b->d2fx*(b->fx - fx);
+  ASSERT(d0 && d1);
+
+  y0 = (-a->dfx + sqrt(d0)) / (a->d2fx);
+  y1 = (-b->dfx + sqrt(d1)) / (b->d2fx);
   x0 = a->x + y0;
   x1 = b->x + y1;
 
-  return x0 - a->x > b->x - x1 ? x1 : x0;
+  x = x0 - a->x > b->x - x1 ? x1 : x0;
+  ASSERT(a->x <= x && x <= b->x);
+
+  return x;
 }
 
 /*******************************************************************************