rnatm

commit e72758ff719ce7368fc42b267c0aa1342f81f733
parent 902b3c458f884fb12adcd0cd28313539b28a295b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  3 Nov 2022 16:17:00 +0100

Add the rnatm_get_radcoef function

Diffstat:
M
src/rnatm.h
 | 
23
+++++++++++++++++++++++
M
src/rnatm_properties.c
 | 
83
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

2 files changed, 106 insertions(+), 0 deletions(-)
diff --git a/src/rnatm.h b/src/rnatm.h
@@ -103,6 +103,23 @@ struct rnatm_band_desc {
 static const struct rnatm_band_desc RNATM_BAND_DESC_NULL =
   RNATM_BAND_DESC_NULL__;
 
+struct rnatm_get_radcoef_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;
+
+  /* For debug: check that the retrieved radcoef is in [k_min, k_max] */
+  double k_min;
+  double k_max;
+};
+#define RNATM_GET_RADCOEF_ARGS_NULL__ {                                        \
+  {0,0,0}, 0, 0, RNATM_RADCOEFS_COUNT__, -DBL_MAX, DBL_MAX                     \
+}
+static const struct rnatm_get_radcoef_args
+RNATM_GET_RADCOEF_ARGS_NULL = RNATM_GET_RADCOEF_ARGS_NULL__;
+
 struct rnatm_cell_get_radcoef_args {
   struct rnatm_cell cell; /* Cell to query */
   double barycentric_coords[4]; /* Position into and relative to the cell */
@@ -248,6 +265,12 @@ rnatm_get_k_svx_voxel
    const enum rnatm_svx_op op);
 
 RNATM_API res_T
+rnatm_get_radcoef
+  (const struct rnatm* atm,
+   const struct rnatm_get_radcoef_args* args,
+   double* k);
+
+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
@@ -29,6 +29,7 @@
 #include <star/sbuf.h>
 #include <star/sck.h>
 #include <star/ssf.h>
+#include <star/suvm.h>
 
 #include <rsys/clock_time.h>
 #include <rsys/cstr.h>
@@ -37,6 +38,38 @@
  * Helper functions
  ******************************************************************************/
 static res_T
+check_rnatm_get_radcoef_args
+  (const struct rnatm* atm,
+   const struct rnatm_get_radcoef_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 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)
@@ -890,6 +923,56 @@ error:
  * Exported functions
  ******************************************************************************/
 res_T
+rnatm_get_radcoef
+  (const struct rnatm* atm,
+   const struct rnatm_get_radcoef_args* args,
+   double* out_k)
+{
+  struct rnatm_cell_get_radcoef_args cell_args = RNATM_CELL_GET_RADCOEF_ARGS_NULL;
+  size_t cpnt = RNATM_GAS;
+  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;
+
+  /* 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;
+
+  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);
+    if(res != RES_OK) goto error;
+
+    /* This component does not exist here */
+    if(SUVM_PRIMITIVE_NONE(&cell_args.cell.prim)) continue;
+
+    /* Add the component's contribution to the radiative coefficient */
+    res = rnatm_cell_get_radcoef(atm, &cell_args, &per_cell_k);
+    if(res != RES_OK) goto error;
+
+    k += per_cell_k;
+
+  } while(++cpnt < darray_aerosol_size_get(&atm->aerosols));
+
+  ASSERT(args->k_min <= k && k <= args->k_max);
+
+exit:
+  *out_k = k;
+  return res;
+error:
+  goto exit;
+}
+
+res_T
 rnatm_cell_get_radcoef
   (const struct rnatm* atm,
    const struct rnatm_cell_get_radcoef_args* args,