star-wf

commit 7207a3746c7a45734ab515b474aa5bd9076af351
parent c465bdbb86537b914653ccaa7d7073179f4c5c7b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Sun, 10 Mar 2024 11:45:24 +0100

Add quadratic interpolation to improve H inversion from tabs

Diffstat:
M
src/swf.h
 | 
6
++++++
M
src/swf_tabulation.c
 | 
38
+++++++++++++++++++++++++++++++++++---

2 files changed, 41 insertions(+), 3 deletions(-)
diff --git a/src/swf.h b/src/swf.h
@@ -31,6 +31,11 @@
   #define SWF(Func) swf_ ## Func
 #endif
 
+enum swf_interpolation {
+  SWF_LINEAR,
+  SWF_QUADRATIC
+};
+
 /* Forward declarations of external data types */
 struct mem_allocator;
 
@@ -93,6 +98,7 @@ swf_tabulation_ref_put
 SWF_API double /* x */
 swf_tabulation_inverse
   (const struct swf_tabulation* tab,
+   const enum swf_interpolation interpolation,
    const double y); /* in [0, 1[ */
 
 END_DECLS
diff --git a/src/swf_tabulation.c b/src/swf_tabulation.c
@@ -52,6 +52,25 @@ linear_interpolation
   return x;
 }
 
+/* H(x) = H(x0) + H'(x0)*(x-x0) + H"(x0)*(x-x0)^2
+ * x = x0 + (-H'(x0) +/-) sqrt(H'(x0)^2 - 2*H"(x0)*(H(x0) - H(x)))) / H"(x0)) */
+static INLINE double
+quadratic_extrapolation
+  (const double fx,
+   const struct item* a,
+   const struct item* b)
+{
+  double 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);
+  x0 = a->x + y0;
+  x1 = b->x + y1;
+
+  return x0 - a->x > b->x - x1 ? x1 : x0;
+}
+
 /*******************************************************************************
  * Exported symbols
  ******************************************************************************/
@@ -72,13 +91,17 @@ swf_tabulation_ref_put(struct swf_tabulation* tab)
 }
 
 double
-swf_tabulation_inverse(const struct swf_tabulation* tab, const double y)
+swf_tabulation_inverse
+  (const struct swf_tabulation* tab,
+   const enum swf_interpolation interpolation,
+   const double y)
 {
   const struct item* items = NULL;
   const struct item* find = NULL;
   double x = 0; /* Argument corresponding to input value y */
   size_t nitems = 0;
   ASSERT(tab && y >= 0 && y < 1);
+  ASSERT(interpolation == SWF_LINEAR || interpolation == SWF_QUADRATIC);
 
   items = darray_item_cdata_get(&tab->items);
   nitems = darray_item_size_get(&tab->items);
@@ -99,8 +122,17 @@ swf_tabulation_inverse(const struct swf_tabulation* tab, const double y)
 
   /* Linear interpolation of tabulated arguments containing the argument
    * corresponding to the input value */
-  ASSERT(find > items && find[-1].f_x < y);
-  x = linear_interpolation(y, &find[-1], &find[0]);
+  ASSERT(find > items && find[-1].fx < y);
+
+  switch(interpolation) {
+    case SWF_LINEAR:
+      x = linear_interpolation(y, &find[-1], &find[0]);
+      break;
+    case SWF_QUADRATIC:
+      x = quadratic_extrapolation(y, &find[-1], &find[0]);
+      break;
+    default: FATAL("Unreachable code\n"); break;
+  }
 
   return x;
 }