commit ca72a7c2c598c9abbb2221d7be00e036529a575a
parent 1903657bfc8786ea103eb99303b9c6e690bc520d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 26 Sep 2022 16:43:55 +0200
Add the rnatm_cell_get_radcoef function
Internal refactoring to use the same code to get the radiative
coefficients whether it was when building octree or when calling the
public API.
Diffstat:
6 files changed, 393 insertions(+), 143 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -70,6 +70,7 @@ set(RNATM_FILES_SRC
rnatm_octree.c
rnatm_octrees_storage.c
rnatm_properties.c
+ rnatm_radcoef.c
rnatm_voxel_partition.c
rnatm_write_vtk.c)
set(RNATM_FILES_INC
diff --git a/src/rnatm.h b/src/rnatm.h
@@ -91,6 +91,31 @@ struct rnatm_cell {
* the RNATM_GAS constant */
size_t component;
};
+#define RNATM_CELL_NULL__ {SUVM_PRIMITIVE_NULL__, RNATM_GAS}
+static const struct rnatm_cell RNATM_CELL_NULL = RNATM_CELL_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 */
+ size_t iband; /* Index of the spectral band to query */
+ size_t iquad; /* Index of the quadrature point to query in the band */
+ double wavelength; /* Wavenlength in the queried band (in nanometers) */
+
+ 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_CELL_GET_RADCOEF_ARGS_NULL__ { \
+ RNATM_CELL_NULL__, \
+ {0, 0, 0, 0}, /* Barycentric coordinates */ \
+ 0, 0, 0, /* Spectral data (band, quadrature point, wavelength) */ \
+ RNATM_RADCOEFS_COUNT__, \
+ -DBL_MAX, DBL_MAX /* For debug: Radcoef range */ \
+}
+static const struct rnatm_cell_get_radcoef_args
+RNATM_CELL_GET_RADCOEF_ARGS_NULL = RNATM_CELL_GET_RADCOEF_ARGS_NULL__;
struct rnatm_create_args {
struct rnatm_gas_args gas;
@@ -189,6 +214,12 @@ rnatm_fetch_cell
struct rnatm_cell* cell,
double barycentric_coords[4]); /* `pos' relative to the returned `cell' */
+RNATM_API res_T
+rnatm_cell_get_radcoef
+ (const struct rnatm* atm,
+ const struct rnatm_cell_get_radcoef_args* args,
+ double* k);
+
RNATM_API size_t
rnatm_get_aerosols_count
(const struct rnatm* atm);
diff --git a/src/rnatm_c.h b/src/rnatm_c.h
@@ -21,12 +21,13 @@
#ifndef RNATM_C_H
#define RNATM_C_H
+#include "rnatm.h"
+
#include <rsys/dynamic_array_size_t.h>
#include <rsys/logger.h>
#include <rsys/ref_count.h>
#include <rsys/str.h>
-struct rnatm_create_args;
struct rnsf;
/* Generate the dynamic array of dynamic array of size_t */
@@ -251,4 +252,18 @@ unload_octree
(struct rnatm* atm,
const size_t iaccel_struct);
+/* Get the radiative coefficient of the input tetrahedron. The output radiative
+ * coefficients per vertex are defined even if the tetrahedron is invalid or if
+ * no spectral data is defined for this component for the considered band. In
+ * such cases, the k returned are zero and the function returns 0 instead of 1*/
+extern LOCAL_SYM int
+tetra_get_radcoef
+ (const struct rnatm* atm,
+ const struct suvm_primitive* tetra,
+ const size_t component,
+ const size_t iband,
+ const size_t iquad_pt,
+ const enum rnatm_radcoef radcoef,
+ float k[4]);
+
#endif /* RNATM_C_H */
diff --git a/src/rnatm_octree.c b/src/rnatm_octree.c
@@ -18,8 +18,7 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
-/* lround, nextafter */
-#define _POSIX_C_SOURCE 200112L
+#define _POSIX_C_SOURCE 200112L /* lround */
#include "rnatm_c.h"
#include "rnatm_log.h"
@@ -386,160 +385,35 @@ compute_grid_definition(struct rnatm* atm, const struct rnatm_create_args* args)
return RES_OK;
}
+/* Compute the radiative coefficients range of the tetrahedron */
static INLINE void
-setup_tetra_radcoefs_gas
+setup_tetra_radcoefs
(struct rnatm* atm,
const struct suvm_primitive* tetra,
+ const size_t cpnt,
const size_t iband,
- const size_t iquad_pt,
+ const size_t iquad,
struct radcoefs* radcoefs)
{
float ka[4];
float ks[4];
- struct sck_band band;
- struct sck_quad_pt quad_pt;
ASSERT(atm && tetra && radcoefs);
- ASSERT(tetra->nvertices == 4);
-
- /* Compute the scattering coefficient range of the tetrahedron */
- SCK(get_band(atm->gas.ck, iband, &band));
- ks[0] = band.ks_list[tetra->indices[0]];
- ks[1] = band.ks_list[tetra->indices[1]];
- ks[2] = band.ks_list[tetra->indices[2]];
- ks[3] = band.ks_list[tetra->indices[3]];
- radcoefs->ks_min = MMIN(MMIN(ks[0], ks[1]), MMIN(ks[2], ks[3]));
- radcoefs->ks_max = MMAX(MMAX(ks[0], ks[1]), MMAX(ks[2], ks[3]));
-
- /* Compute the absorption coefficient range of the tetrahedron */
- SCK(band_get_quad_pt(&band, iquad_pt, &quad_pt));
- ka[0] = quad_pt.ka_list[tetra->indices[0]];
- ka[1] = quad_pt.ka_list[tetra->indices[1]];
- ka[2] = quad_pt.ka_list[tetra->indices[2]];
- ka[3] = quad_pt.ka_list[tetra->indices[3]];
- radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
- radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
-}
-
-static void
-setup_tetra_radcoefs_aerosol
- (struct rnatm* atm,
- const struct suvm_primitive* tetra,
- const size_t iband,
- const struct aerosol* aerosol,
- struct radcoefs* radcoefs)
-{
- struct sck_band gas_band;
- struct sars_band ars_band;
- double gas_spectral_range[2]; /* In nm */
- size_t ars_ibands[2];
- float ka[4], ks[4];
- ASSERT(atm && aerosol && tetra && radcoefs);
-
- /* Look for the aerosol bands covered by the gas band */
- SCK(get_band(atm->gas.ck, iband, &gas_band));
- gas_spectral_range[0] = gas_band.lower;
- gas_spectral_range[1] = nextafter(gas_band.upper, 0); /* inclusive */
- SARS(find_bands(aerosol->sars, gas_spectral_range, ars_ibands));
-
- /* No aerosol band is overlaid by the gas band: ignore the aerosol */
- if(ars_ibands[0] > ars_ibands[1]) {
- radcoefs->ks_min = FLT_MAX;
- radcoefs->ks_max =-FLT_MAX;
- radcoefs->ka_min = FLT_MAX;
- radcoefs->ka_max =-FLT_MAX;
- return;
- }
-
- SARS(get_band(aerosol->sars, ars_ibands[0], &ars_band));
-
- /* The gas band is included in an aerosol band: use the aerosol radiative
- * coefficients */
- if(ars_ibands[0] == ars_ibands[1]
- && ars_band.lower <= gas_band.lower
- && ars_band.upper >= gas_band.upper) {
+ ASSERT(cpnt == RNATM_GAS || cpnt < darray_aerosol_size_get(&atm->aerosols));
- /* Compute the absorption coefficient range of the tetrahedron */
- ka[0] = sars_band_get_ka(&ars_band, tetra->indices[0]);
- ka[1] = sars_band_get_ka(&ars_band, tetra->indices[1]);
- ka[2] = sars_band_get_ka(&ars_band, tetra->indices[2]);
- ka[3] = sars_band_get_ka(&ars_band, tetra->indices[3]);
+ if(tetra_get_radcoef(atm, tetra, cpnt, iband, iquad, RNATM_RADCOEF_Ka, ka)
+ && tetra_get_radcoef(atm, tetra, cpnt, iband, iquad, RNATM_RADCOEF_Ks, ks)) {
+ /* Radiative coefficients are defined for this component at the considered
+ * spectral band */
radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
-
- /* Compute the scattering coefficient range of the tetrahedron */
- ks[0] = sars_band_get_ks(&ars_band, tetra->indices[0]);
- ks[1] = sars_band_get_ks(&ars_band, tetra->indices[1]);
- ks[2] = sars_band_get_ks(&ars_band, tetra->indices[2]);
- ks[3] = sars_band_get_ks(&ars_band, tetra->indices[3]);
radcoefs->ks_min = MMIN(MMIN(ks[0], ks[1]), MMIN(ks[2], ks[3]));
radcoefs->ks_max = MMAX(MMAX(ks[0], ks[1]), MMAX(ks[2], ks[3]));
-
- /* The gas band overlaid N aerosol bands (N >= 1) */
} else {
- float tau_ka[4] = {0, 0, 0, 0};
- float tau_ks[4] = {0, 0, 0, 0};
- double lambda_min;
- double lambda_max;
- float rcp_gas_band_len;
- size_t iars_band;
-
- FOR_EACH(iars_band, ars_ibands[0], ars_ibands[1]+1) {
- float lambda_len;
-
- SARS(get_band(aerosol->sars, iars_band, &ars_band));
- lambda_min = MMAX(gas_band.lower, ars_band.lower);
- lambda_max = MMIN(gas_band.upper, ars_band.upper); /* exclusive */
- lambda_max = nextafter(lambda_max, 0); /* inclusive */
-
- /* Shrink lambda_<min|max> to the spectral range */
- lambda_min = MMAX(lambda_min, atm->spectral_range[0]);
- lambda_max = MMIN(lambda_max, atm->spectral_range[1]);
- lambda_len = (float)(lambda_max - lambda_min);
- ASSERT(lambda_len > 0);
-
- tau_ka[0] += sars_band_get_ka(&ars_band, tetra->indices[0]) * lambda_len;
- tau_ka[1] += sars_band_get_ka(&ars_band, tetra->indices[1]) * lambda_len;
- tau_ka[2] += sars_band_get_ka(&ars_band, tetra->indices[2]) * lambda_len;
- tau_ka[3] += sars_band_get_ka(&ars_band, tetra->indices[3]) * lambda_len;
-
- tau_ks[0] += sars_band_get_ks(&ars_band, tetra->indices[0]) * lambda_len;
- tau_ks[1] += sars_band_get_ks(&ars_band, tetra->indices[1]) * lambda_len;
- tau_ks[2] += sars_band_get_ks(&ars_band, tetra->indices[2]) * lambda_len;
- tau_ks[3] += sars_band_get_ks(&ars_band, tetra->indices[3]) * lambda_len;
- }
-
- /* Compute the radiative coefficients of the tetrahedron */
- lambda_min = gas_band.lower;
- lambda_max = nextafter(gas_band.upper, 0); /* inclusive */
- rcp_gas_band_len = 1.f/(float)(lambda_max - lambda_min);
- f4_mulf(ks, tau_ks, rcp_gas_band_len);
- f4_mulf(ka, tau_ka, rcp_gas_band_len);
-
- /* Compute the radiative coefficients range of the tetrahedron */
- radcoefs->ks_min = MMIN(MMIN(ks[0], ks[1]), MMIN(ks[2], ks[3]));
- radcoefs->ks_max = MMAX(MMAX(ks[0], ks[1]), MMAX(ks[2], ks[3]));
- radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
- radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
- }
-}
-
-static INLINE void
-setup_tetra_radcoefs
- (struct rnatm* atm,
- const struct suvm_primitive* tetra,
- const size_t iband,
- const size_t iquad_pt,
- const size_t cpnt,
- struct radcoefs* radcoefs)
-{
- ASSERT(atm && tetra && radcoefs);
- ASSERT(cpnt == RNATM_GAS || cpnt < darray_aerosol_size_get(&atm->aerosols));
-
- if(cpnt == RNATM_GAS) {
- setup_tetra_radcoefs_gas(atm, tetra, iband, iquad_pt, radcoefs);
- } else {
- const struct aerosol* aerosol = darray_aerosol_cdata_get(&atm->aerosols)+cpnt;
- setup_tetra_radcoefs_aerosol(atm, tetra, iband, aerosol, radcoefs);
+ /* No valid radiative coefficients are set: ignore the tetrahedron */
+ radcoefs->ka_min = FLT_MAX;
+ radcoefs->ka_max =-FLT_MAX;
+ radcoefs->ks_min = FLT_MAX;
+ radcoefs->ks_max =-FLT_MAX;
}
}
@@ -642,7 +516,7 @@ voxelize_tetra
const size_t iband = args->accel_structs[iitem].iband;
const size_t iquad_pt = args->accel_structs[iitem].iquad_pt;
struct radcoefs* radcoefs = args->per_item_radcoefs + iitem;
- setup_tetra_radcoefs(atm, &tetra, iband, iquad_pt, cpnt, radcoefs);
+ setup_tetra_radcoefs(atm, &tetra, cpnt, iband, iquad_pt, radcoefs);
}
/* Iterate voxels intersected by the AABB of the polyedron */
diff --git a/src/rnatm_properties.c b/src/rnatm_properties.c
@@ -35,6 +35,61 @@
/*******************************************************************************
* Helper functions
******************************************************************************/
+static res_T
+check_rnatm_cell_get_radcoef_args
+ (const struct rnatm* atm,
+ const struct rnatm_cell_get_radcoef_args* args)
+{
+ struct sck_band band;
+ double sum_bcoords;
+ size_t i, n;
+ ASSERT(darray_band_size_get(&atm->bands));
+
+ if(!args) return RES_BAD_ARG;
+
+ /* Invalid geometric primitive */
+ if(SUVM_PRIMITIVE_NONE(&args->cell.prim)
+ || args->cell.prim.nvertices != 4) /* Expect a tetraheron */
+ return RES_BAD_ARG;
+
+ /* Invalid component */
+ if(args->cell.component != RNATM_GAS
+ && args->cell.component >= darray_aerosol_size_get(&atm->aerosols))
+ return RES_BAD_ARG;
+
+ /* Invalid barycentric coordinates */
+ sum_bcoords =
+ args->barycentric_coords[0]
+ + args->barycentric_coords[1]
+ + args->barycentric_coords[2]
+ + args->barycentric_coords[3];
+ if(!eq_eps(sum_bcoords, 1, 1.e-6))
+ 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 wavelength */
+ SCK(get_band(atm->gas.ck, args->iband, &band));
+ if(args->wavelength < band.lower || args->wavelength > band.upper)
+ return RES_BAD_ARG;
+
+ /* Invalid K range */
+ if(args->k_min > args->k_max)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
static INLINE void
reset_gas_properties(struct gas* gas)
{
@@ -498,6 +553,42 @@ error:
}
/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+rnatm_cell_get_radcoef
+ (const struct rnatm* atm,
+ const struct rnatm_cell_get_radcoef_args* args,
+ double* out_k)
+{
+ float per_vertex_k[4];
+ double k = NaN;
+ res_T res = RES_OK;
+
+ if(!atm || !out_k) { res = RES_BAD_ARG; goto error; }
+ res = check_rnatm_cell_get_radcoef_args(atm, args);
+ if(res != RES_OK) goto error;
+
+ /* Get the radiative coefficients of the tetrahedron */
+ tetra_get_radcoef(atm, &args->cell.prim, args->cell.component, args->iband,
+ args->iquad, args->radcoef, per_vertex_k);
+
+ /* Calculate the radiative coefficient by linearly interpolating the
+ * coefficients defined by tetrahedron vertex */
+ k = per_vertex_k[0] * args->barycentric_coords[0]
+ + per_vertex_k[1] * args->barycentric_coords[1]
+ + per_vertex_k[2] * args->barycentric_coords[2]
+ + per_vertex_k[3] * args->barycentric_coords[3];
+ ASSERT(args->k_min <= k && k <= args->k_max);
+
+exit:
+ *out_k = k;
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
* Local functions
******************************************************************************/
res_T
diff --git a/src/rnatm_radcoef.c b/src/rnatm_radcoef.c
@@ -0,0 +1,238 @@
+/* Copyright (C) 2022 Centre National de la Recherche Scientifique
+ * Copyright (C) 2022 Institut de Physique du Globe de Paris
+ * Copyright (C) 2022 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2022 Université de Reims Champagne-Ardenne
+ * Copyright (C) 2022 Université de Versaille Saint-Quentin
+ * Copyright (C) 2022 Université Paul Sabatier (contact@laplace.univ-tlse.fr)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "rnatm_c.h"
+
+#include <star/sars.h>
+#include <star/sck.h>
+
+#include <rsys/float4.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static int
+tetra_get_radcoef_gas
+ (const struct rnatm* atm,
+ const struct suvm_primitive* tetra,
+ const size_t iband,
+ const size_t iquad_pt,
+ const enum rnatm_radcoef radcoef,
+ float k[4])
+{
+ struct sck_band band;
+ float ka[4];
+ float ks[4];
+ ASSERT(atm && tetra && k);
+ ASSERT(iband < sck_get_bands_count(atm->gas.ck));
+ ASSERT(!SUVM_PRIMITIVE_NONE(tetra));
+ ASSERT(tetra->nvertices == 4);
+ ASSERT(tetra->indices[0] < sck_get_nodes_count(atm->gas.ck));
+ ASSERT(tetra->indices[1] < sck_get_nodes_count(atm->gas.ck));
+ ASSERT(tetra->indices[2] < sck_get_nodes_count(atm->gas.ck));
+ ASSERT(tetra->indices[3] < sck_get_nodes_count(atm->gas.ck));
+
+ /* Retrieve the band */
+ SCK(get_band(atm->gas.ck, iband, &band));
+ ASSERT(iquad_pt < band.quad_pts_count);
+
+ /* Get the absorption coefficients */
+ if(radcoef == RNATM_RADCOEF_Ka || radcoef == RNATM_RADCOEF_Kext) {
+ /* Retrieve the quadrature point */
+ struct sck_quad_pt quad;
+ SCK(band_get_quad_pt(&band, iquad_pt, &quad));
+ ka[0] = quad.ka_list[tetra->indices[0]];
+ ka[1] = quad.ka_list[tetra->indices[1]];
+ ka[2] = quad.ka_list[tetra->indices[2]];
+ ka[3] = quad.ka_list[tetra->indices[3]];
+ }
+
+ /* Get the scattering coefficients */
+ if(radcoef == RNATM_RADCOEF_Ks || radcoef == RNATM_RADCOEF_Kext) {
+ ks[0] = band.ks_list[tetra->indices[0]];
+ ks[1] = band.ks_list[tetra->indices[1]];
+ ks[2] = band.ks_list[tetra->indices[2]];
+ ks[3] = band.ks_list[tetra->indices[3]];
+ }
+
+ switch(radcoef) {
+ case RNATM_RADCOEF_Ka: f4_set(k, ka); break;
+ case RNATM_RADCOEF_Ks: f4_set(k, ks); break;
+ case RNATM_RADCOEF_Kext: f4_add(k, ka, ks); break;
+ default: FATAL("Unreachable code\n");
+ }
+
+ return 1;
+}
+
+static int
+tetra_get_radcoef_aerosol
+ (const struct rnatm* atm,
+ const struct suvm_primitive* tetra,
+ const size_t iband,
+ const struct aerosol* aerosol,
+ const enum rnatm_radcoef radcoef,
+ float k[4])
+{
+ struct sck_band gas_band;
+ struct sars_band ars_band;
+ double gas_spectral_range[2]; /* In nm */
+ size_t ars_ibands[2];
+ float ka[4], ks[4];
+ ASSERT(atm && aerosol && tetra && k);
+
+ ASSERT(!SUVM_PRIMITIVE_NONE(tetra));
+ ASSERT(tetra->nvertices == 4);
+ ASSERT(tetra->indices[0] < sars_get_nodes_count(aerosol->sars));
+ ASSERT(tetra->indices[1] < sars_get_nodes_count(aerosol->sars));
+ ASSERT(tetra->indices[2] < sars_get_nodes_count(aerosol->sars));
+ ASSERT(tetra->indices[3] < sars_get_nodes_count(aerosol->sars));
+
+ /* Look for the aerosol bands covered by the gas band */
+ SCK(get_band(atm->gas.ck, iband, &gas_band));
+ gas_spectral_range[0] = gas_band.lower;
+ gas_spectral_range[1] = nextafter(gas_band.upper, 0); /* inclusive */
+ SARS(find_bands(aerosol->sars, gas_spectral_range, ars_ibands));
+
+ /* No aerosol band is overlaid by the gas band */
+ if(ars_ibands[0] > ars_ibands[1]) {
+ k[0] = 0;
+ k[1] = 0;
+ k[2] = 0;
+ k[3] = 0;
+ return 0;
+ }
+
+ SARS(get_band(aerosol->sars, ars_ibands[0], &ars_band));
+
+ /* The gas band is included in an aerosol band: use the aerosol radiative
+ * coefficients */
+ if(ars_ibands[0] == ars_ibands[1]
+ && ars_band.lower <= gas_band.lower
+ && ars_band.upper >= gas_band.upper) {
+
+ /* Get the absorption coefficient */
+ if(radcoef == RNATM_RADCOEF_Ka || radcoef == RNATM_RADCOEF_Kext) {
+ ka[0] = sars_band_get_ka(&ars_band, tetra->indices[0]);
+ ka[1] = sars_band_get_ka(&ars_band, tetra->indices[1]);
+ ka[2] = sars_band_get_ka(&ars_band, tetra->indices[2]);
+ ka[3] = sars_band_get_ka(&ars_band, tetra->indices[3]);
+ }
+
+ /* Get the scattering coefficient */
+ if(radcoef == RNATM_RADCOEF_Ks || radcoef == RNATM_RADCOEF_Kext) {
+ ks[0] = sars_band_get_ks(&ars_band, tetra->indices[0]);
+ ks[1] = sars_band_get_ks(&ars_band, tetra->indices[1]);
+ ks[2] = sars_band_get_ks(&ars_band, tetra->indices[2]);
+ ks[3] = sars_band_get_ks(&ars_band, tetra->indices[3]);
+ }
+
+ /* The gas band overlaid N aerosol bands (N >= 1) */
+ } else {
+ float tau_ka[4] = {0, 0, 0, 0};
+ float tau_ks[4] = {0, 0, 0, 0};
+ double lambda_min;
+ double lambda_max;
+ float rcp_gas_band_len;
+ size_t iars_band;
+
+ FOR_EACH(iars_band, ars_ibands[0], ars_ibands[1]+1) {
+ float lambda_len;
+
+ SARS(get_band(aerosol->sars, iars_band, &ars_band));
+ lambda_min = MMAX(gas_band.lower, ars_band.lower);
+ lambda_max = MMIN(gas_band.upper, ars_band.upper); /* exclusive */
+ lambda_max = nextafter(lambda_max, 0); /* inclusive */
+
+ /* Shrink lambda_<min|max> to the spectral range */
+ lambda_min = MMAX(lambda_min, atm->spectral_range[0]);
+ lambda_max = MMIN(lambda_max, atm->spectral_range[1]);
+ lambda_len = (float)(lambda_max - lambda_min);
+ ASSERT(lambda_len > 0);
+
+ /* Get the absorption coefficient */
+ if(radcoef == RNATM_RADCOEF_Ka || radcoef == RNATM_RADCOEF_Kext) {
+ tau_ka[0] += sars_band_get_ka(&ars_band, tetra->indices[0]) * lambda_len;
+ tau_ka[1] += sars_band_get_ka(&ars_band, tetra->indices[1]) * lambda_len;
+ tau_ka[2] += sars_band_get_ka(&ars_band, tetra->indices[2]) * lambda_len;
+ tau_ka[3] += sars_band_get_ka(&ars_band, tetra->indices[3]) * lambda_len;
+ }
+
+ /* Get the scattering coefficient */
+ if(radcoef == RNATM_RADCOEF_Ks || radcoef == RNATM_RADCOEF_Kext) {
+ tau_ks[0] += sars_band_get_ks(&ars_band, tetra->indices[0]) * lambda_len;
+ tau_ks[1] += sars_band_get_ks(&ars_band, tetra->indices[1]) * lambda_len;
+ tau_ks[2] += sars_band_get_ks(&ars_band, tetra->indices[2]) * lambda_len;
+ tau_ks[3] += sars_band_get_ks(&ars_band, tetra->indices[3]) * lambda_len;
+ }
+ }
+
+ /* Compute the radiative coefficients of the tetrahedron */
+ lambda_min = gas_band.lower;
+ lambda_max = nextafter(gas_band.upper, 0); /* inclusive */
+ rcp_gas_band_len = 1.f/(float)(lambda_max - lambda_min);
+ f4_mulf(ks, tau_ks, rcp_gas_band_len);
+ f4_mulf(ka, tau_ka, rcp_gas_band_len);
+ }
+
+ switch(radcoef) {
+ case RNATM_RADCOEF_Ka: f4_set(k, ka); break;
+ case RNATM_RADCOEF_Ks: f4_set(k, ks); break;
+ case RNATM_RADCOEF_Kext: f4_add(k, ka, ks); break;
+ default: FATAL("Unreachable code\n");
+ }
+
+ return 1;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+int
+tetra_get_radcoef
+ (const struct rnatm* atm,
+ const struct suvm_primitive* tetra,
+ const size_t component,
+ const size_t iband,
+ const size_t iquad_pt,
+ const enum rnatm_radcoef radcoef,
+ float k[4])
+{
+ ASSERT(atm);
+
+ if(SUVM_PRIMITIVE_NONE(tetra)) {
+ k[0] = 0;
+ k[1] = 0;
+ k[2] = 0;
+ k[3] = 0;
+ return 0;
+ }
+
+ if(component == RNATM_GAS) {
+ return tetra_get_radcoef_gas(atm, tetra, iband, iquad_pt, radcoef, k);
+ } else {
+ const struct aerosol* aerosol = NULL;
+ ASSERT(component < darray_aerosol_size_get(&atm->aerosols));
+ aerosol = darray_aerosol_cdata_get(&atm->aerosols) + component;
+ return tetra_get_radcoef_aerosol(atm, tetra, iband, aerosol, radcoef, k);
+ }
+}
