commit 3ff881a2accdc907ad9b5d7de412e8ab88d366a0
parent 7e84fcd1788b967a6af17d953ed911a10febb949
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 2 Mar 2016 14:24:56 +0100
Updates the Star-Schiff API
The integration ensures than the submitted wavelengths are sorted in
ascending order.
Use the index of the wavelength rather than its real value to retrieve
its estimated data.
Replace the function that returns the number of wavelengths by
sschiff_estimator_get_wavelengths that optionally returns the number
and/or the list of wavelengths.
Add a function that returns the list of scattering angles.
Add functions to retrieve the estimated phase function [cumulative].
Diffstat:
5 files changed, 150 insertions(+), 67 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -148,7 +148,8 @@ sschiff_integrate
(struct sschiff_device* dev,
struct ssp_rng* rng,
struct sschiff_geometry_distribution* distribution,
- const double* wavelengths, /* list of wavelengths to estimate in micron */
+ /* Wavelengths to estimate in micron. Must be sorted in ascending order. */
+ const double* wavelengths,
const size_t wavelengths_count, /* # wavelengths to estimate */
const sschiff_scattering_angles_distribution_T angles, /* angles distrib */
const size_t scattering_angles_count, /* # scattering angles. Must be >= 3 */
@@ -164,11 +165,19 @@ SSCHIFF_API res_T
sschiff_estimator_ref_put
(struct sschiff_estimator* estimator);
-/* Return the number of estimated wavelength */
+/* Return the list of estimated wavelengths. */
SSCHIFF_API res_T
-sschiff_estimator_get_wavelengths_count
+sschiff_estimator_get_wavelengths
(const struct sschiff_estimator* estimator,
- size_t* count);
+ const double** wavelengths, /* May be NULL, i.e. do not return this list */
+ size_t* count); /* May be NULL, i.e. do not return the # wavelengths */
+
+/* Return the list of scattering angles. */
+SSCHIFF_API res_T
+sschiff_estimator_get_scattering_angles
+ (const struct sschiff_estimator* estimator,
+ const double** angles, /* May be NULL, i.e. do not return the angles */
+ size_t* count); /* May be NULL, i.e. do not return the # scattering angles */
/* Return the number of Monte Carlo realisations */
SSCHIFF_API res_T
@@ -176,27 +185,51 @@ sschiff_estimator_get_realisations_count
(const struct sschiff_estimator* estimator,
size_t* count);
-/* Retrieve the estimation state of a given wavelength */
+/* Retrieve the estimation state of a given wavelength. */
SSCHIFF_API res_T
-sschiff_estimator_get_wavelength_cross_section
+sschiff_estimator_get_cross_section
(const struct sschiff_estimator* estimator,
- const double wavelength, /* In micron */
+ const size_t wavelength_index, /* Id of the wavelengths */
struct sschiff_cross_section* cross_section);
/* Retrieve the estimated cross sections of all wavelengths. The length of
* `cross_sections' must be at least equal to the count of integrated
- * wavelengths. One can use the sschiff_estimator_get_wavelengths_count
- * function to retrieve this information. */
+ * wavelengths. One can use the sschiff_estimator_get_wavelengths function to
+ * retrieve this information. */
SSCHIFF_API res_T
sschiff_estimator_get_cross_sections
(const struct sschiff_estimator* estimator,
struct sschiff_cross_section cross_sections[]);
+/* Retrieve a pointer onto the estimated phase function for the scattering
+ * angles of the estimator. Use the sschiff_estimator_get_scattering_angles
+ * function to retrieve these scattering angles. */
+SSCHIFF_API res_T
+sschiff_estimator_get_phase_function
+ (const struct sschiff_estimator* estimator,
+ const size_t wavelength_index,
+ const struct sschiff_state** states);
+
+/* Retrieve a pointer onto the estimated phase function cumulative for the
+ * scattering angles of the estimator. Use the
+ * sschiff_estimator_get_scattering_angles function to retrieve these
+ * scattering angles */
+SSCHIFF_API res_T
+sschiff_estimator_get_phase_function_cumulative
+ (const struct sschiff_estimator* estimator,
+ const size_t wavelength_index,
+ const struct sschiff_state** states);
+
+/* Compute the inverse cumulative of the estimated phase function. The inverse
+ * cumulative is computed for a set of `nthetas' values distributed in [0, 1]
+ * as [0, 1*S, 2*S, ..., i*S, ..., (nthetas-1)*S] with S=1/(nthetas-1). The
+ * ouput array `thetas' stored the angles corresponding to these cumulative
+ * values. Its length must be at least equal to `nthetas'. */
SSCHIFF_API res_T
sschiff_estimator_inverse_cumulative_phase_function
(const struct sschiff_estimator* estimator,
- const size_t iwlen,
- double* thetas,
+ const size_t wavelength_index,
+ double thetas[],
const size_t nthetas); /* Must be >= 2 */
END_DECLS
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -1226,11 +1226,19 @@ estimator_create
memset(Array, 0, sizeof(Array[0])*Count); \
} (void)0
- /* Setup and sort the wavelengths */
+ /* Check the wavelengths */
+ FOR_EACH(i, 1, nwavelengths) {
+ if(wavelengths[i] <= wavelengths[i-1]) {
+ log_error(dev,
+ "The submitted wavelengths are not sorted in ascending order.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ /* Copy the wavelengths */
RESIZE(estimator->wavelengths, nwavelengths,
"Couldn't allocate the list of estimated wavelengths.\n");
memcpy(estimator->wavelengths, wavelengths, nwavelengths*sizeof(double));
- qsort(estimator->wavelengths, nwavelengths, sizeof(double), cmp_double);
/* Generate the scattering angles */
RESIZE(estimator->angles, nangles,
@@ -1506,11 +1514,26 @@ sschiff_estimator_ref_put(struct sschiff_estimator* estimator)
}
res_T
-sschiff_estimator_get_wavelengths_count
- (const struct sschiff_estimator* estimator, size_t* count)
+sschiff_estimator_get_wavelengths
+ (const struct sschiff_estimator* estimator,
+ const double** wavelengths,
+ size_t* count)
{
- if(!estimator || !count) return RES_BAD_ARG;
- *count = sa_size(estimator->wavelengths);
+ if(!estimator) return RES_BAD_ARG;
+ if(wavelengths) *wavelengths = estimator->wavelengths;
+ if(count) *count = sa_size(estimator->wavelengths);
+ return RES_OK;
+}
+
+res_T
+sschiff_estimator_get_scattering_angles
+ (const struct sschiff_estimator* estimator,
+ const double** angles,
+ size_t* count)
+{
+ if(!estimator) return RES_BAD_ARG;
+ if(angles) *angles = estimator->angles;
+ if(count) *count = sa_size(estimator->angles);
return RES_OK;
}
@@ -1524,30 +1547,18 @@ sschiff_estimator_get_realisations_count
}
res_T
-sschiff_estimator_get_wavelength_cross_section
+sschiff_estimator_get_cross_section
(const struct sschiff_estimator* estimator,
- const double wavelength,
+ const size_t iwlen,
struct sschiff_cross_section* cross_section)
{
const struct mc_accum* acc;
- const double* wavelengths;
- const double* find;
- size_t iwavelength;
- size_t nwavelengths;
size_t N;
- if(!estimator || !cross_section)
+ if(!estimator || !cross_section || iwlen >= sa_size(estimator->wavelengths))
return RES_BAD_ARG;
- wavelengths = estimator->wavelengths;
- nwavelengths = sa_size(estimator->wavelengths);
- find = bsearch
- (&wavelength, wavelengths, nwavelengths, sizeof(double), cmp_double);
- if(!find) return RES_BAD_ARG;
-
- iwavelength = (size_t)(find - wavelengths);
- acc = estimator->accums + iwavelength;
-
+ acc = estimator->accums + iwlen;
N = estimator->nrealisations;
get_mc_value(&acc->extinction_cross_section, N, &cross_section->extinction);
get_mc_value(&acc->absorption_cross_section, N, &cross_section->absorption);
@@ -1581,6 +1592,31 @@ sschiff_estimator_get_cross_sections
}
res_T
+sschiff_estimator_get_phase_function
+ (const struct sschiff_estimator* estimator,
+ const size_t iwlen,
+ const struct sschiff_state** states)
+{
+ if(!estimator || !states || iwlen >= sa_size(estimator->wavelengths))
+ return RES_BAD_ARG;
+ *states = estimator->phase_functions[iwlen].values;
+ return RES_OK;
+}
+
+/* Retrieve a pointer onto the estimated phase function cumulative. */
+SSCHIFF_API res_T
+sschiff_estimator_get_phase_function_cumulative
+ (const struct sschiff_estimator* estimator,
+ const size_t iwlen,
+ const struct sschiff_state** states)
+{
+ if(!estimator || !states || iwlen >= sa_size(estimator->wavelengths))
+ return RES_BAD_ARG;
+ *states = estimator->phase_functions[iwlen].cumulative;
+ return RES_OK;
+}
+
+res_T
sschiff_estimator_inverse_cumulative_phase_function
(const struct sschiff_estimator* estimator,
const size_t iwlen,
@@ -1640,7 +1676,7 @@ sschiff_estimator_inverse_cumulative_phase_function
/* Inverse the phase function cumulative */
thetas[0] = 0.0;
thetas[nthetas-1] = PI;
- cumulative = step = 1.0/(double)nthetas;
+ cumulative = step = 1.0/(double)(nthetas-1);
FOR_EACH(itheta, 1, nthetas-1) {
/* TODO since the submitted cumulative are strictly increasing, one can
* speed up the inversion by using the previous search result to reduce the
diff --git a/src/test_sschiff_estimator_cylinder.c b/src/test_sschiff_estimator_cylinder.c
@@ -455,8 +455,8 @@ main(int argc, char** argv)
time_sub(&t0, &t1, &t0);
time_dump(&t0, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (estimator, wavelength, &cross_section), RES_OK);
+ CHECK(sschiff_estimator_get_cross_section
+ (estimator, 0, &cross_section), RES_OK);
printf("%u - x = %g - Wavelength = %g micron - %s\n",
(unsigned)i, x[i], wavelength, buf);
diff --git a/src/test_sschiff_estimator_rhodo.c b/src/test_sschiff_estimator_rhodo.c
@@ -47,7 +47,7 @@ get_material_property
const double wavelength,
struct sschiff_material_properties* props)
{
- (void)mtl;
+ (void)mtl, (void)wavelength;
#if 1
props->medium_refractive_index = 1.341779;
props->relative_imaginary_refractive_index = 2.20382600910679e-03/1.341779;
@@ -113,8 +113,8 @@ main(int argc, char** argv)
sschiff_uniform_scattering_angles, nscatt_angles, ngeoms, ndirs,
&estimator), RES_OK);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (estimator, wavelength, &cross_section), RES_OK);
+ CHECK(sschiff_estimator_get_cross_section
+ (estimator, 0, &cross_section), RES_OK);
printf("Wavelength = %g micron\n", wavelength);
printf(" Extinction ~ %9.3g +/- %9.3g\n",
diff --git a/src/test_sschiff_estimator_sphere.c b/src/test_sschiff_estimator_sphere.c
@@ -202,7 +202,9 @@ main(int argc, char** argv)
{ 21.0, 1350, 1580, 2930 }
};
const double wlen = 0.6; /* Micron */
- size_t count;
+ double wlens[3] = { 0, 0, 0 };
+ const double* dbls = NULL;
+ size_t i, count;
(void)argc, (void)argv;
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
@@ -294,36 +296,50 @@ main(int argc, char** argv)
CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1,
sschiff_uniform_scattering_angles, 1, 1, 1, &estimator), RES_BAD_ARG);
CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1,
- sschiff_uniform_scattering_angles, 3, 1, 1, &estimator), RES_OK);
+ sschiff_uniform_scattering_angles, 10, 1000, 1, &estimator), RES_OK);
- CHECK(sschiff_estimator_get_wavelengths_count(NULL, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelengths_count(estimator, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelengths_count(NULL, &count), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelengths_count(estimator, &count), RES_OK);
+ CHECK(sschiff_estimator_get_wavelengths(NULL, NULL, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelengths(estimator, NULL, NULL), RES_OK);/*Useless*/
+ CHECK(sschiff_estimator_get_wavelengths(NULL, NULL, &count), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelengths(estimator, NULL, &count), RES_OK);
+ CHECK(count, 1);
+ CHECK(sschiff_estimator_get_wavelengths(NULL, &dbls, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelengths(estimator, &dbls, NULL), RES_OK);
+ CHECK(dbls[0], wlen);
+ CHECK(sschiff_estimator_get_wavelengths(NULL, &dbls, &count), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_wavelengths(estimator, &dbls, &count), RES_OK);
+ CHECK(dbls[0], wlen);
CHECK(count, 1);
+ CHECK(sschiff_estimator_get_scattering_angles(NULL, NULL, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_scattering_angles(estimator, NULL, NULL), RES_OK); /* Useless */
+ CHECK(sschiff_estimator_get_scattering_angles(NULL, NULL, &count), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_scattering_angles(estimator, NULL, &count), RES_OK);
+ CHECK(count, 10);
+ CHECK(sschiff_estimator_get_scattering_angles(NULL, &dbls, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_scattering_angles(estimator, &dbls, NULL), RES_OK); /* Useless */
+ FOR_EACH(i, 0, count)
+ CHECK(eq_eps(dbls[i], (double)i*PI/((double)count-1), 1.e-6), 1);
+ CHECK(sschiff_estimator_get_scattering_angles(NULL, &dbls, &count), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_scattering_angles(estimator, &dbls, &count), RES_OK);
+ FOR_EACH(i, 0, count)
+ CHECK(eq_eps(dbls[i], (double)i*PI/((double)count-1), 1.e-6), 1);
+ CHECK(count, 10);
+
CHECK(sschiff_estimator_get_realisations_count(NULL, NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_get_realisations_count(estimator, NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_get_realisations_count(NULL, &count), RES_BAD_ARG);
CHECK(sschiff_estimator_get_realisations_count(estimator, &count), RES_OK);
- CHECK(count, 1);
+ CHECK(count, 1000);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (NULL, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (estimator, 0, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (NULL, wlen, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (estimator, wlen, NULL), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (NULL, 0, &cross_section), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (estimator, 0, &cross_section), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (NULL, wlen, &cross_section), RES_BAD_ARG);
- CHECK(sschiff_estimator_get_wavelength_cross_section
- (estimator, wlen, &cross_section), RES_OK);
+ CHECK(sschiff_estimator_get_cross_section(NULL, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_cross_section(estimator, 1, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_cross_section(NULL, 1, &cross_section), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_cross_section(estimator, 1, &cross_section), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_cross_section(NULL, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_cross_section(estimator, 0, NULL), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_cross_section(NULL, 0, &cross_section), RES_BAD_ARG);
+ CHECK(sschiff_estimator_get_cross_section(estimator, 0, &cross_section), RES_OK);
CHECK(sschiff_estimator_ref_get(NULL), RES_BAD_ARG);
CHECK(sschiff_estimator_ref_get(estimator), RES_OK);
@@ -331,11 +347,9 @@ main(int argc, char** argv)
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
- CHECK(sschiff_integrate(dev, rng, &distrib, &wlen, 1,
- sschiff_uniform_scattering_angles, 3, 10, 1, &estimator), RES_OK);
- CHECK(sschiff_estimator_get_realisations_count(estimator, &count), RES_OK);
- CHECK(count, 10);
- CHECK(sschiff_estimator_ref_put(estimator), RES_OK);
+ wlens[0] = 0.2, wlens[1] = 0.6, wlens[2] = 0.4;
+ CHECK(sschiff_integrate(dev, rng, &distrib, wlens, 3,
+ sschiff_uniform_scattering_angles, 10, 1000, 1, &estimator), RES_BAD_ARG);
sampler_ctx.relative_real_refractive_index = 1.1;
check_schiff_estimation(dev, rng, &sampler_ctx, results_n_r_1_1,
