rnatm

commit 599b0523dc28a4b2ba9550ff80e5222b9cefab2e
parent df5128e1d8db7d669924720c6aeb4281bba93aa7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue,  8 Nov 2022 20:40:19 +0100

Fix accuracy issues when getting k

To ensure that the radiative coefficients lie whithin k_min/k_max of
the acceleration structures, both calculations must be performed with
the same floating point precision, i.e. single precision

Diffstat:
M
src/rnatm_properties.c
 | 
20
++++++++++++++------

1 file changed, 14 insertions(+), 6 deletions(-)
diff --git a/src/rnatm_properties.c b/src/rnatm_properties.c
@@ -35,7 +35,7 @@
 #include <rsys/cstr.h>
 
 #define CUMULATIVE_SIZE_MAX 32
-typedef double (cumulative_T)[CUMULATIVE_SIZE_MAX];
+typedef float (cumulative_T)[CUMULATIVE_SIZE_MAX];
 
 /*******************************************************************************
  * Helper functions
@@ -981,7 +981,7 @@ compute_unnormalized_cumulative_radcoef
   size_t cpnt = RNATM_GAS;
   size_t naerosols = 0;
   size_t icumul = 0;
-  double k = 0;
+  float k = 0;
   res_T res = RES_OK;
   ASSERT(atm && (unsigned)radcoef < RNATM_RADCOEFS_COUNT__);
   ASSERT(cumulative && cumulative_sz);
@@ -1014,15 +1014,15 @@ compute_unnormalized_cumulative_radcoef
     res = rnatm_cell_get_radcoef(atm, &cell_args, &per_cell_k);
     if(res != RES_OK) goto error;
 
-    k += per_cell_k;
-    ASSERT(k <= k_max);
+    k += (float)per_cell_k;
     cumulative[icumul] = k;
     ++icumul;
 
   } while(++cpnt < naerosols);
   *cumulative_sz = icumul;
 
-  ASSERT(!icumul || k_min <= cumulative[icumul-1]);
+  ASSERT(!icumul || (float)k_min <= cumulative[icumul-1]);
+  ASSERT(!icumul || (float)k_max >= cumulative[icumul-1]);
 
 exit:
   return res;
@@ -1076,7 +1076,7 @@ rnatm_sample_component
    size_t* cpnt)
 {
   cumulative_T cumul;
-  double norm;
+  float norm;
   size_t cumul_sz;
   size_t i;
   res_T res = RES_OK;
@@ -1138,12 +1138,20 @@ rnatm_cell_get_radcoef
     /* Avoid precision issue when iterpolating the same value */
     k = vtx_k[0];
   } else {
+    float min_vtx_k;
+    float max_vtx_k;
+
     /* Calculate the radiative coefficient by linearly interpolating the
      * coefficients defined by tetrahedron vertex */
     k = vtx_k[0] * args->barycentric_coords[0]
       + vtx_k[1] * args->barycentric_coords[1]
       + vtx_k[2] * args->barycentric_coords[2]
       + vtx_k[3] * args->barycentric_coords[3];
+
+    /* Fix interpolation accuracy issues */
+    min_vtx_k = MMIN(MMIN(vtx_k[0], vtx_k[1]), MMIN(vtx_k[2], vtx_k[3]));
+    max_vtx_k = MMAX(MMAX(vtx_k[0], vtx_k[1]), MMAX(vtx_k[2], vtx_k[3]));
+    k = CLAMP((float)k, min_vtx_k, max_vtx_k);
   }
   ASSERT(k <= args->k_max);