schiff

commit 045f481d02b694fc64747d4aea61628ddfa946b2
parent 83054901a0cecd3849f6184144ed82396c1b9401
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 24 Nov 2015 11:06:21 +0100

Improve the search of the optical properties

Use binary search rather than linear search to retrieve the optical
properties of a given wavelength.

Diffstat:
M
src/schiff_optical_properties.c
 | 
47
+++++++++++++++++++++++++++++------------------

1 file changed, 29 insertions(+), 18 deletions(-)
diff --git a/src/schiff_optical_properties.c b/src/schiff_optical_properties.c
@@ -29,6 +29,7 @@
 #include "schiff_optical_properties.h"
 #include "schiff_streamline.h"
 
+#include <rsys/algorithm.h>
 #include <rsys/cstr.h>
 #include <rsys/stretchy_array.h>
 
@@ -88,7 +89,7 @@ parse_optical_properties
   return RES_OK;
 }
 
-static int
+static FINLINE int
 cmp_properties(const void* op0, const void* op1)
 {
   const struct schiff_optical_properties* prop0 = op0;
@@ -96,6 +97,14 @@ cmp_properties(const void* op0, const void* op1)
   return prop0->W < prop1->W ? -1 : (prop0->W > prop1->W ? 1 : 0);
 }
 
+static FINLINE int
+cmp_wavelength_to_properties(const void* wavelength, const void* op)
+{
+  const double* W = wavelength;
+  const struct schiff_optical_properties* prop = op;
+  return *W < prop->W ? -1 : (*W > prop->W ? 1 : 0);
+}
+
 /*******************************************************************************
  * Local function
  ******************************************************************************/
@@ -173,29 +182,31 @@ schiff_optical_properties_fetch
    struct sschiff_material_properties* props)
 {
   const size_t nproperties = sa_size(properties);
+  struct schiff_optical_properties* find;
   double Ne, Kr, Nr;
-  size_t i;
   ASSERT(properties && props && wavelength > 0.0);
 
   /* Assume that properties are sorted in ascending order with respect to the
-   * wavelength. TODO use a binary search of the upper bound */
-  FOR_EACH(i, 0, nproperties) {
-    if(properties[i].W >= wavelength) break;
-  }
-
-  if(eq_eps(properties[i].W, wavelength, 1.e-12)
-  || i==0 || i==nproperties) { /* Clamp to wavelength */
-    Ne = properties[i].Ne;
-    Kr = properties[i].Kr;
-    Nr = properties[i].Nr;
+   * wavelength. */
+  find = search_lower_bound(&wavelength, properties, nproperties,
+    sizeof(struct schiff_optical_properties), cmp_wavelength_to_properties);
+  ASSERT(!find || find->W >= wavelength);
+
+  if(!find) /* Clamp to upper wavelength */
+    find = properties + (nproperties - 1);
+
+  if(find == properties || find == (properties + (nproperties - 1))
+  || eq_eps(find->W, wavelength, 1.e-12)) {
+    Ne = find[0].Ne;
+    Kr = find[0].Kr;
+    Nr = find[0].Nr;
   } else { /* Linear interpolation of optical properties */
-    const size_t i_prev = i - 1;
-    const double d = properties[i].W - properties[i_prev].W;
-    const double u = (wavelength - properties[i_prev].W) / d;
+    const double d = find[0].W - find[-1].W;
+    const double u = (wavelength - find[-1].W) / d;
     const double v = 1.0 - u;
-    Ne = u * properties[i].Ne + v * properties[i_prev].Ne;
-    Kr = u * properties[i].Kr + v * properties[i_prev].Kr;
-    Nr = u * properties[i].Nr + v * properties[i_prev].Nr;
+    Ne = u * find[0].Ne + v * find[-1].Ne;
+    Kr = u * find[0].Kr + v * find[-1].Kr;
+    Nr = u * find[0].Nr + v * find[-1].Nr;
   }
 
   props->medium_refractive_index = Ne;