star-wf

commit 07bd6358729991cca9640f6c793d1f880558ef77
parent e324a89de5800d68a10b84a5146c5dc059de7fc0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Sun, 10 Mar 2024 11:55:23 +0100

Add function swf_tabulation_max_x_error

Calculates the maximum tabulation error when interpolating x.
Currently, it seems to be buggy, since the maximum error with quadratic
interpolation would be greater than the maximum error with linear
interpolation.

Diffstat:
M
src/swf.h
 | 
6
++++++
M
src/swf_tabulation.c
 | 
32
++++++++++++++++++++++++++++++++
M
src/test_swf_H3d.c
 | 
4
++++

3 files changed, 42 insertions(+), 0 deletions(-)
diff --git a/src/swf.h b/src/swf.h
@@ -101,6 +101,12 @@ swf_tabulation_inverse
    const enum swf_interpolation interpolation,
    const double y); /* in [0, 1[ */
 
+/* Estimation of maximum tabulation error on X */
+SWF_API double
+swf_tabulation_max_x_error
+  (const struct swf_tabulation* tab,
+   const enum swf_interpolation interpolation);
+
 END_DECLS
 
 #endif /* SWF_H */
diff --git a/src/swf_tabulation.c b/src/swf_tabulation.c
@@ -137,6 +137,38 @@ swf_tabulation_inverse
   return x;
 }
 
+double
+swf_tabulation_max_x_error
+  (const struct swf_tabulation* tab,
+   const enum swf_interpolation interpolation)
+{
+  const struct item* items = NULL;
+  size_t i, n;
+  double max_error = 0;
+  ASSERT(tab);
+  ASSERT(interpolation == SWF_LINEAR || interpolation == SWF_QUADRATIC);
+
+  items = darray_item_cdata_get(&tab->items);
+  n = darray_item_size_get(&tab->items);
+
+  FOR_EACH(i, 1, n) {
+    const double delta = items[i].x - items[i-1].x;
+    double err = INF;
+
+    switch(interpolation) {
+      case SWF_LINEAR:
+        err = (items[i].d2fx * (delta*delta)/2.0) / items[i].dfx;
+        break;
+      case SWF_QUADRATIC:
+        err = (items[i].d3fx * (delta*delta*delta)/6.0) / items[i].dfx;
+        break;
+      default: FATAL("Unreachable code\n"); break;
+    }
+    max_error = MMAX(max_error, err);
+  }
+  return max_error;
+}
+
 /*******************************************************************************
  * Local functions
  ******************************************************************************/
diff --git a/src/test_swf_H3d.c b/src/test_swf_H3d.c
@@ -80,6 +80,10 @@ check_tabulation_inversion(void)
 
   CHK(swf_H3d_tabulate(&args, &tab) == RES_OK);
 
+  /*errl = swf_tabulation_max_x_error(tab, SWF_LINEAR);
+  errq = swf_tabulation_max_x_error(tab, SWF_QUADRATIC);
+  CHK(errq < errl);*/
+
   n = (size_t)ceil(1.0/delta_Hx);
 
   FOR_EACH(i, 0, n) {