commit 909e743399ed2d7d5fc512e0496b5d8420644c2d
parent 6c0a54868890874ed5518a695fd206c81ed1b55d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 4 Nov 2022 12:31:25 +0100
Add the rnatm_sample_component function
Diffstat:
2 files changed, 164 insertions(+), 18 deletions(-)
diff --git a/src/rnatm.h b/src/rnatm.h
@@ -120,6 +120,20 @@ struct rnatm_get_radcoef_args {
static const struct rnatm_get_radcoef_args
RNATM_GET_RADCOEF_ARGS_NULL = RNATM_GET_RADCOEF_ARGS_NULL__;
+struct rnatm_sample_component_args {
+ double pos[3]; /* Position to be queried */
+ size_t iband; /* Index of the spectral band to consider */
+ size_t iquad; /* Index of the quadrature point to consider */
+
+ enum rnatm_radcoef radcoef;
+
+ double r; /* Random number uniformaly distributed in [0, 1[ */
+
+ /* For debug: check that the retrieved radcoef is in [k_min, k_max] */
+ double k_min;
+ double k_max;
+};
+
struct rnatm_cell_get_radcoef_args {
struct rnatm_cell cell; /* Cell to query */
double barycentric_coords[4]; /* Position into and relative to the cell */
@@ -271,6 +285,12 @@ rnatm_get_radcoef
double* k);
RNATM_API res_T
+rnatm_sample_component
+ (const struct rnatm* atm,
+ const struct rnatm_sample_component_args* args,
+ size_t* cpnt);
+
+RNATM_API res_T
rnatm_fetch_cell
(const struct rnatm* atm,
const double pos[3],
diff --git a/src/rnatm_properties.c b/src/rnatm_properties.c
@@ -34,6 +34,9 @@
#include <rsys/clock_time.h>
#include <rsys/cstr.h>
+#define CUMULATIVE_SIZE_MAX 32
+typedef double (cumulative_T)[CUMULATIVE_SIZE_MAX];
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -70,6 +73,42 @@ check_rnatm_get_radcoef_args
}
static res_T
+check_rnatm_sample_component_args
+ (const struct rnatm* atm,
+ const struct rnatm_sample_component_args* args)
+{
+ size_t i, n;
+
+ if(!args) return RES_BAD_ARG;
+
+ /* Invalid band index */
+ n = darray_band_size_get(&atm->bands) - 1;
+ if(args->iband < darray_band_cdata_get(&atm->bands)[0].index
+ || args->iband > darray_band_cdata_get(&atm->bands)[n].index)
+ return RES_BAD_ARG;
+
+ /* Invalid quadrature point index */
+ i = args->iband - darray_band_cdata_get(&atm->bands)[0].index;
+ ASSERT(i <= n && args->iband == darray_band_cdata_get(&atm->bands)[i].index);
+ if(args->iquad >= darray_band_cdata_get(&atm->bands)[i].nquad_pts)
+ return RES_BAD_ARG;
+
+ /* Invalid radiative coefficient */
+ if((unsigned)args->radcoef >= RNATM_RADCOEFS_COUNT__)
+ return RES_BAD_ARG;
+
+ /* Invalid random number */
+ if(args->r < 0 || args->r >= 1)
+ return RES_BAD_ARG;
+
+ /* Invalid K range */
+ if(args->k_min > args->k_max)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+static res_T
check_rnatm_cell_get_radcoef_args
(const struct rnatm* atm,
const struct rnatm_cell_get_radcoef_args* args)
@@ -923,37 +962,42 @@ error:
goto exit;
}
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-res_T
-rnatm_get_radcoef
+static res_T
+compute_unnormalized_cumulative_radcoef
(const struct rnatm* atm,
- const struct rnatm_get_radcoef_args* args,
- double* out_k)
+ const enum rnatm_radcoef radcoef,
+ const double pos[3],
+ const size_t iband,
+ const size_t iquad,
+ cumulative_T cumulative,
+ size_t* cumulative_sz)
{
struct rnatm_cell_get_radcoef_args cell_args = RNATM_CELL_GET_RADCOEF_ARGS_NULL;
size_t cpnt = RNATM_GAS;
+ size_t naerosols = 0;
+ size_t icumul = 0;
double k = 0;
res_T res = RES_OK;
+ ASSERT(atm && (unsigned)radcoef < RNATM_RADCOEFS_COUNT__);
+ ASSERT(cumulative && cumulative_sz);
- if(!atm || !out_k) { res = RES_BAD_ARG; goto error; }
- res = check_rnatm_get_radcoef_args(atm, args);
- if(res != RES_OK) goto error;
+ naerosols = darray_aerosol_size_get(&atm->aerosols);
+ if(naerosols+1/*gas*/ > CUMULATIVE_SIZE_MAX) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
- /* Setup the arguments common to all component */
- cell_args.iband = args->iband;
- cell_args.iquad = args->iquad;
- cell_args.radcoef = args->radcoef;
- cell_args.k_min = args->k_min;
- cell_args.k_max = args->k_max;
+ /* Setup the arguments common to all components */
+ cell_args.iband = iband;
+ cell_args.iquad = iquad;
+ cell_args.radcoef = radcoef;
do {
double per_cell_k;
/* Retrieve the cell of the component */
res = rnatm_fetch_cell
- (atm, args->pos, cpnt, &cell_args.cell, cell_args.barycentric_coords);
+ (atm, pos, cpnt, &cell_args.cell, cell_args.barycentric_coords);
if(res != RES_OK) goto error;
/* This component does not exist here */
@@ -965,8 +1009,45 @@ rnatm_get_radcoef
k += per_cell_k;
- } while(++cpnt < darray_aerosol_size_get(&atm->aerosols));
+ cumulative[icumul] = k;
+ ++icumul;
+
+ } while(++cpnt < naerosols);
+ *cumulative_sz = naerosols + 1;
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+rnatm_get_radcoef
+ (const struct rnatm* atm,
+ const struct rnatm_get_radcoef_args* args,
+ double* out_k)
+{
+ cumulative_T cumul;
+ size_t cumul_sz;
+ double k = 0;
+ res_T res = RES_OK;
+
+ if(!atm || !out_k) { res = RES_BAD_ARG; goto error; }
+ res = check_rnatm_get_radcoef_args(atm, args);
+ if(res != RES_OK) goto error;
+
+ /* Calculate the cumulative (unormalized) of radiative coefficients. Its last
+ * entry is the sum of the radiative coefficients of each component, which is
+ * the atmospheric radiative coefficient to be returned. */
+ res = compute_unnormalized_cumulative_radcoef
+ (atm, args->radcoef, args->pos, args->iband, args->iquad, cumul, &cumul_sz);
+ if(res != RES_OK) goto error;
+ ASSERT(cumul_sz >= 1);
+
+ k = cumul[cumul_sz-1];
ASSERT(args->k_min <= k && k <= args->k_max);
exit:
@@ -977,6 +1058,51 @@ error:
}
res_T
+rnatm_sample_component
+ (const struct rnatm* atm,
+ const struct rnatm_sample_component_args* args,
+ size_t* cpnt)
+{
+ cumulative_T cumul;
+ double rcp_norm;
+ size_t cumul_sz;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!atm || !cpnt) { res = RES_BAD_ARG; goto error; }
+ res = check_rnatm_sample_component_args(atm, args);
+ if(res != RES_OK) goto error;
+
+ /* Discard the calculation of the cumulative if there is only gas */
+ if(!darray_aerosol_size_get(&atm->aerosols)) {
+ *cpnt = RNATM_GAS;
+ goto exit;
+ }
+
+ res = compute_unnormalized_cumulative_radcoef
+ (atm, args->radcoef, args->pos, args->iband, args->iquad, cumul, &cumul_sz);
+ if(res != RES_OK) goto error;
+ ASSERT(cumul_sz >= 1);
+
+ /* Normalize the cumulative */
+ rcp_norm = 1.0/cumul[cumul_sz-1];
+ FOR_EACH(i, 0, cumul_sz) cumul[i] *= rcp_norm;
+ cumul[cumul_sz-1] = 1.0; /* Handle precision issues */
+
+ /* Use a simple linear search to sample the component since there are
+ * usually very few aerosols to consider */
+ FOR_EACH(i, 0, cumul_sz) if(args->r < cumul[i]) break;
+
+ *cpnt = i-1;
+ ASSERT(*cpnt == RNATM_GAS || *cpnt < darray_aerosol_size_get(&atm->aerosols));
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
rnatm_cell_get_radcoef
(const struct rnatm* atm,
const struct rnatm_cell_get_radcoef_args* args,
