commit bbff16dc0dafbf4f2b76e0b278d7b9948d835692
parent ae3f5f8178fb7c601a061671d03a1a7ada6f092f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 3 Mar 2016 16:16:41 +0100
Add several API functions
sschiff_estimator_get_limit_scattering_angle_index: returns the index of
the last small scattering angle
sschiff_estimator_get_wide_scattering_angles_model_parameter: returns
the estimated `n' parameter of the wide scattering angle model used to
analytically compute the differential cross section [cumulative].
sschiff_estimator_get_differential_cross_section: returns the
estimated differential cross section for a wavelength at a given
scattering angle.
sschiff_estimator_get_differential_cross_section_cumulative: returns the
estimated differential cross section cumulative for a wavelength at a
given scattering angle.
Diffstat:
2 files changed, 139 insertions(+), 6 deletions(-)
diff --git a/src/sschiff.h b/src/sschiff.h
@@ -232,6 +232,46 @@ sschiff_estimator_inverse_cumulative_phase_function
double thetas[],
const size_t nthetas); /* Must be >= 2 */
+
+/* Retrieve, for the submitted wavelength, the index of the last small
+ * scattering angle, i.e. the last angle from which the phase function
+ * [cumulative] is estimated by Monte Carlo. One can use this index to get the
+ * value of the limit angle from the scattering angles returned by the
+ * sschiff_estimator_get_scattering_angles function. Return RES_BAD_OP if the
+ * phase function is not computable for this wavelength. */
+SSCHIFF_API res_T
+sschiff_estimator_get_limit_scattering_angle_index
+ (const struct sschiff_estimator* estimator,
+ const size_t wavelength_index,
+ size_t* scattering_angle_index);
+
+/* Return the `n' parameter of the model used to compute the phase function of
+ * the wide scattering angles for a given wavelength. Return RES_BAD_OP if the
+ * phase function is not computable for this wavelength */
+SSCHIFF_API res_T
+sschiff_estimator_get_wide_scattering_angle_model_parameter
+ (const struct sschiff_estimator* estimator,
+ const size_t wavelength_index,
+ double* n);
+
+/* Return the estimated differential cross section for a wavelength at a given
+ * scattering angle. Return RES_BAD_OP if it was not computable. */
+SSCHIFF_API res_T
+sschiff_estimator_get_differential_cross_section
+ (const struct sschiff_estimator* estimator,
+ const size_t wavelength_index,
+ const size_t scattering_angle_index,
+ struct sschiff_state* state);
+
+/* Return the estimated differential cross section cumulative for a wavelength
+ * at a given scattering angle. Return RES_BAD_OP if it was not computable. */
+SSCHIFF_API res_T
+sschiff_estimator_get_differential_cross_section_cumulative
+ (const struct sschiff_estimator* estimator,
+ const size_t wavelength_index,
+ const size_t scanttering_angle_index,
+ struct sschiff_state* state);
+
END_DECLS
#endif /* SSCHIFF_H */
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -112,6 +112,7 @@ struct sschiff_estimator {
/* Per wavelength index of the scattering angle from which the MC estimation
* of the phase function is no more valid. Must be < nangles. */
size_t* limit_angles;
+ double* wide_angles_parameter;
double* wavenumbers;
struct sschiff_material_properties* properties; /* Material properties */
@@ -855,7 +856,7 @@ compute_phase_function
cos_limit_angle = cos(limit_angle);
cos_2_limit_angle = cos(2.0*limit_angle);
- /* Find the connector values between small and whide angles */
+ /* Find the connector values between small and wide angles */
get_mc_value
(&estimator->accums[iwlen].scattering_cross_section,
estimator->nrealisations,
@@ -883,7 +884,7 @@ compute_phase_function
log_error(estimator->dev,
"Couldn't estimate the parameters of the phase function for wide scattering\n"
"angles. The phase function is thus not computed.\n"
-"Wavelength = %g microns\n",
+"Wavelength = %g micron\n",
estimator->wavelengths[iwlen]);
goto error;
}
@@ -892,10 +893,12 @@ compute_phase_function
log_warning(estimator->dev,
"The wide scattering angles phase function model parameter `n' is not in\n"
"[2, 4]. One might increase the number of realisations.\n"
+" wavelength = %g micron\n"
" n = %g\n"
" scattering cross section = %g +/- %g\n"
" differential cross section = %g +/- %g\n"
" differential cross section cumulative = %g +/- %g\n",
+ estimator->wavelengths[iwlen],
n,
scattering_cross_section.E,
scattering_cross_section.SE,
@@ -905,6 +908,9 @@ compute_phase_function
limit_differential_cross_section_cumulative.SE);
}
+ /* Save the `n' wide angle parameter */
+ estimator->wide_angles_parameter[iwlen] = n;
+
/* Precomputed vlaue */
coef_limit = compute_differential_cross_section_cumulative_term
(sin_half_limit_angle, cos_limit_angle, cos_2_limit_angle, n);
@@ -1193,6 +1199,7 @@ estimator_release(ref_T* ref)
RELEASE(estimator->wavelengths);
RELEASE(estimator->angles);
RELEASE(estimator->limit_angles);
+ RELEASE(estimator->wide_angles_parameter);
RELEASE(estimator->wavenumbers);
RELEASE(estimator->properties);
if(estimator->accums) {
@@ -1290,6 +1297,9 @@ estimator_create
/* Allocate miscellaneous array of temporary/precomputed data */
RESIZE(estimator->limit_angles, nwavelengths,
"Couldn't allocate the per wavelength indices of the limit scattering angle.\n");
+ RESIZE(estimator->wide_angles_parameter, nwavelengths,
+ "Couldn't allocate the per wavelength `n' parameter of the wide scattering\n"
+ "angles model.\n");
RESIZE(estimator->wavenumbers, nwavelengths,
"Couldn't allocate the list of wavenumbers.\n");
RESIZE(estimator->properties, nwavelengths,
@@ -1343,9 +1353,12 @@ estimator_create
}
estimator->no_phase_function = !hold_valid_limit_angle;
if(estimator->no_phase_function) {
- /* No phase function => The scattering angles are useless */
+ /* No phase function => scattering angles and wide angle parameter are
+ * useless */
sa_release(estimator->angles);
+ sa_release(estimator->wide_angles_parameter);
estimator->angles = NULL;
+ estimator->wide_angles_parameter = NULL;
}
/* Allocate the estimator accumulators */
@@ -1377,7 +1390,6 @@ estimator_create
}
#undef RESIZE
-
exit:
*out_estimator = estimator;
return res;
@@ -1681,7 +1693,10 @@ sschiff_estimator_inverse_cumulative_phase_function
size_t iangle, nsmall_angles;
res_T res = RES_OK;
- if(!estimator || !nthetas < 2 || !thetas) {
+ if(!estimator
+ || nthetas < 2
+ || !thetas
+ || iwlen >= sa_size(estimator->wavelengths)) {
res = RES_BAD_ARG;
goto error;
}
@@ -1705,7 +1720,7 @@ sschiff_estimator_inverse_cumulative_phase_function
cumul[0] = &estimator->phase_functions[iwlen].cumulative[0];
cumulative_small_angles[0] = cumul[0]->E;
FOR_EACH(iangle, 1, nsmall_angles) {
- cumul[1] = &estimator->phase_functions[iwlen].cumulative[1];
+ cumul[1] = &estimator->phase_functions[iwlen].cumulative[iangle];
if(cumul[0]->E <= cumul[1]->E) {
cumulative_small_angles[iangle] = cumul[1]->E;
@@ -1737,6 +1752,7 @@ sschiff_estimator_inverse_cumulative_phase_function
* search domain in the cumulative arrays */
inverse_cumulative_phase_function
(estimator, cumulative_small_angles, iwlen, cumulative, &thetas[itheta]);
+ cumulative += step;
}
exit:
@@ -1745,3 +1761,80 @@ exit:
error:
goto exit;
}
+
+res_T
+sschiff_estimator_get_limit_scattering_angle_index
+ (const struct sschiff_estimator* estimator,
+ const size_t iwlen,
+ size_t* iangle)
+{
+ size_t limit_angle;
+ if(!estimator || iwlen >= sa_size(estimator->wavelengths) || !iangle)
+ return RES_BAD_ARG;
+ limit_angle = estimator->limit_angles[iwlen];
+ if(limit_angle == INVALID_LIMIT_ANGLE)
+ return RES_BAD_OP;
+ *iangle = limit_angle - 1/*<=>index of the last small angle*/;
+ return RES_OK;
+}
+
+res_T
+sschiff_estimator_get_wide_scattering_angle_model_parameter
+ (const struct sschiff_estimator* estimator,
+ const size_t iwlen,
+ double* out_n)
+{
+ if(!estimator || iwlen >= sa_size(estimator->wavelengths) || !out_n)
+ return RES_BAD_ARG;
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
+ return RES_BAD_OP;
+ *out_n = estimator->wide_angles_parameter[iwlen];
+ return RES_OK;
+}
+
+res_T
+sschiff_estimator_get_differential_cross_section
+ (const struct sschiff_estimator* estimator,
+ const size_t iwlen,
+ const size_t iangle,
+ struct sschiff_state* state)
+{
+ if(!estimator
+ || iwlen >= sa_size(estimator->wavelengths)
+ || iangle >= sa_size(estimator->angles)
+ || !state)
+ return RES_BAD_ARG;
+
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
+ return RES_BAD_OP;
+
+ get_mc_value
+ (&estimator->accums[iwlen].differential_cross_section[iangle],
+ estimator->nrealisations,
+ state);
+ return RES_OK;
+}
+
+res_T
+sschiff_estimator_get_differential_cross_section_cumulative
+ (const struct sschiff_estimator* estimator,
+ const size_t iwlen,
+ const size_t iangle,
+ struct sschiff_state* state)
+{
+ if(!estimator
+ || iwlen >= sa_size(estimator->wavelengths)
+ || iangle >= sa_size(estimator->angles)
+ || !state)
+ return RES_BAD_ARG;
+
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE)
+ return RES_BAD_OP;
+
+ get_mc_value
+ (&estimator->accums[iwlen].differential_cross_section_cumulative[iangle],
+ estimator->nrealisations,
+ state);
+ return RES_OK;
+}
+
