star-wf

commit 7bf3e3b60734ff481fb2141d51f27ea9ce9de5e3
parent 84e35d25d679a5a7c2e5f7c9b5b76a9e9d347ee6
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed,  6 Mar 2024 17:21:05 +0100

Implement the swf_tabulation_inverse function

Diffstat:
M
src/swf_tabulation.c
 | 
37
+++++++++++++++++++++++++++++++++++++

1 file changed, 37 insertions(+), 0 deletions(-)
diff --git a/src/swf_tabulation.c b/src/swf_tabulation.c
@@ -16,6 +16,7 @@
 #include "swf.h"
 #include "swf_tabulation.h"
 
+#include <rsys/algorithm.h>
 #include <rsys/mem_allocator.h>
 
 /*******************************************************************************
@@ -30,6 +31,14 @@ release_tabulation(ref_T* ref)
   MEM_RM(tab->allocator, tab);
 }
 
+static FINLINE int
+cmp_item(const void* a, const void* b)
+{
+  const double f_a = *((const double*)a);
+  const double f_b = ((const struct item*)b)->f_x;
+  return f_a < f_b ? -1 : (f_a > f_b ? 1 : 0);
+}
+
 /*******************************************************************************
  * Exported symbols
  ******************************************************************************/
@@ -49,6 +58,34 @@ swf_tabulation_ref_put(struct swf_tabulation* tab)
   return RES_OK;
 }
 
+double
+swf_tabulation_inverse(const struct swf_tabulation* tab, const double y)
+{
+  const struct item* items = NULL;
+  const struct item* find = NULL;
+  double u = 0; /* Linear interpolation parameter */
+  double x = 0; /* Argument corresponding to input value y */
+  size_t nitems = 0;
+  ASSERT(tab && y >= 0 && y < 1);
+
+  items = darray_item_cdata_get(&tab->items);
+  nitems = darray_item_size_get(&tab->items);
+  ASSERT(nitems);
+
+  find = search_lower_bound(&y, items, nitems, sizeof(*items), cmp_item);
+  ASSERT(find && find->f_x >= y);
+
+  /* The input y correspond exactly to a tabulated argument */
+  if(find->f_x == y) return find->x;
+
+  /* Linear interpolation of tabulated arguments containing the argument
+   * corresponding to the input value */
+  ASSERT(find > items && find[-1].f_x < y);
+  u = (y - find[-1].f_x) / (find[0].f_x - find[-1].f_x);
+  x = u*(find[0].x - find[-1].x) + find[-1].x;
+  return x;
+}
+
 /*******************************************************************************
  * Local functions
  ******************************************************************************/