commit 055c2e753c48994b92ef0f847fbeb2e4635ebe42
parent a0d29fbba1e21a8d36a381d5e67c4e9def08c05e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 12 Jul 2018 16:57:30 +0200
Implement compute functions on cross sections
Compute the average or the boundaries for a given spectral band.
Diffstat:
| M | src/htmie.c | | | 193 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------- |
| M | src/htmie.h | | | 32 | ++++++++++++++++++++++++++++++++ |
2 files changed, 209 insertions(+), 16 deletions(-)
diff --git a/src/htmie.c b/src/htmie.c
@@ -449,6 +449,43 @@ massic_to_per_particle_xsections
}
static FINLINE double
+interpolate_data
+ (const struct htmie* htmie,
+ const double wavelength, /* Input wavelength */
+ const enum htmie_filter_type type, /* Interpolation type */
+ const size_t idata, /* Id of the upper bound data wrt `wavelength' */
+ const double* data) /* Data to interpolate */
+{
+ const double* wlens;
+ double a, b, u;
+ double val;
+ ASSERT(data);
+
+ wlens = htmie_get_wavelengths(htmie);
+
+ ASSERT(idata < htmie_get_wavelengths_count(htmie));
+ ASSERT(wlens[idata-1] < wavelength && wavelength <= wlens[idata-0]);
+
+ a = data[idata-1];
+ b = data[idata-0];
+ u = (wavelength - wlens[idata-1]) / (wlens[idata] - wlens[idata-1]);
+ ASSERT(eq_eps(u, 1, 1.e-6) || u < 1);
+ ASSERT(eq_eps(u, 0, 1.e-6) || u > 0);
+
+ switch(type) {
+ case HTMIE_FILTER_NEAREST:
+ val = u < 0.5 ? a : b;
+ break;
+ case HTMIE_FILTER_LINEAR:
+ CLAMP(u, 0, 1); /* Handle numerical inaccuracy */
+ val = u*b + (1-u)*a;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return val;
+}
+
+static FINLINE double
fetch_data
(const struct htmie* htmie,
const double wavelength,
@@ -474,26 +511,89 @@ fetch_data
val = data[0];
} else {
const size_t i = (size_t)(wlen_upp - wlens);
- const double a = data[i-1];
- const double b = data[i-0];
- const double u = (wavelength - wlens[i-1]) / (wlens[i] - wlens[i-1]);
- ASSERT(i < nwlens && (wlens[i] >= wavelength || wlens[i-1] <= wavelength));
- ASSERT(eq_eps(u, 1, 1.e-6) || u < 1);
-
- switch(type) {
- case HTMIE_FILTER_NEAREST:
- val = u < 0.5 ? a : b;
- break;
- case HTMIE_FILTER_LINEAR:
- CLAMP(u, 0, 1); /* Handle numerical inaccuracy */
- val = u*b + (1-u)*a;
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
+ val = interpolate_data(htmie, wavelength, type, i, data);
}
return val;
}
+static FINLINE void
+min_max(const double wavelength, const double data, void* ctx)
+{
+ double* bounds = ctx;
+ ASSERT(ctx);
+ (void)wavelength;
+ bounds[0] = MMIN(bounds[0], data);
+ bounds[1] = MMIN(bounds[1], data);
+}
+
+static FINLINE void
+add(const double wavelength, const double data, void* ctx)
+{
+ double* accum = ctx;
+ ASSERT(ctx);
+ (void)wavelength;
+ *accum += data;
+}
+
+static INLINE void
+for_each_wavelength_in_spectral_band
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of the spectral band in nanometer */
+ const enum htmie_filter_type type,
+ const double* data, /* Input data of the spectral band*/
+ void (*op)(const double wavelength, const double data, void* ctx),
+ void* context) /* Sent as the last argument of the 'op' functor */
+{
+ const double* wlens;
+ const double* wlen_low;
+ const double* wlen_upp;
+ double data_low;
+ double data_upp;
+ size_t ilow;
+ size_t iupp;
+ size_t nwlens;
+ size_t i;
+ ASSERT(htmie && band[0] <= band[1]);
+
+ wlens = htmie_get_wavelengths(htmie);
+ nwlens = htmie_get_wavelengths_count(htmie);
+ ASSERT(nwlens);
+
+ wlen_low = search_lower_bound
+ (&band[0], wlens, nwlens, sizeof(double), cmp_dbl);
+ wlen_upp = search_lower_bound
+ (&band[1], wlens, nwlens, sizeof(double), cmp_dbl);
+
+ if(!wlen_low) {
+ ilow = nwlens;
+ data_low = data[nwlens - 1];
+ } else if(wlen_low == wlens) {
+ ilow = 1;
+ data_low = data[0];
+ } else {
+ ilow = (size_t)(wlen_low - wlens);
+ data_low = interpolate_data(htmie, band[0], type, ilow, data);
+ }
+
+ if(!wlen_upp) {
+ iupp = nwlens;
+ data_upp = data[nwlens - 1];
+ } else if(wlen_upp == wlens) {
+ iupp = 1;
+ data_upp = data[0];
+ } else {
+ iupp = (size_t)(wlen_upp - wlens);
+ data_upp = interpolate_data(htmie, band[1], type, iupp, data);
+ }
+
+ op(band[0], data_low, context);
+ FOR_EACH(i, ilow, iupp) {
+ op(wlens[i], data[i], context);
+ }
+ if(band[0] != band[1])
+ op(band[1], data_upp, context);
+}
+
static void
release_htmie(ref_T* ref)
{
@@ -676,6 +776,67 @@ htmie_fetch_xsection_scattering
(htmie, wavelength, type, htmie_get_xsections_scattering(htmie));
}
+
+void
+htmie_compute_xsection_absorption_bounds
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of the spectral band in nanometer */
+ const enum htmie_filter_type type,
+ double bounds[2]) /* Min and Max scattering cross sections in m^2.kg^-1 */
+{
+ const double* xsections;
+ ASSERT(htmie);
+ bounds[0] = DBL_MAX;
+ bounds[1] =-DBL_MAX;
+ xsections = darray_double_cdata_get(&htmie->xsections_absorption);
+ for_each_wavelength_in_spectral_band
+ (htmie, band, type, xsections, min_max, bounds);
+}
+
+void
+htmie_compute_xsection_scattering_bounds
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of the spectral band in nanometer */
+ const enum htmie_filter_type type,
+ double bounds[2]) /* Min and Max scattering cross sections in m^2.kg^-1 */
+{
+ const double* xsections;
+ ASSERT(htmie);
+ bounds[0] = DBL_MAX;
+ bounds[1] =-DBL_MAX;
+ xsections = darray_double_cdata_get(&htmie->xsections_scattering);
+ for_each_wavelength_in_spectral_band
+ (htmie, band, type, xsections, min_max, bounds);
+}
+
+double
+htmie_compute_xsection_absorption_average
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of of the spectral band in nanometer */
+ const enum htmie_filter_type type)
+{
+ const double* xsections;
+ double sum = 0;
+ ASSERT(htmie && band[0] < band[1]);
+ xsections = darray_double_cdata_get(&htmie->xsections_absorption);
+ for_each_wavelength_in_spectral_band(htmie, band, type, xsections, add, &sum);
+ return sum / (band[1] - band[0]);
+}
+
+double
+htmie_compute_xsection_scattering_average
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of of the spectral band in nanometer */
+ const enum htmie_filter_type type)
+{
+ const double* xsections;
+ double sum = 0;
+ ASSERT(htmie && band[0] < band[1]);
+ xsections = darray_double_cdata_get(&htmie->xsections_scattering);
+ for_each_wavelength_in_spectral_band(htmie, band, type, xsections, add, &sum);
+ return sum / (band[1] - band[0]);
+}
+
double
htmie_fetch_asymmetric_parameter
(const struct htmie* htmie,
diff --git a/src/htmie.h b/src/htmie.h
@@ -111,6 +111,38 @@ htmie_fetch_xsection_scattering
const double wavelength,
const enum htmie_filter_type type);
+/* Compute the maximum and the minimum value of the absorption cross section in
+ * a given spectral interval */
+HTMIE_API void
+htmie_compute_xsection_absorption_bounds
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of the spectral band in nanometer */
+ const enum htmie_filter_type type,
+ double bounds[2]); /* Min and Max scattering cross sections in m^2.kg^-1 */
+
+/* Compute the maximum and the minimum value of the scattering cross section in
+ * a given spectral interval */
+HTMIE_API void
+htmie_compute_xsection_scattering_bounds
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of the spectral band in nanometer */
+ const enum htmie_filter_type type,
+ double bounds[2]); /* Min and Max scattering cross sections in m^2.kg^-1 */
+
+/* Compute the average absorption cross section in a given spectral interval */
+HTMIE_API double
+htmie_compute_xsection_absorption_average
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of of the spectral band in nanometer */
+ const enum htmie_filter_type type);
+
+/* Compute the average scattering cross section in a given spectral interval */
+HTMIE_API double
+htmie_compute_xsection_scattering_average
+ (const struct htmie* htmie,
+ const double band[2], /* Boundaries of of the spectral band in nanometer */
+ const enum htmie_filter_type type);
+
HTMIE_API double
htmie_fetch_asymmetric_parameter
(const struct htmie* htmie,
