commit 747d6686186a368619dbf1313e0e9bf047b6b6a0
parent 900372082fa4dc19fd75d7d76078cf57df796dc6
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 29 Mar 2016 12:11:12 +0200
Major update of the schiff output
Write the cross sections, the phase function descriptors, the phase
functions, the cumulative phase functions and the inverse cumulative
phase functions. The [[inverse] cumulative] phase functions are
preceded by their associated value (an angle or a probability for the
inverse cumulative).
Diffstat:
2 files changed, 208 insertions(+), 56 deletions(-)
diff --git a/src/schiff.c b/src/schiff.c
@@ -100,53 +100,22 @@ error:
goto exit;
}
-static res_T
-run_integration
- (const struct schiff_args* args,
- struct sschiff_geometry_distribution* distrib,
- struct ssp_rng* rng,
- FILE* stream)
+static void
+write_cross_sections
+ (FILE* stream,
+ struct sschiff_estimator* estimator,
+ const double* wlens, /* List of wavelengths in microns */
+ const size_t nwlens) /* #wavelengths */
{
- 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);
-
- res = sschiff_device_create(NULL, NULL, args->nthreads, 1, NULL, &sschiff);
- if(res != RES_OK) {
- fprintf(stderr, "Couldn't create the Star Schiff device.\n");
- 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;
- }
+ size_t iwlen;
+ ASSERT(stream && estimator && wlens && nwlens);
- /* Invoke the Schiff integration */
- nwlens = sa_size(args->wavelengths);
- res = sschiff_integrate(sschiff, rng, distrib, args->wavelengths, nwlens,
- sschiff_uniform_scattering_angles, args->nangles, args->nrealisations,
- args->ndirs, &estimator);
- if(res != RES_OK) {
- fprintf(stderr, "Schiff integration error.\n");
- goto error;
- }
-
- /* Print the estimated cross sections */
FOR_EACH(iwlen, 0, nwlens) {
const struct sschiff_state* val;
struct sschiff_cross_section cross_section;
SSCHIFF(estimator_get_cross_section(estimator, iwlen, &cross_section));
-
- fprintf(stream, "%g ", args->wavelengths[iwlen]);
+ fprintf(stream, "%g ", wlens[iwlen]);
val = &cross_section.extinction;
fprintf(stream, "%g %g ", val->E, val->SE);
@@ -158,9 +127,24 @@ run_integration
fprintf(stream, "%g %g\n", val->E, val->SE);
}
fprintf(stream, "\n");
+}
+
+static void
+write_phase_function_descriptors
+ (FILE* stream,
+ struct sschiff_estimator* estimator,
+ const double* wlens, /* List of wavelengths in microns */
+ const size_t nwlens, /* #wavelenghts */
+ const size_t nangles_inv) /* Number of inversed cumulative entries */
+{
+ size_t iwlen;
+ const double* angles;
+ size_t nangles;
+ res_T res = RES_OK;
+ ASSERT(stream && estimator && wlens && nwlens && nangles_inv > 1);
- /* 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;
@@ -168,13 +152,16 @@ run_integration
res = sschiff_estimator_get_limit_scattering_angle_index
(estimator, iwlen, &iangle);
- if(res == RES_BAD_OP) { /* No valid limit angle => no phase function */
+
+ if(res == RES_BAD_OP) {
+ /* No valid limit angle => no phase function. Print -1 for all the data */
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");
+ "wavelength = %g microns\n", wlens[iwlen]);
+ fprintf(stream, "-1 -1 -1 -1 -1 -1 -1 -1 -1\n");
+
} else {
+ /* Retrieve the phase function descriptor data */
ASSERT(res == RES_OK);
SSCHIFF(estimator_get_wide_scattering_angle_model_parameter
(estimator, iwlen, &n));
@@ -183,30 +170,184 @@ run_integration
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], pf.E, pf.SE, cdf.E, cdf.SE, n);
+ /* Write the phase function descriptor */
+ fprintf(stream, "%g %g %g %g %g %g %g %lu %lu\n",
+ wlens[iwlen], angles[iangle], pf.E, pf.SE, cdf.E, cdf.SE, n,
+ (unsigned long)nangles,
+ (unsigned long)nangles_inv);
+ }
+ }
+ fprintf(stream, "\n"); /* Empty line */
+}
+
+static void
+write_phase_functions
+ (FILE* stream,
+ const struct sschiff_estimator* estimator,
+ const double* wlens,
+ const size_t nwlens)
+{
+ size_t iwlen;
+ const double* angles;
+ const struct sschiff_state* phase_funcs;
+ size_t iangle, nangles;
+ ASSERT(stream && estimator && wlens && nwlens);
+ (void)wlens;
+
+ SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles));
+ FOR_EACH(iwlen, 0, nwlens) {
+ const res_T res = sschiff_estimator_get_phase_function
+ (estimator, iwlen, &phase_funcs);
+ if(res == RES_BAD_OP) {
+ /* The phase function couldn't be computed. Write default data. */
+ FOR_EACH(iangle, 0, nangles) {
+ fprintf(stream, "-1 -1 -1\n");
+ }
+ } else {
+ ASSERT(res == RES_OK);
+ FOR_EACH(iangle, 0, nangles) {
+ fprintf(stream, "%g %g %g\n",
+ angles[iangle],
+ phase_funcs[iangle].E,
+ phase_funcs[iangle].SE);
+ }
+ }
+ fprintf(stream, "\n");
+ }
+}
+
+static void
+write_phase_functions_cum
+ (FILE* stream,
+ struct sschiff_estimator* estimator,
+ const double* wlens,
+ const size_t nwlens)
+{
+ size_t iwlen;
+ const double* angles;
+ const struct sschiff_state* phase_funcs_cum;
+ size_t iangle, nangles;
+ ASSERT(stream && estimator && wlens && nwlens);
+ (void)wlens;
+
+ SSCHIFF(estimator_get_scattering_angles(estimator, &angles, &nangles));
+ FOR_EACH(iwlen, 0, nwlens) {
+ const res_T res = sschiff_estimator_get_phase_function_cumulative
+ (estimator, iwlen, &phase_funcs_cum);
+ if(res == RES_BAD_OP) {
+ /* The cumulative phase func couldn't be computed. Write default data */
+ FOR_EACH(iangle, 0, nangles) {
+ fprintf(stream, "-1 -1 -1\n");
+ }
+ } else {
+ ASSERT(res == RES_OK);
+ FOR_EACH(iangle, 0, nangles) {
+ fprintf(stream, "%g %g %g\n",
+ angles[iangle],
+ phase_funcs_cum[iangle].E,
+ phase_funcs_cum[iangle].SE);
+ }
+ }
+ fprintf(stream, "\n");
+ }
+}
+
+static void
+write_phase_functions_invcum
+ (FILE* stream,
+ struct sschiff_estimator* estimator,
+ const double* wlens,
+ const size_t nwlens,
+ const size_t nangles_inv)
+{
+ double* invcum = NULL;
+ size_t iangle;
+ size_t iwlen;
+ double step;
+ ASSERT(stream && estimator && wlens && nwlens && nangles_inv > 1);
+
+ step = 1.0 / (double)(nangles_inv - 1);
- res = sschiff_estimator_inverse_cumulative_phase_function
- (estimator, iwlen, thetas, args->nangles_inv);
+ if(!sa_add(invcum, nangles_inv)) {
+ fprintf(stderr,
+ "Couldn't allocate the list of inverse cumulative phase function angles.\n");
+
+ /* Print the default -1 value */
+ FOR_EACH(iwlen, 0, nwlens) {
+ FOR_EACH(iangle, 0, nangles_inv) {
+ fprintf(stream, "-1 -1\n");
+ }
+ fprintf(stream, "\n");
+ }
+ } else {
+ FOR_EACH(iwlen, 0, nwlens) {
+ const res_T res = sschiff_estimator_inverse_cumulative_phase_function
+ (estimator, iwlen, invcum, nangles_inv);
if(res != RES_OK) {
+ /* The inverse cumulative cannot be computed. Write -1 */
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");
+ "wavelength = %g microns\n", wlens[iwlen]);
+ FOR_EACH(iangle, 0, nangles_inv) {
+ fprintf(stream, "-1 -1\n");
+ }
} else {
- FOR_EACH(iangle, 0, args->nangles_inv) {
- fprintf(stream, "%g\n", thetas[iangle]);
+ /* Write the inverse cumulative values */
+ FOR_EACH(iangle, 0, nangles_inv-1) {
+ fprintf(stream, "%g %g\n", (double)iangle*step, invcum[iangle]);
}
+ /* Handle precision issues */
+ fprintf(stream, "1 %g\n", invcum[nangles_inv-1]);
}
+ fprintf(stream, "\n"); /* Empty line */
}
- fprintf(stream, "\n");
}
+ if(invcum)
+ sa_release(invcum);
+}
+
+static res_T
+run_integration
+ (const struct schiff_args* args,
+ struct sschiff_geometry_distribution* distrib,
+ struct ssp_rng* rng,
+ FILE* stream)
+{
+ struct sschiff_device* sschiff = NULL;
+ struct sschiff_estimator* estimator = NULL;
+ size_t nwlens;
+ res_T res = RES_OK;
+ ASSERT(args && sa_size(args->wavelengths) && rng && stream);
+
+ res = sschiff_device_create(NULL, NULL, args->nthreads, 1, NULL, &sschiff);
+ if(res != RES_OK) {
+ fprintf(stderr, "Couldn't create the Star Schiff device.\n");
+ goto error;
+ }
+
+ /* Invoke the Schiff integration */
+ nwlens = sa_size(args->wavelengths);
+ res = sschiff_integrate(sschiff, rng, distrib, args->wavelengths, nwlens,
+ sschiff_uniform_scattering_angles, args->nangles, args->nrealisations,
+ args->ndirs, &estimator);
+ if(res != RES_OK) {
+ fprintf(stderr, "Schiff integration error.\n");
+ goto error;
+ }
+
+ /* Write results */
+ write_cross_sections(stream, estimator, args->wavelengths, nwlens);
+ write_phase_function_descriptors
+ (stream, estimator, args->wavelengths, nwlens, args->nangles_inv);
+ write_phase_functions(stream, estimator, args->wavelengths, nwlens);
+ write_phase_functions_cum(stream, estimator, args->wavelengths, nwlens);
+ write_phase_functions_invcum
+ (stream, estimator, args->wavelengths, nwlens, args->nangles_inv);
+
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.c b/src/schiff_args.c
@@ -86,6 +86,14 @@ print_help(const char* binary)
" integrate.\n\n");
}
+static int
+cmp_double(const void* op0, const void* op1)
+{
+ const double* a = op0;
+ const double* b = op1;
+ return *a < *b ? -1 : (*a > *b ? 1 : 0);
+}
+
static INLINE void
param_release(struct schiff_param* param)
{
@@ -1466,6 +1474,9 @@ schiff_args_init
res = RES_BAD_ARG;
goto error;
}
+ /* Sort the submitted wavelengths in ascending order */
+ qsort(args->wavelengths, sa_size(args->wavelengths),
+ sizeof(args->wavelengths[0]), cmp_double);
if(args->characteristic_length < 0.0) {
fprintf(stderr,
