star-wf

commit 70ebbd4523e821b7b944a37a4de682b81688e8cc
parent f4c603c2a4cd4782f56f3d911c403aec9d3fba9a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Sun, 10 Mar 2024 11:02:35 +0100

Add function swf_H3d_inverse

Calculates H(x)^-1 by dichotomy. This function will be used as a
reference to compare the inversion calculated from the tabulation of H.

Diffstat:
M
src/swf.h
 | 
4
++++
M
src/swf_H.c
 | 
23
+++++++++++++++++++++++

2 files changed, 27 insertions(+), 0 deletions(-)
diff --git a/src/swf.h b/src/swf.h
@@ -73,6 +73,10 @@ SWF_API double
 swf_H3d_eval
   (const double x);
 
+SWF_API double
+swf_H3d_inverse
+  (const double y); /* in [0, 1[ */
+
 SWF_API res_T
 swf_H3d_tabulate
   (const struct swf_H_tabulate_args* args,
diff --git a/src/swf_H.c b/src/swf_H.c
@@ -76,6 +76,29 @@ swf_H3d_eval(const double x)
   return 1.0 + 2.0 * sum;
 }
 
+double
+swf_H3d_inverse(const double y)
+{
+  const double epsilon = 1.0e-12;
+  double x_min = SWF_H_TABULATE_ARGS_DEFAULT.x_min;
+  double x_max = SWF_H_TABULATE_ARGS_DEFAULT.x_max;
+  ASSERT(y >= 0 && y < 1);
+
+  if(y == 0) return 0;
+
+  while(x_max - x_min > epsilon) {
+    const double x = (x_min + x_max) * 0.5;
+    const double f_x = swf_H3d_eval(x);
+
+    if(f_x < y) {
+      x_min = x;
+    } else {
+      x_max = x;
+    }
+  }
+  return (x_min + x_max) * 0.5;
+}
+
 res_T
 swf_H3d_tabulate
   (const struct swf_H_tabulate_args* args,