star-schiff

commit 831d21d76495d54c465f5c3fcdcc47c679b93af1
parent c150b1f843df36f2520f9edc13f8b430b8a94972
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 21 Oct 2015 11:10:49 +0200

Update the cylinder Schiff estimation test

Take into account the standard error of the reference value in the
comparison of the estimated result.

Diffstat:
M
src/test_sschiff_estimator_cylinder.c
 | 
89
++++++++++++++++++++++---------------------------------------------------------

1 file changed, 24 insertions(+), 65 deletions(-)
diff --git a/src/test_sschiff_estimator_cylinder.c b/src/test_sschiff_estimator_cylinder.c
@@ -38,7 +38,8 @@
 
 struct result {
   double wavelength;
-  double extinction_cross_section;
+  double extinction_E; /* Expected value */
+  double extinction_SE; /* Standard error */
 };
 
 struct optical_properties { double wavelength, n_r, kappa_r, n_medium; };
@@ -812,62 +813,13 @@ main(int argc, char** argv)
   struct ssp_rng* rng;
   struct time t0, t1;
   struct result results[] = {
-    { 4.00e-07, 2.76750904574817e-12 },
-    { 4.10e-07, 2.73672907058880e-12 },
-    { 4.20e-07, 2.66637662867280e-12 },
-    { 4.30e-07, 2.59077868096684e-12 },
-    { 4.40e-07, 2.54198659470604e-12 },
-    { 4.50e-07, 2.50099694990826e-12 },
-    { 4.60e-07, 2.46942367528726e-12 },
-    { 4.70e-07, 2.44558265209247e-12 },
-    { 4.80e-07, 2.41111486015226e-12 },
-    { 4.90e-07, 2.37268038342598e-12 },
-    { 5.00e-07, 2.32782435081181e-12 },
-    { 5.10e-07, 2.31478533603249e-12 },
-    { 5.20e-07, 2.30848723160716e-12 },
-    { 5.30e-07, 2.26268600507108e-12 },
-    { 5.40e-07, 2.21637927974485e-12 },
-    { 5.50e-07, 2.21022013899174e-12 },
-    { 5.60e-07, 2.18986515715339e-12 },
-    { 5.70e-07, 2.13806095885549e-12 },
-    { 5.80e-07, 2.06975243849105e-12 },
-    { 5.90e-07, 2.01854621247124e-12 },
-    { 6.00e-07, 1.96661661205101e-12 },
-    { 6.10e-07, 1.90708640327711e-12 },
-    { 6.20e-07, 1.85044620857265e-12 },
-    { 6.30e-07, 1.79722475451741e-12 },
-    { 6.40e-07, 1.75148807605901e-12 },
-    { 6.50e-07, 1.71149696366194e-12 },
-    { 6.60e-07, 1.67387968010696e-12 },
-    { 6.70e-07, 1.63918775517145e-12 },
-    { 6.80e-07, 1.60549326749567e-12 },
-    { 6.90e-07, 1.57349750025259e-12 },
-    { 7.00e-07, 1.54229472112306e-12 },
-    { 7.10e-07, 1.51255220254738e-12 },
-    { 7.20e-07, 1.48275984511749e-12 },
-    { 7.30e-07, 1.45374742411889e-12 },
-    { 7.40e-07, 1.42499165124167e-12 },
-    { 7.50e-07, 1.39877652964399e-12 },
-    { 7.60e-07, 1.37346471683876e-12 },
-    { 7.70e-07, 1.34367558859022e-12 },
-    { 7.80e-07, 1.33277978493000e-12 },
-    { 7.90e-07, 1.32305511913364e-12 },
-    { 8.00e-07, 1.31249380149019e-12 },
-    { 8.10e-07, 1.29863360991060e-12 },
-    { 8.20e-07, 1.26716950339541e-12 },
-    { 8.30e-07, 1.23935510053982e-12 },
-    { 8.40e-07, 1.23737221612295e-12 },
-    { 8.50e-07, 1.25066511976589e-12 },
-    { 8.60e-07, 1.28717131701198e-12 },
-    { 8.70e-07, 1.34490237420400e-12 },
-    { 8.80e-07, 1.39983308298635e-12 },
-    { 8.90e-07, 1.41749581570064e-12 },
-    { 9.00e-07, 1.35842722128322e-12 },
-    { 9.10e-07, 1.25760593510938e-12 },
-    { 9.20e-07, 1.17176942304411e-12 },
-    { 9.30e-07, 1.12091282607288e-12 },
-    { 9.40e-07, 1.08845962290538e-12 },
-    { 9.50e-07, 1.06375988557876e-12 }
+    { 4.00E-7, 2.76750904574817e-12, 2.34364058801497e-15 },
+    { 4.50e-7, 2.50099694990826e-12, 2.20010054373943e-15 },
+    { 5.00e-7, 2.32782435081181e-12, 2.05630225764867e-15 },
+    { 6.00e-7, 1.96661661205101e-12, 1.82091931178789e-15 },
+    { 6.50e-7, 1.71149696366194e-12, 1.63845385406345e-15 },
+    { 7.00e-7, 1.54229472112306e-12, 1.49426946826458e-15 },
+    { 7.50e-7, 1.39877652964399e-12, 1.36729344819050e-15 }
   };
   const size_t nresults = sizeof(results) / sizeof(struct result);
   size_t i;
@@ -888,10 +840,10 @@ main(int argc, char** argv)
 
   FOR_EACH(i, 0, nresults) {
     const size_t nscatt_angles = 1;
-    const size_t ngeoms = 300;
-    const size_t ndirs = 30;
+    const size_t ngeoms = 1000;
+    const size_t ndirs = 10;
     const double wavelength = results[i].wavelength;
-    double dst;
+    double interval0[2], interval1[2], interval2[2];
 
     time_current(&t0);
     CHECK(sschiff_integrate(dev, rng, &distrib, wavelength, nscatt_angles,
@@ -902,19 +854,26 @@ main(int argc, char** argv)
 
     CHECK(sschiff_estimator_get_extinction_cross_section(estimator, &status), RES_OK);
     CHECK(status.nsteps, ngeoms);
-    dst = status.E - results[i].extinction_cross_section;
+
+    interval0[0] = results[i].extinction_E - 3*results[i].extinction_SE;
+    interval0[1] = results[i].extinction_E + 3*results[i].extinction_SE;
+    interval1[0] = status.E - 3*status.SE;
+    interval1[1] = status.E + 3*status.SE;
+    interval2[0] = MMAX(interval0[0], interval1[0]);
+    interval2[1] = MMIN(interval0[1], interval1[1]);
 
     printf(
 "%2.1d - Wavelength: %7.2g - "
-"Extinction cross section = %12.7g ~ %12.7g +/- %12.7g (%6.2g) - "
+"Extinction cross section ~ %12.7g +/- %12.7g ~ %12.7g +/- %12.7g (%6.2g) - "
 "%s\n",
       (int)i,
       results[i].wavelength,
-      results[i].extinction_cross_section,
-      status.E, status.SE, dst / status.SE,
+      results[i].extinction_E,
+      results[i].extinction_SE,
+      status.E, status.SE, interval2[1] - interval2[0],
       buf);
 
-    CHECK(eq_eps(results[i].extinction_cross_section, status.E, 3*status.SE), 1);
+    CHECK(interval2[0] <= interval2[1], 1);
     CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
   }