schiff

commit 096dd4482775ec2698a415909ccf0e79523941c9
parent a2d8927486f97b3a6c9cba491701776c18b7e47a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 15 Mar 2016 16:40:32 +0100

Output the inverse cumulative phase function

Diffstat:
M
src/schiff.c
 | 
55
+++++++++++++++++++++++++++++++++++++++++++++++++++++++
M
src/schiff_args.h
 | 
2
++

2 files changed, 57 insertions(+), 0 deletions(-)
diff --git a/src/schiff.c b/src/schiff.c
@@ -122,7 +122,10 @@ run_integration
 {
   struct sschiff_device* sschiff = NULL;
   struct sschiff_estimator* estimator = NULL;
+  const double* angles = NULL;
+  double* thetas = NULL;
   size_t iwlen, nwlens;
+  size_t nangles;
   res_T res = RES_OK;
   ASSERT(args && sa_size(args->wavelengths) && rng && stream);
 
@@ -132,6 +135,12 @@ run_integration
     goto error;
   }
 
+  if(!sa_add(thetas, args->nangles_inv)) {
+    fprintf(stderr,
+"Couldn't allocate the list of inverse cumulative phase function angles.\n");
+    res = RES_MEM_ERR;
+    goto error;
+  }
 
   /* Invoke the Schiff integration */
   nwlens = sa_size(args->wavelengths);
@@ -161,10 +170,56 @@ run_integration
     val = &cross_section.average_projected_area;
     fprintf(stream, "%g %g\n", val->E, val->SE);
   }
+  fprintf(stream, "\n");
+
+  /* Print the phase functions estimations */
+  SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles));
+  FOR_EACH(iwlen, 0, nwlens) {
+    size_t iangle;
+    struct sschiff_state cdf, pf;
+    double n;
+
+    res = sschiff_estimator_get_limit_scattering_angle_index
+      (estimator, iwlen, &iangle);
+    if(res == RES_BAD_OP) { /* No valid limit angle => no phase function */
+      fprintf(stderr,
+        "The phase function couldn't be computed.\n"
+        "wavelength = %g microns\n", args->wavelengths[iwlen]);
+      fprintf(stream, "-1 -1 -1 -1 -1 -1 -1\n");
+      FOR_EACH(iangle, 0, args->nangles_inv) fprintf(stream, "-1\n");
+    } else {
+      ASSERT(res == RES_OK);
+      SSCHIFF(estimator_get_wide_scattering_angle_model_parameter
+        (estimator, iwlen, &n));
+      SSCHIFF(estimator_get_differential_cross_section
+        (estimator, iwlen, iangle, &pf));
+      SSCHIFF(estimator_get_differential_cross_section_cumulative
+        (estimator, iwlen, iangle, &cdf));
+
+      fprintf(stream, "%g %g %g %g %g %g %g\n",
+        args->wavelengths[iwlen], angles[iangle], cdf.E, cdf.SE, pf.E, pf.SE, n);
+
+      res = sschiff_estimator_inverse_cumulative_phase_function
+        (estimator, iwlen, thetas, args->nangles_inv);
+      if(res != RES_OK) {
+        fprintf(stderr,
+          "Couldn't inverse the cumulative phase function.\n"
+          "wavelength = %g microns, n = %g\n",
+          args->wavelengths[iwlen], n);
+        FOR_EACH(iangle, 0, args->nangles_inv) fprintf(stream, "-1\n");
+      } else {
+        FOR_EACH(iangle, 0, args->nangles_inv) {
+          fprintf(stream, "%g\n", thetas[iangle]);
+        }
+      }
+    }
+    fprintf(stream, "\n");
+  }
 
 exit:
   if(estimator) SSCHIFF(estimator_ref_put(estimator));
   if(sschiff) SSCHIFF(device_ref_put(sschiff));
+  if(thetas) sa_release(thetas);
   return res;
 error:
   goto exit;
diff --git a/src/schiff_args.h b/src/schiff_args.h
@@ -49,6 +49,7 @@ struct schiff_args {
   unsigned nrealisations; /* # realisation */
   unsigned ndirs; /* Number of directions to sample per realisation */
   unsigned nangles; /* Number of scattering angles */
+  unsigned nangles_inv; /* Number of angles in the the inverse phase function */
 
   unsigned nthreads; /* Hint on the number of thread to use */
 };
@@ -64,6 +65,7 @@ static const struct schiff_args SCHIFF_ARGS_NULL = {
   1000, /* # Sampled geometries */
   100, /* # Sampled directions per geometry */
   1000, /* # Scattering angles */
+  2000, /* # angles in the inv cumulative phase function */
   SSCHIFF_NTHREADS_DEFAULT
 };