star-sf

Set of surface and volume scattering functions
git clone git://git.meso-star.fr/star-sf.git
Log | Files | Refs | README | LICENSE

commit 3e2fe000de4c03b4a7c2be7eaa0ada2270ddb689
parent 1f42db3f055ebd8f8db174dd0f7ba2a6e8dc6360
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  5 Feb 2021 13:18:59 +0100

Add the ssf_phase_rdgfa_get_<desc|interval> functions

Diffstat:
Msrc/ssf.h | 40++++++++++++++++++++++++++++++++++++++--
Msrc/ssf_phase_rdgfa.c | 102+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------
2 files changed, 125 insertions(+), 17 deletions(-)

diff --git a/src/ssf.h b/src/ssf.h @@ -17,6 +17,7 @@ #define SSF_H #include <rsys/rsys.h> +#include <float.h> /* Library symbol management */ #if defined(SSF_SHARED_BUILD) /* Build shared library */ @@ -194,9 +195,30 @@ struct ssf_phase_rdgfa_setup_args { /* #intervals used to discretized the both angular domains */ size_t nintervals[2]; }; -#define SSF_PHASE_RDGFA_SETUP_DEFAULT_ARGS__ {0, 0, 0, {1000, 1000}} +#define SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT__ {0,0,0,{1000,1000}} static const struct ssf_phase_rdgfa_setup_args -SSF_PHASE_RDGFA_SETUP_DEFAULT_ARGS = SSF_PHASE_RDGFA_SETUP_DEFAULT_ARGS__; +SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT = SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT__; + +struct ssf_phase_rdgfa_desc { + double wavelength; /* In nm */ + double fractal_dimension; /* No unit */ + double gyration_radius; /* In nm */ + + double theta_limit; /* Angle that splits the 2 angular domains. In rad */ + double normalization_factor; /* Normalization factor of the cumulative */ + size_t nintervals; /* #intervals used to discretized the cumulative */ +}; +#define SSF_PHASE_RDGFA_DESC_NULL__ {0,0,0,0,0,0} +static const struct ssf_phase_rdgfa_desc SSF_PHASE_RDGFA_DESC_NULL = + SSF_PHASE_RDGFA_DESC_NULL__; + +struct ssf_phase_rdgfa_interval { + double range[2]; /* Angular range of the interval. In rad */ + double cumulative; /* Value of the cumulative of the interval */ +}; +#define SSF_PHASE_RDGFA_INTERVAL_NULL__ {{DBL_MAX,-DBL_MAX}, 0} +static const struct ssf_phase_rdgfa_interval SSF_PHASE_RDGFA_INTERVAL_NULL = + SSF_PHASE_RDGFA_INTERVAL_NULL__; BEGIN_DECLS @@ -464,11 +486,25 @@ ssf_phase_hg_setup (struct ssf_phase* phase, const double g); /* Asymmetric parameter in [-1, 1] */ +/******************************************************************************* + * RDGFA phase function + ******************************************************************************/ SSF_API res_T ssf_phase_rdgfa_setup (struct ssf_phase* phase, const struct ssf_phase_rdgfa_setup_args* args); +SSF_API res_T +ssf_phase_rdgfa_get_desc + (const struct ssf_phase* phase, + struct ssf_phase_rdgfa_desc* desc); + +SSF_API res_T +ssf_phase_rdgfa_get_interval + (const struct ssf_phase* phase, + const size_t interval_id, /* In [0, desc.nintervals[ */ + struct ssf_phase_rdgfa_interval* interval); + /******************************************************************************* * Fresnel API - Define the equation of the fresnel term ******************************************************************************/ diff --git a/src/ssf_phase_rdgfa.c b/src/ssf_phase_rdgfa.c @@ -78,6 +78,29 @@ cmp_dbl(const void* a, const void* b) } } +/* Compute the angular range of an interval */ +static INLINE void +compute_interval_angular_range + (const struct rdgfa* rdgfa, + const size_t interval, + double angular_range[2]) +{ + ASSERT(rdgfa && interval < rdgfa->nintervals[0] + rdgfa->nintervals[1]); + ASSERT(angular_range); + + /* First angular domain in [0, theta_limit[ */ + if(interval < rdgfa->nintervals[0]) { + angular_range[0] = rdgfa->dtheta[0] * (double)(interval + 0); + angular_range[1] = rdgfa->dtheta[0] * (double)(interval + 1); + + /* Second angular domain in [theta_limit, PI[ */ + } else { + const size_t i = interval - rdgfa->nintervals[0]; + angular_range[0] = rdgfa->dtheta[1] * (double)(i + 0) + rdgfa->theta_limit; + angular_range[1] = rdgfa->dtheta[1] * (double)(i + 1) + rdgfa->theta_limit; + } +} + #if 0 /* Not normalized */ static INLINE double @@ -267,8 +290,7 @@ rdgfa_sample double* find = NULL; double frame[9]; double w[3]; - double theta1; - double theta2; + double thetas[2]; double theta; double sin_theta; double phi; @@ -288,16 +310,8 @@ rdgfa_sample ASSERT(find && (*find) >= r); i = (size_t)(find - cdf); - /* Compute the angle interval into which the sample lies */ - if(i < rdgfa->nintervals[0]) { - theta1 = rdgfa->dtheta[0] * (double)(i + 0); - theta2 = rdgfa->dtheta[0] * (double)(i + 1); - } else { - theta1 = rdgfa->dtheta[1] * (double)(i - rdgfa->nintervals[0] + 0); - theta2 = rdgfa->dtheta[1] * (double)(i - rdgfa->nintervals[0] + 1); - theta1 = theta1 + rdgfa->theta_limit; - theta2 = theta2 + rdgfa->theta_limit; - } + /* Compute the angular range into which the sample lies */ + compute_interval_angular_range(rdgfa, i, thetas); /* Compute the sampled theta angle by linearly interpolate it from the * boundaries of its interval */ @@ -307,7 +321,7 @@ rdgfa_sample u = (r - cdf[i-1]) / (cdf[i] - cdf[i-1]); } ASSERT(0 <= u && u < 1); - theta = u * (theta2 - theta1) + theta1; + theta = u * (thetas[1] - thetas[0]) + thetas[0]; /* Uniformly sample a phi angle in [0, 2PI[ */ phi = ssp_rng_uniform_double(rng, 0, 2*PI); @@ -351,7 +365,7 @@ ssf_phase_rdgfa_setup double lambda, Df, Rg, k, k2; res_T res = RES_OK; - if(!phase + if(!phase || !PHASE_TYPE_EQ(&phase->type, &ssf_phase_rdgfa) || !check_phase_rdgfa_setup_args(args)) { res = RES_BAD_ARG; @@ -363,13 +377,14 @@ ssf_phase_rdgfa_setup rdgfa->fractal_dimension = args->fractal_dimension; rdgfa->gyration_radius = args->gyration_radius; rdgfa->nintervals[0] = args->nintervals[0]; + rdgfa->nintervals[1] = args->nintervals[1]; /* Fetch input data */ lambda = rdgfa->wavelength; Df = rdgfa->fractal_dimension; Rg = rdgfa->gyration_radius; - /* Precompute constants */ + /* Precompute constants */ rdgfa->Rg2 = Rg*Rg; /* gyration_radius^2 */ rdgfa->cst_4pi_div_lambda = 4*PI / lambda; rdgfa->cst_3Df_div_2E = 3*Df/(2.0*EXP1); @@ -388,3 +403,60 @@ exit: error: goto exit; } + +res_T +ssf_phase_rdgfa_get_desc + (const struct ssf_phase* phase, + struct ssf_phase_rdgfa_desc* desc) +{ + struct rdgfa* rdgfa = NULL; + res_T res = RES_OK; + + if(!phase || !PHASE_TYPE_EQ(&phase->type, &ssf_phase_rdgfa) || !desc) { + res = RES_BAD_ARG; + goto error; + } + + rdgfa = phase->data; + desc->wavelength = rdgfa->wavelength; + desc->gyration_radius = rdgfa->gyration_radius; + desc->fractal_dimension = rdgfa->fractal_dimension; + desc->theta_limit = rdgfa->theta_limit; + desc->normalization_factor = 1.0/rdgfa->rcp_normalize_factor; + desc->nintervals = rdgfa->nintervals[0] + rdgfa->nintervals[1]; + +exit: + return res; +error: + goto exit; +} + +res_T +ssf_phase_rdgfa_get_interval + (const struct ssf_phase* phase, + const size_t interval_id, /* In [0, #intervals[ */ + struct ssf_phase_rdgfa_interval* interval) +{ + struct rdgfa* rdgfa = NULL; + res_T res = RES_OK; + + if(!phase || !PHASE_TYPE_EQ(&phase->type, &ssf_phase_rdgfa) || !interval) { + res = RES_BAD_ARG; + goto error; + } + + rdgfa = phase->data; + if(interval_id >= darray_double_size_get(&rdgfa->cdf)) { + res = RES_BAD_ARG; + goto error; + } + + compute_interval_angular_range(rdgfa, interval_id, interval->range); + interval->cumulative = darray_double_cdata_get(&rdgfa->cdf)[interval_id]; + +exit: + return res; +error: + goto exit; +} +