commit cc23f79b77a6aec67d501013396c3ce941d6c975
parent 6afe4cf69bc00ebea1c7ed2a10dbca428964752d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 21 Jan 2021 11:44:47 +0100
Add the atrstm_fetch_radcoefs function
Rename property/optprop in radcoef
Diffstat:
14 files changed, 1068 insertions(+), 437 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -56,16 +56,16 @@ set(ATRSTM_FILES_SRC
atrstm.c
atrstm_cache.c
atrstm_log.c
- atrstm_optprops.c
atrstm_partition.c
+ atrstm_radcoefs.c
atrstm_setup_octrees.c
atrstm_setup_uvm.c)
set(ATRSTM_FILES_INC
atrstm_c.h
atrstm_cache.h
atrstm_log.h
- atrstm_optprops.h
atrstm_partition.h
+ atrstm_radcoefs.h
atrstm_setup_octrees.h)
set(ATRSTM_FILES_INC_API
@@ -110,7 +110,7 @@ if(NOT NO_TEST)
endfunction()
build_test(test_atrstm)
- new_test(test_atrstm_optprops)
+ new_test(test_atrstm_radcoefs)
endif()
################################################################################
diff --git a/src/atrstm.c b/src/atrstm.c
@@ -21,7 +21,6 @@
#include "atrstm_log.h"
#include "atrstm_setup_octrees.h"
-#include <astoria/atrri.h>
#include <astoria/atrtp.h>
#include <star/suvm.h>
@@ -166,6 +165,21 @@ atrstm_create
res = atrri_load(atrstm->atrri, args->atrri_filename);
if(res != RES_OK) goto error;
+ /* In shortwave, pre-fetched the refractive index */
+ if(atrstm->spectral_type != ATRSTM_SPECTRAL_SW) {
+ atrstm->refract_id = ATRRI_REFRACTIVE_INDEX_NULL;
+ } else {
+ const double wlen = atrstm->wlen_range[0];
+ ASSERT(atrstm->wlen_range[0] == atrstm->wlen_range[1]);
+ res = atrri_fetch_refractive_index(atrstm->atrri, wlen, &atrstm->refract_id);
+ if(res != RES_OK) {
+ log_err(atrstm,
+ "Could not fetch the refractive index of the shortwave wavelength %g "
+ "-- %s.\n", wlen, res_to_cstr(res));
+ goto error;
+ }
+ }
+
/* Create the Star-UnstructuredVolumetricMesh device */
res = suvm_device_create
(atrstm->logger, atrstm->allocator, atrstm->verbose, &atrstm->suvm);
@@ -222,6 +236,13 @@ atrstm_ref_put(struct atrstm* atrstm)
return RES_OK;
}
+const char*
+atrstm_get_name(const struct atrstm* atrstm)
+{
+ ASSERT(atrstm);
+ return str_cget(&atrstm->name);
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
diff --git a/src/atrstm.h b/src/atrstm.h
@@ -38,11 +38,22 @@
#define ATRSTM(Func) atrstm_ ## Func
#endif
-enum atrstm_property {
- ATRSTM_Ks, /* Scattering coefficient */
- ATRSTM_Ka, /* Absorption coefficient */
- ATRSTM_Kext, /* Extinction coefficient = Ks + Ka */
- ATRSTM_PROPS_COUNT__
+/* Optical properties */
+enum atrstm_radcoef {
+ ATRSTM_RADCOEF_Ks, /* Scattering coefficient */
+ ATRSTM_RADCOEF_Ka, /* Absorption coefficient */
+ ATRSTM_RADCOEF_Kext, /* Extinction coefficient = Ks + Ka */
+ ATRSTM_RADCOEFS_COUNT__
+};
+
+enum atrstm_radcoef_flag {
+ ATRSTM_RADCOEF_FLAG_Ks = BIT(ATRSTM_RADCOEF_Ks),
+ ATRSTM_RADCOEF_FLAG_Ka = BIT(ATRSTM_RADCOEF_Ka),
+ ATRSTM_RADCOEF_FLAG_Kext = BIT(ATRSTM_RADCOEF_Kext),
+ ATRSTM_RADCOEFS_MASK_ALL =
+ ATRSTM_RADCOEF_FLAG_Ks
+ | ATRSTM_RADCOEF_FLAG_Ka
+ | ATRSTM_RADCOEF_FLAG_Kext
};
/* List of medium components */
@@ -52,6 +63,12 @@ enum atrstm_component {
ATRSTM_CPNTS_COUNT__
};
+enum atrstm_component_flag {
+ ATRSTM_CPNT_FLAG_SOOT = BIT(ATRSTM_CPNT_SOOT),
+ ATRSTM_CPNT_FLAG_GAS = BIT(ATRSTM_CPNT_GAS),
+ ATRSTM_CPNTS_MASK_ALL = ATRSTM_CPNT_FLAG_SOOT | ATRSTM_CPNT_FLAG_GAS
+};
+
enum atrstm_svx_op {
ATRSTM_SVX_MIN,
ATRSTM_SVX_MAX,
@@ -106,6 +123,35 @@ struct atrstm_args {
}
static const struct atrstm_args ATRSTM_ARGS_DEFAULT = ATRSTM_ARGS_DEFAULT__;
+struct atrstm_fetch_radcoefs_args {
+ double pos[3]; /* Position to query */
+ size_t iband; /* Spectral band index. Not use in shortwave */
+ size_t iquad; /* Quadrature point index. Not use in shortwave */
+ double wavelength; /* In nm */
+
+ int radcoefs_mask; /* Combination of atrstm_radcoef_flag */
+ int components_mask; /* Combination of atrstm_component_flag */
+
+ /* For debug: assert if the fetched property is not in [k_min, k_max] */
+ double k_min[ATRSTM_RADCOEFS_COUNT__];
+ double k_max[ATRSTM_RADCOEFS_COUNT__];
+};
+#define ATRSTM_FETCH_RADCOEFS_ARGS_DEFAULT__ { \
+ {0, 0, 0}, /* Position */ \
+ SIZE_MAX, /* Index of the spectral band */ \
+ SIZE_MAX, /* Index of the quadrature point */ \
+ DBL_MAX, /* Wavelength */ \
+ \
+ ATRSTM_RADCOEFS_MASK_ALL, /* Mask of radiative properties to fetch */ \
+ ATRSTM_CPNTS_MASK_ALL, /* Mask of comonent to handle */ \
+ \
+ /* For debug */ \
+ {-DBL_MAX,-DBL_MAX,-DBL_MAX}, /* Kmin */ \
+ { DBL_MAX, DBL_MAX, DBL_MAX}, /* Kmax */ \
+}
+static const struct atrstm_fetch_radcoefs_args
+ATRSTM_FETCH_RADCOEFS_ARGS_DEFAULT = ATRSTM_FETCH_RADCOEFS_ARGS_DEFAULT__;
+
/* Forward declaration of extern data types */
struct logger;
struct mem_allocator;
@@ -129,11 +175,20 @@ ATRSTM_API res_T
atrstm_ref_get
(struct atrstm* atrstm);
-
ATRSTM_API res_T
atrstm_ref_put
(struct atrstm* atrstm);
+ATRSTM_API const char*
+atrstm_get_name
+ (const struct atrstm* atrstm);
+
+ATRSTM_API res_T
+atrstm_fetch_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args,
+ double radcoefs[ATRSTM_RADCOEFS_COUNT__]);
+
END_DECLS
#endif /* ATRSTM_H */
diff --git a/src/atrstm_c.h b/src/atrstm_c.h
@@ -18,6 +18,8 @@
#include "atrstm.h"
+#include <astoria/atrri.h>
+
#include <rsys/logger.h>
#include <rsys/mem_allocator.h>
#include <rsys/ref_count.h>
@@ -34,7 +36,7 @@ struct atrstm {
struct atrri* atrri; /* Handle toward the Astoria Refractive Indices */
/* Unstructured volumetric mesh */
- struct suvm_device* suvm;
+ struct suvm_device* suvm;
struct suvm_volume* volume;
struct svx_device* svx;
@@ -42,7 +44,7 @@ struct atrstm {
unsigned nthreads; /* #nthreads */
- struct str name; /* Name of the gas mixture */
+ struct str name; /* Name of the semi-transparent medium */
double gyration_radius_prefactor;
double fractal_dimension;
@@ -51,7 +53,10 @@ struct atrstm {
double wlen_range[2]; /* Spectral range in nm */
unsigned grid_max_definition[3];
double optical_thickness;
-
+
+ /* Pre-fetched refractive id in shortwave */
+ struct atrri_refractive_index refract_id;
+
struct cache* cache; /* Cache of the SVX data structures */
struct mem_allocator* allocator;
diff --git a/src/atrstm_optprops.c b/src/atrstm_optprops.c
@@ -1,159 +0,0 @@
-/* Copyright (C) 2020 CNRS
- *
- * 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/>. */
-
-#include "atrstm_c.h"
-#include "atrstm_log.h"
-#include "atrstm_optprops.h"
-
-#include <astoria/atrri.h>
-#include <astoria/atrtp.h>
-
-#include <star/suvm.h>
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-primitive_compute_optprops
- (const struct atrstm* atrstm,
- const struct atrri_refractive_index* refract_id,
- const struct suvm_primitive* prim,
- struct optprops optprops[])
-{
- struct optprops_compute_args args = OPTPROPS_COMPUTE_ARGS_NULL;
- struct atrtp_desc desc = ATRTP_DESC_NULL;
- size_t inode;
- res_T res = RES_OK;
- ASSERT(atrstm && prim && prim->nvertices == 4 && refract_id && optprops);
-
- res = atrtp_get_desc(atrstm->atrtp, &desc);
- if(res != RES_OK) goto error;
-
- /* Setup the constant input arguments of the optical properties computation
- * routine */
- args.lambda = refract_id->wavelength;
- args.n = refract_id->n;
- args.kappa = refract_id->kappa;
- args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
- args.fractal_dimension = atrstm->fractal_dimension;
-
- FOR_EACH(inode, 0, 4) {
- const double* node;
-
- /* Fetch the thermodynamic properties of the node */
- node = atrtp_desc_get_node_properties(&desc, prim->indices[inode]);
-
- /* Setup the per node input arguments of the optical properties computation
- * routine */
- args.soot_volumic_fraction =
- node[ATRTP_SOOT_VOLFRAC];
- args.soot_primary_particles_count =
- node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT];
- args.soot_primary_particles_diameter =
- node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER];
-
- /* Compute the node optical properties */
- optprops_compute(&optprops[inode], &args);
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-res_T
-dump_optprops
- (const struct atrstm* atrstm,
- const struct atrri_refractive_index* refract_id,
- FILE* stream)
-{
- struct optprops_compute_args args = OPTPROPS_COMPUTE_ARGS_NULL;
- struct atrtp_desc desc = ATRTP_DESC_NULL;
- struct optprops props_min;
- struct optprops props_max;
- size_t inode;
-
- res_T res = RES_OK;
- ASSERT(atrstm && stream);
-
- res = atrtp_get_desc(atrstm->atrtp, &desc);
- if(res != RES_OK) goto error;
-
- /* Setup the constant input params of the optical properties computation */
- args.lambda = refract_id->wavelength;
- args.n = refract_id->n;
- args.kappa = refract_id->kappa;
- args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
- args.fractal_dimension = atrstm->fractal_dimension;
-
- props_min.ka = DBL_MAX;
- props_max.ka =-DBL_MAX;
- props_min.ks = DBL_MAX;
- props_max.ks =-DBL_MAX;
- props_min.kext = DBL_MAX;
- props_max.kext =-DBL_MAX;
-
- #define FPRINTF(Text, Args) { \
- const int n = fprintf(stream, Text COMMA_##Args LIST_##Args); \
- if(n < 0) { \
- log_err(atrstm, "Error writing optical properties.\n"); \
- res = RES_IO_ERR; \
- goto error; \
- } \
- } (void)0
-
- FPRINTF("# Ka Ks Kext\n", ARG0());
-
- FOR_EACH(inode, 0, desc.nnodes) {
- struct optprops props;
- const double* node;
-
- /* Fetch the thermodynamic properties of the node */
- node = atrtp_desc_get_node_properties(&desc, inode);
-
- /* Setup the per node input args of the optical properties computation */
- args.soot_volumic_fraction =
- node[ATRTP_SOOT_VOLFRAC];
- args.soot_primary_particles_count =
- node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT];
- args.soot_primary_particles_diameter =
- node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER];
-
- /* Compute the node optical properties */
- optprops_compute(&props, &args);
-
- FPRINTF("%g %g %g\n", ARG3(props.ka, props.ks, props.kext));
-
- props_min.ka = MMIN(props_min.ka, props.ka);
- props_max.ka = MMAX(props_max.ka, props.ka);
- props_min.ks = MMIN(props_min.ks, props.ks);
- props_max.ks = MMAX(props_max.ks, props.ks);
- props_min.kext = MMIN(props_min.kext, props.kext);
- props_max.kext = MMAX(props_max.kext, props.kext);
-
- }
-
- FPRINTF("# Bounds = [%g %g %g], [%g, %g, %g]\n", ARG6
- (props_min.ka, props_min.ks, props_min.kext,
- props_max.ka, props_max.ks, props_max.kext));
-
- #undef FPRINTF
-
-exit:
- return res;
-error:
- goto exit;
-}
diff --git a/src/atrstm_optprops.h b/src/atrstm_optprops.h
@@ -1,165 +0,0 @@
-/* Copyright (C) 2020 CNRS
- *
- * 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/>. */
-
-#ifndef ATRSTM_OPTPROPS_H
-#define ATRSTM_OPTPROPS_H
-
-#include <rsys/math.h>
-#include <rsys/rsys.h>
-
-/* Optical properties */
-struct optprops {
- double ka; /* Absorption coefficient */
- double ks; /* Scattering coefficient */
- double kext; /* Extinction coefficient */
-};
-#define OPTPROPS_NULL__ {0}
-static const struct optprops OPTPROPS_NULL = OPTPROPS_NULL__;
-
-struct optprops_compute_args {
- double lambda; /* wavelength [nm] */
- double n; /* Refractive index (real component) */
- double kappa; /* Refractive index (imaginary component) */
- double gyration_radius_prefactor;
- double fractal_dimension;
- double soot_volumic_fraction; /* In m^3(soot) / m^3 */
- double soot_primary_particles_count;
- double soot_primary_particles_diameter; /* In nm */
-};
-#define OPTPROPS_COMPUTE_ARGS_NULL__ {0}
-static const struct optprops_compute_args OPTPROPS_COMPUTE_ARGS_NULL =
- OPTPROPS_COMPUTE_ARGS_NULL__;
-
-/* Forward declarations */
-struct atrri_refractive_index;
-struct atrstm;
-struct suvm_primitive;
-
-extern LOCAL_SYM res_T
-primitive_compute_optprops
- (const struct atrstm* atrstm,
- const struct atrri_refractive_index* refract_id,
- const struct suvm_primitive* prim,
- struct optprops optprops[]); /* Per primitive node optical properties */
-
-/* Write per node optical properties regarding the submitted refracted index */
-extern LOCAL_SYM res_T
-dump_optprops
- (const struct atrstm* atrstm,
- const struct atrri_refractive_index* refract_id,
- FILE* stream);
-
-static INLINE res_T
-primitive_compute_optprops_range
- (const struct atrstm* atrstm,
- const struct atrri_refractive_index* refract_id,
- const struct suvm_primitive* prim,
- struct optprops* optprops_min,
- struct optprops* optprops_max)
-{
- struct optprops props[4];
- res_T res = RES_OK;
- ASSERT(optprops_min && optprops_max);
-
- res = primitive_compute_optprops(atrstm, refract_id, prim, props);
- if(res != RES_OK) return res;
-
- optprops_min->ka = MMIN
- (MMIN(props[0].ka, props[1].ka),
- MMIN(props[2].ka, props[3].ka));
- optprops_min->ks = MMIN
- (MMIN(props[0].ks, props[1].ks),
- MMIN(props[2].ks, props[3].ks));
- optprops_min->kext = MMIN
- (MMIN(props[0].kext, props[1].kext),
- MMIN(props[2].kext, props[3].kext));
-
- optprops_max->ka = MMAX
- (MMAX(props[0].ka, props[1].ka),
- MMAX(props[2].ka, props[3].ka));
- optprops_max->ks = MMAX
- (MMAX(props[0].ks, props[1].ks),
- MMAX(props[2].ks, props[3].ks));
- optprops_max->kext = MMAX
- (MMAX(props[0].kext, props[1].kext),
- MMAX(props[2].kext, props[3].kext));
-
- return RES_OK;
-}
-
-/* Compute the optical properties of the semi transparent medium as described
- * in "Effects of multiple scattering on radiative properties of soot fractal
- * aggregates" J. Yon et al., Journal of Quantitative Spectroscopy and
- * Radiative Transfer Vol 133, pp. 374-381, 2014 */
-static INLINE void
-optprops_compute
- (struct optprops* optprops,
- const struct optprops_compute_args* args)
-{
- /* Input variables */
- const double lambda = args->lambda; /* [nm] */
- const double n = args->n;
- const double kappa = args->kappa;
- const double Fv = args->soot_volumic_fraction; /* [nm^3(soot)/nm] */
- const double Np = args->soot_primary_particles_count;
- const double Dp = args->soot_primary_particles_diameter; /* [nm] */
- const double kf = args->gyration_radius_prefactor;
- const double Df = args->fractal_dimension;
-
- if(Np == 0 || Fv == 0) {
- optprops->ka = 0;
- optprops->ks = 0;
- optprops->kext = 0;
- } else {
- /* Precomputed values */
- const double n2 = n*n;
- const double kappa2 = kappa*kappa;
- const double four_n2_kappa2 = 4*n2*kappa2;
- const double xp = (PI*Dp) / lambda;
- const double xp3 = xp*xp*xp;
- const double k = (2*PI) / lambda;
- const double k2 = k*k;
-
- /* E(m), m = n + i*kappa */
- const double tmp0 = (n2 - kappa2 + 2);
- const double denom = tmp0*tmp0 + four_n2_kappa2;
- const double Em = (6*n*kappa) / denom;
-
- /* F(m), m = n + i*kappa */
- const double tmp1 = ((n2-kappa2 - 1)*(n2-kappa2+2) + four_n2_kappa2)/denom;
- const double Fm = tmp1*tmp1 + Em*Em;
-
- /* Gyration radius */
- const double Rg = 0.5 * Dp * pow(Np/kf, 1.0/Df); /* [nm] */
- const double Rg2 = Rg*Rg;
- const double g = pow(1 + 4*k2*Rg2/(3*Df), -Df*0.5);
-
- /* Cross sections, [nm^3/aggrefate] */
- const double Cabs = Np * (4*PI*xp3)/(k*k) * Em;
- const double Csca = Np*Np * (8*PI*xp3*xp3) / (3*k2) * Fm * g;
-
- /* Optical properties */
- const double Va = (Np*PI*Dp*Dp*Dp) / 6; /* [nm^3] */
- const double rho = Fv / Va; /* [#aggregates / nm^3] */
- const double tmp2 = 1.e-9 * rho;
-
- ASSERT(denom != 0 && Df != 0 && Dp != 0);
- optprops->ka = tmp2 * Cabs;
- optprops->ks = tmp2 * Csca;
- optprops->kext = optprops->ka + optprops->ks;
- }
-}
-
-#endif /* ATRSTM_OPTPROPS_H */
diff --git a/src/atrstm_radcoefs.c b/src/atrstm_radcoefs.c
@@ -0,0 +1,271 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * 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/>. */
+
+#include "atrstm_c.h"
+#include "atrstm_log.h"
+#include "atrstm_radcoefs.h"
+
+#include <astoria/atrri.h>
+#include <astoria/atrtp.h>
+
+#include <star/suvm.h>
+
+#include <rsys/cstr.h>
+
+#include <string.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static INLINE int
+check_fetch_radcoefs_args
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args)
+{
+ int i;
+ ASSERT(atrstm);
+
+ if(!args
+ || args->wavelength < atrstm->wlen_range[0]
+ || args->wavelength > atrstm->wlen_range[1]
+ || !(args->radcoefs_mask & ATRSTM_RADCOEFS_MASK_ALL)
+ || !(args->components_mask & ATRSTM_CPNTS_MASK_ALL))
+ return 0;
+
+ FOR_EACH(i, 0, ATRSTM_RADCOEFS_COUNT__) {
+ if(args->k_min[i] > args->k_max[i]) return 0;
+ }
+
+ return 1;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+atrstm_fetch_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args,
+ double radcoefs[ATRSTM_RADCOEFS_COUNT__])
+{
+ struct radcoefs prim_radcoefs[4];
+ double bcoords[4];
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ res_T res = RES_OK;
+
+ if(!atrstm || !check_fetch_radcoefs_args(atrstm, args) || !radcoefs) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ memset(radcoefs, 0, sizeof(double[ATRSTM_RADCOEFS_COUNT__]));
+
+ /* Find the primitive into which pos lies */
+ res = suvm_volume_at(atrstm->volume, args->pos, &prim, bcoords);
+ if(res != RES_OK) {
+ log_err(atrstm,
+ "Error accessing the volumetric mesh at {%g, %g, %g} -- %s.\n",
+ SPLIT3(args->pos), res_to_cstr(res));
+ goto error;
+ }
+
+ /* No primitive found */
+ if(SUVM_PRIMITIVE_NONE(&prim)) {
+ log_warn(atrstm,
+ "%s: the position {%g, %g, %g} does not belongs to the semi-transparent "
+ "medium.\n", FUNC_NAME, SPLIT3(args->pos));
+ goto exit;
+ }
+
+ /* Compute the radiative properties of the primitive nodes */
+ res = primitive_compute_radcoefs
+ (atrstm, &atrstm->refract_id, &prim, prim_radcoefs);
+ if(res != RES_OK) {
+ log_err(atrstm,
+ "Error computing the radiative properties for the primitive %lu -- %s.\n",
+ (unsigned long)prim.iprim,
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Linearly interpolate the node radiative properties */
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) {
+ radcoefs[ATRSTM_RADCOEF_Ka] =
+ prim_radcoefs[0].ka * bcoords[0]
+ + prim_radcoefs[1].ka * bcoords[1]
+ + prim_radcoefs[2].ka * bcoords[2]
+ + prim_radcoefs[3].ka * bcoords[3];
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] >= args->k_min[ATRSTM_RADCOEF_Ka]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] <= args->k_max[ATRSTM_RADCOEF_Ka]);
+ }
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) {
+ radcoefs[ATRSTM_RADCOEF_Ks] =
+ prim_radcoefs[0].ks * bcoords[0]
+ + prim_radcoefs[1].ks * bcoords[1]
+ + prim_radcoefs[2].ks * bcoords[2]
+ + prim_radcoefs[3].ks * bcoords[3];
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] >= args->k_min[ATRSTM_RADCOEF_Ks]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] <= args->k_max[ATRSTM_RADCOEF_Ks]);
+ }
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) {
+ radcoefs[ATRSTM_RADCOEF_Kext] =
+ prim_radcoefs[0].kext * bcoords[0]
+ + prim_radcoefs[1].kext * bcoords[1]
+ + prim_radcoefs[2].kext * bcoords[2]
+ + prim_radcoefs[3].kext * bcoords[3];
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] >= args->k_min[ATRSTM_RADCOEF_Kext]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] <= args->k_max[ATRSTM_RADCOEF_Kext]);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+primitive_compute_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct radcoefs radcoefs[])
+{
+ struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL;
+ struct atrtp_desc desc = ATRTP_DESC_NULL;
+ size_t inode;
+ res_T res = RES_OK;
+ ASSERT(atrstm && prim && prim->nvertices == 4 && refract_id && radcoefs);
+
+ res = atrtp_get_desc(atrstm->atrtp, &desc);
+ if(res != RES_OK) goto error;
+
+ /* Setup the constant input arguments of the optical properties computation
+ * routine */
+ args.lambda = refract_id->wavelength;
+ args.n = refract_id->n;
+ args.kappa = refract_id->kappa;
+ args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
+ args.fractal_dimension = atrstm->fractal_dimension;
+
+ FOR_EACH(inode, 0, 4) {
+ const double* node;
+
+ /* Fetch the thermodynamic properties of the node */
+ node = atrtp_desc_get_node_properties(&desc, prim->indices[inode]);
+
+ /* Setup the per node input arguments of the optical properties computation
+ * routine */
+ args.soot_volumic_fraction =
+ node[ATRTP_SOOT_VOLFRAC];
+ args.soot_primary_particles_count =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT];
+ args.soot_primary_particles_diameter =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER];
+
+ /* Compute the node optical properties */
+ radcoefs_compute(&radcoefs[inode], &args);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+dump_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ FILE* stream)
+{
+ struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL;
+ struct atrtp_desc desc = ATRTP_DESC_NULL;
+ struct radcoefs props_min;
+ struct radcoefs props_max;
+ size_t inode;
+
+ res_T res = RES_OK;
+ ASSERT(atrstm && stream);
+
+ res = atrtp_get_desc(atrstm->atrtp, &desc);
+ if(res != RES_OK) goto error;
+
+ /* Setup the constant input params of the optical properties computation */
+ args.lambda = refract_id->wavelength;
+ args.n = refract_id->n;
+ args.kappa = refract_id->kappa;
+ args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
+ args.fractal_dimension = atrstm->fractal_dimension;
+
+ props_min.ka = DBL_MAX;
+ props_max.ka =-DBL_MAX;
+ props_min.ks = DBL_MAX;
+ props_max.ks =-DBL_MAX;
+ props_min.kext = DBL_MAX;
+ props_max.kext =-DBL_MAX;
+
+ #define FPRINTF(Text, Args) { \
+ const int n = fprintf(stream, Text COMMA_##Args LIST_##Args); \
+ if(n < 0) { \
+ log_err(atrstm, "Error writing optical properties.\n"); \
+ res = RES_IO_ERR; \
+ goto error; \
+ } \
+ } (void)0
+
+ FPRINTF("# Ka Ks Kext\n", ARG0());
+
+ FOR_EACH(inode, 0, desc.nnodes) {
+ struct radcoefs props;
+ const double* node;
+
+ /* Fetch the thermodynamic properties of the node */
+ node = atrtp_desc_get_node_properties(&desc, inode);
+
+ /* Setup the per node input args of the optical properties computation */
+ args.soot_volumic_fraction =
+ node[ATRTP_SOOT_VOLFRAC];
+ args.soot_primary_particles_count =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT];
+ args.soot_primary_particles_diameter =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER];
+
+ /* Compute the node optical properties */
+ radcoefs_compute(&props, &args);
+
+ FPRINTF("%g %g %g\n", ARG3(props.ka, props.ks, props.kext));
+
+ props_min.ka = MMIN(props_min.ka, props.ka);
+ props_max.ka = MMAX(props_max.ka, props.ka);
+ props_min.ks = MMIN(props_min.ks, props.ks);
+ props_max.ks = MMAX(props_max.ks, props.ks);
+ props_min.kext = MMIN(props_min.kext, props.kext);
+ props_max.kext = MMAX(props_max.kext, props.kext);
+
+ }
+
+ FPRINTF("# Bounds = [%g %g %g], [%g, %g, %g]\n", ARG6
+ (props_min.ka, props_min.ks, props_min.kext,
+ props_max.ka, props_max.ks, props_max.kext));
+
+ #undef FPRINTF
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/atrstm_radcoefs.h b/src/atrstm_radcoefs.h
@@ -0,0 +1,166 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * 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/>. */
+
+#ifndef ATRSTM_RADCOEFS_H
+#define ATRSTM_RADCOEFS_H
+
+#include <rsys/math.h>
+#include <rsys/rsys.h>
+
+/* Radiative coefficients */
+struct radcoefs {
+ double ka; /* Absorption coefficient */
+ double ks; /* Scattering coefficient */
+ double kext; /* Extinction coefficient */
+};
+#define RADCOEFS_NULL__ {0}
+static const struct radcoefs RADCOEFS_NULL = RADCOEFS_NULL__;
+
+struct radcoefs_compute_args {
+ double lambda; /* wavelength [nm] */
+ double n; /* Refractive index (real component) */
+ double kappa; /* Refractive index (imaginary component) */
+ double gyration_radius_prefactor;
+ double fractal_dimension;
+ double soot_volumic_fraction; /* In m^3(soot) / m^3 */
+ double soot_primary_particles_count;
+ double soot_primary_particles_diameter; /* In nm */
+};
+#define RADCOEFS_COMPUTE_ARGS_NULL__ {0}
+static const struct radcoefs_compute_args RADCOEFS_COMPUTE_ARGS_NULL =
+ RADCOEFS_COMPUTE_ARGS_NULL__;
+
+/* Forward declarations */
+struct atrri_refractive_index;
+struct atrstm;
+struct suvm_primitive;
+
+extern LOCAL_SYM res_T
+primitive_compute_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct radcoefs radcoefs[]); /* Per primitive node radiative coefficients */
+
+/* Write per node radiative coefficients regarding the submitted refracted
+ * index */
+extern LOCAL_SYM res_T
+dump_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ FILE* stream);
+
+static INLINE res_T
+primitive_compute_radcoefs_range
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct radcoefs* radcoefs_min,
+ struct radcoefs* radcoefs_max)
+{
+ struct radcoefs props[4];
+ res_T res = RES_OK;
+ ASSERT(radcoefs_min && radcoefs_max);
+
+ res = primitive_compute_radcoefs(atrstm, refract_id, prim, props);
+ if(res != RES_OK) return res;
+
+ radcoefs_min->ka = MMIN
+ (MMIN(props[0].ka, props[1].ka),
+ MMIN(props[2].ka, props[3].ka));
+ radcoefs_min->ks = MMIN
+ (MMIN(props[0].ks, props[1].ks),
+ MMIN(props[2].ks, props[3].ks));
+ radcoefs_min->kext = MMIN
+ (MMIN(props[0].kext, props[1].kext),
+ MMIN(props[2].kext, props[3].kext));
+
+ radcoefs_max->ka = MMAX
+ (MMAX(props[0].ka, props[1].ka),
+ MMAX(props[2].ka, props[3].ka));
+ radcoefs_max->ks = MMAX
+ (MMAX(props[0].ks, props[1].ks),
+ MMAX(props[2].ks, props[3].ks));
+ radcoefs_max->kext = MMAX
+ (MMAX(props[0].kext, props[1].kext),
+ MMAX(props[2].kext, props[3].kext));
+
+ return RES_OK;
+}
+
+/* Compute the radiative coefficients of the semi transparent medium as described
+ * in "Effects of multiple scattering on radiative properties of soot fractal
+ * aggregates" J. Yon et al., Journal of Quantitative Spectroscopy and
+ * Radiative Transfer Vol 133, pp. 374-381, 2014 */
+static INLINE void
+radcoefs_compute
+ (struct radcoefs* radcoefs,
+ const struct radcoefs_compute_args* args)
+{
+ /* Input variables */
+ const double lambda = args->lambda; /* [nm] */
+ const double n = args->n;
+ const double kappa = args->kappa;
+ const double Fv = args->soot_volumic_fraction; /* [nm^3(soot)/nm] */
+ const double Np = args->soot_primary_particles_count;
+ const double Dp = args->soot_primary_particles_diameter; /* [nm] */
+ const double kf = args->gyration_radius_prefactor;
+ const double Df = args->fractal_dimension;
+
+ if(Np == 0 || Fv == 0) {
+ radcoefs->ka = 0;
+ radcoefs->ks = 0;
+ radcoefs->kext = 0;
+ } else {
+ /* Precomputed values */
+ const double n2 = n*n;
+ const double kappa2 = kappa*kappa;
+ const double four_n2_kappa2 = 4*n2*kappa2;
+ const double xp = (PI*Dp) / lambda;
+ const double xp3 = xp*xp*xp;
+ const double k = (2*PI) / lambda;
+ const double k2 = k*k;
+
+ /* E(m), m = n + i*kappa */
+ const double tmp0 = (n2 - kappa2 + 2);
+ const double denom = tmp0*tmp0 + four_n2_kappa2;
+ const double Em = (6*n*kappa) / denom;
+
+ /* F(m), m = n + i*kappa */
+ const double tmp1 = ((n2-kappa2 - 1)*(n2-kappa2+2) + four_n2_kappa2)/denom;
+ const double Fm = tmp1*tmp1 + Em*Em;
+
+ /* Gyration radius */
+ const double Rg = 0.5 * Dp * pow(Np/kf, 1.0/Df); /* [nm] */
+ const double Rg2 = Rg*Rg;
+ const double g = pow(1 + 4*k2*Rg2/(3*Df), -Df*0.5);
+
+ /* Cross sections, [nm^3/aggrefate] */
+ const double Cabs = Np * (4*PI*xp3)/(k*k) * Em;
+ const double Csca = Np*Np * (8*PI*xp3*xp3) / (3*k2) * Fm * g;
+
+ /* Radiative coefficients */
+ const double Va = (Np*PI*Dp*Dp*Dp) / 6; /* [nm^3] */
+ const double rho = Fv / Va; /* [#aggregates / nm^3] */
+ const double tmp2 = 1.e-9 * rho;
+
+ ASSERT(denom != 0 && Df != 0 && Dp != 0);
+ radcoefs->ka = tmp2 * Cabs;
+ radcoefs->ks = tmp2 * Csca;
+ radcoefs->kext = radcoefs->ka + radcoefs->ks;
+ }
+}
+
+#endif /* ATRSTM_RADCOEFS_H */
diff --git a/src/atrstm_radprops.c b/src/atrstm_radprops.c
@@ -0,0 +1,271 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * 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/>. */
+
+#include "atrstm_c.h"
+#include "atrstm_log.h"
+#include "atrstm_radcoefs.h"
+
+#include <astoria/atrri.h>
+#include <astoria/atrtp.h>
+
+#include <star/suvm.h>
+
+#include <rsys/cstr.h>
+
+#include <string.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static INLINE int
+check_fetch_radcoefs_args
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args)
+{
+ int i;
+ ASSERT(atrstm);
+
+ if(!args
+ || args->wavelength < atrstm->wlen_range[0]
+ || args->wavelength > atrstm->wlen_range[1]
+ || !(args->radcoefs_mask & ATRSTM_RADCOEFS_MASK_ALL)
+ || !(args->components_mask & ATRSTM_CPNTS_MASK_ALL))
+ return 0;
+
+ FOR_EACH(i, 0, ATRSTM_RADCOEFS_COUNT__) {
+ if(args->k_min[i] > args->k_max[i]) return 0;
+ }
+
+ return 1;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+atrstm_fetch_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args,
+ double radcoefs[ATRSTM_RADCOEFS_COUNT__])
+{
+ struct radcoefs prim_radcoefs[4];
+ double bcoords[4];
+ struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
+ res_T res = RES_OK;
+
+ if(!atrstm || !check_fetch_radcoefs_args(atrstm, args) || !radcoefs) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ memset(radcoefs, 0, sizeof(double[ATRSTM_RADCOEFS_COUNT__]));
+
+ /* Find the primitive into which pos lies */
+ res = suvm_volume_at(atrstm->volume, args->pos, &prim, bcoords);
+ if(res != RES_OK) {
+ log_err(atrstm,
+ "Error accessing the volumetric mesh at {%g, %g, %g} -- %s.\n",
+ SPLIT3(args->pos), res_to_cstr(res));
+ goto error;
+ }
+
+ /* No primitive found */
+ if(SUVM_PRIMITIVE_NONE(&prim)) {
+ log_warn(atrstm,
+ "%s: the position {%g, %g, %g} does not belongs to the semi-transparent "
+ "medium.\n", FUNC_NAME, SPLIT3(args->pos));
+ goto exit;
+ }
+
+ /* Compute the radiative properties of the primitive nodes */
+ res = primitive_compute_radcoefs
+ (atrstm, &atrstm->refract_id, &prim, prim_radcoefs);
+ if(res != RES_OK) {
+ log_err(atrstm,
+ "Error computing the radiative properties for the primitive %lu -- %s.\n",
+ (unsigned long)prim.iprim,
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Linearly interpolate the node radiative properties */
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) {
+ radcoefs[ATRSTM_RADCOEF_Ka] =
+ prim_radcoefs[0].ka * bcoords[0]
+ + prim_radcoefs[1].ka * bcoords[1]
+ + prim_radcoefs[2].ka * bcoords[2]
+ + prim_radcoefs[3].ka * bcoords[3];
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] >= args->k_min[ATRSTM_RADCOEF_Ka]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] <= args->k_max[ATRSTM_RADCOEF_Ka]);
+ }
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) {
+ radcoefs[ATRSTM_RADCOEF_Ks] =
+ prim_radcoefs[0].ks * bcoords[0]
+ + prim_radcoefs[1].ks * bcoords[1]
+ + prim_radcoefs[2].ks * bcoords[2]
+ + prim_radcoefs[3].ks * bcoords[3];
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] >= args->k_min[ATRSTM_RADCOEF_Ks]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] <= args->k_max[ATRSTM_RADCOEF_Ks]);
+ }
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) {
+ radcoefs[ATRSTM_RADCOEF_Kext] =
+ prim_radcoefs[0].kext * bcoords[0]
+ + prim_radcoefs[1].kext * bcoords[1]
+ + prim_radcoefs[2].kext * bcoords[2]
+ + prim_radcoefs[3].kext * bcoords[3];
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] >= args->k_min[ATRSTM_RADCOEF_Kext]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] <= args->k_max[ATRSTM_RADCOEF_Kext]);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+primitive_compute_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct radcoefs radcoefs[])
+{
+ struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL;
+ struct atrtp_desc desc = ATRTP_DESC_NULL;
+ size_t inode;
+ res_T res = RES_OK;
+ ASSERT(atrstm && prim && prim->nvertices == 4 && refract_id && radcoefs);
+
+ res = atrtp_get_desc(atrstm->atrtp, &desc);
+ if(res != RES_OK) goto error;
+
+ /* Setup the constant input arguments of the optical properties computation
+ * routine */
+ args.lambda = refract_id->wavelength;
+ args.n = refract_id->n;
+ args.kappa = refract_id->kappa;
+ args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
+ args.fractal_dimension = atrstm->fractal_dimension;
+
+ FOR_EACH(inode, 0, 4) {
+ const double* node;
+
+ /* Fetch the thermodynamic properties of the node */
+ node = atrtp_desc_get_node_properties(&desc, prim->indices[inode]);
+
+ /* Setup the per node input arguments of the optical properties computation
+ * routine */
+ args.soot_volumic_fraction =
+ node[ATRTP_SOOT_VOLFRAC];
+ args.soot_primary_particles_count =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT];
+ args.soot_primary_particles_diameter =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER];
+
+ /* Compute the node optical properties */
+ radcoefs_compute(&radcoefs[inode], &args);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+dump_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ FILE* stream)
+{
+ struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL;
+ struct atrtp_desc desc = ATRTP_DESC_NULL;
+ struct radcoefs props_min;
+ struct radcoefs props_max;
+ size_t inode;
+
+ res_T res = RES_OK;
+ ASSERT(atrstm && stream);
+
+ res = atrtp_get_desc(atrstm->atrtp, &desc);
+ if(res != RES_OK) goto error;
+
+ /* Setup the constant input params of the optical properties computation */
+ args.lambda = refract_id->wavelength;
+ args.n = refract_id->n;
+ args.kappa = refract_id->kappa;
+ args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
+ args.fractal_dimension = atrstm->fractal_dimension;
+
+ props_min.ka = DBL_MAX;
+ props_max.ka =-DBL_MAX;
+ props_min.ks = DBL_MAX;
+ props_max.ks =-DBL_MAX;
+ props_min.kext = DBL_MAX;
+ props_max.kext =-DBL_MAX;
+
+ #define FPRINTF(Text, Args) { \
+ const int n = fprintf(stream, Text COMMA_##Args LIST_##Args); \
+ if(n < 0) { \
+ log_err(atrstm, "Error writing optical properties.\n"); \
+ res = RES_IO_ERR; \
+ goto error; \
+ } \
+ } (void)0
+
+ FPRINTF("# Ka Ks Kext\n", ARG0());
+
+ FOR_EACH(inode, 0, desc.nnodes) {
+ struct radcoefs props;
+ const double* node;
+
+ /* Fetch the thermodynamic properties of the node */
+ node = atrtp_desc_get_node_properties(&desc, inode);
+
+ /* Setup the per node input args of the optical properties computation */
+ args.soot_volumic_fraction =
+ node[ATRTP_SOOT_VOLFRAC];
+ args.soot_primary_particles_count =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT];
+ args.soot_primary_particles_diameter =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER];
+
+ /* Compute the node optical properties */
+ radcoefs_compute(&props, &args);
+
+ FPRINTF("%g %g %g\n", ARG3(props.ka, props.ks, props.kext));
+
+ props_min.ka = MMIN(props_min.ka, props.ka);
+ props_max.ka = MMAX(props_max.ka, props.ka);
+ props_min.ks = MMIN(props_min.ks, props.ks);
+ props_max.ks = MMAX(props_max.ks, props.ks);
+ props_min.kext = MMIN(props_min.kext, props.kext);
+ props_max.kext = MMAX(props_max.kext, props.kext);
+
+ }
+
+ FPRINTF("# Bounds = [%g %g %g], [%g, %g, %g]\n", ARG6
+ (props_min.ka, props_min.ks, props_min.kext,
+ props_max.ka, props_max.ks, props_max.kext));
+
+ #undef FPRINTF
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/atrstm_radprops.h b/src/atrstm_radprops.h
@@ -0,0 +1,166 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * 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/>. */
+
+#ifndef ATRSTM_RADCOEFS_H
+#define ATRSTM_RADCOEFS_H
+
+#include <rsys/math.h>
+#include <rsys/rsys.h>
+
+/* Radiative coefficients */
+struct radcoefs {
+ double ka; /* Absorption coefficient */
+ double ks; /* Scattering coefficient */
+ double kext; /* Extinction coefficient */
+};
+#define RADCOEFS_NULL__ {0}
+static const struct radcoefs RADCOEFS_NULL = RADCOEFS_NULL__;
+
+struct radcoefs_compute_args {
+ double lambda; /* wavelength [nm] */
+ double n; /* Refractive index (real component) */
+ double kappa; /* Refractive index (imaginary component) */
+ double gyration_radius_prefactor;
+ double fractal_dimension;
+ double soot_volumic_fraction; /* In m^3(soot) / m^3 */
+ double soot_primary_particles_count;
+ double soot_primary_particles_diameter; /* In nm */
+};
+#define RADCOEFS_COMPUTE_ARGS_NULL__ {0}
+static const struct radcoefs_compute_args RADCOEFS_COMPUTE_ARGS_NULL =
+ RADCOEFS_COMPUTE_ARGS_NULL__;
+
+/* Forward declarations */
+struct atrri_refractive_index;
+struct atrstm;
+struct suvm_primitive;
+
+extern LOCAL_SYM res_T
+primitive_compute_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct radcoefs radcoefs[]); /* Per primitive node radiative coefficients */
+
+/* Write per node radiative coefficients regarding the submitted refracted
+ * index */
+extern LOCAL_SYM res_T
+dump_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ FILE* stream);
+
+static INLINE res_T
+primitive_compute_radcoefs_range
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ struct radcoefs* radcoefs_min,
+ struct radcoefs* radcoefs_max)
+{
+ struct radcoefs props[4];
+ res_T res = RES_OK;
+ ASSERT(radcoefs_min && radcoefs_max);
+
+ res = primitive_compute_radcoefs(atrstm, refract_id, prim, props);
+ if(res != RES_OK) return res;
+
+ radcoefs_min->ka = MMIN
+ (MMIN(props[0].ka, props[1].ka),
+ MMIN(props[2].ka, props[3].ka));
+ radcoefs_min->ks = MMIN
+ (MMIN(props[0].ks, props[1].ks),
+ MMIN(props[2].ks, props[3].ks));
+ radcoefs_min->kext = MMIN
+ (MMIN(props[0].kext, props[1].kext),
+ MMIN(props[2].kext, props[3].kext));
+
+ radcoefs_max->ka = MMAX
+ (MMAX(props[0].ka, props[1].ka),
+ MMAX(props[2].ka, props[3].ka));
+ radcoefs_max->ks = MMAX
+ (MMAX(props[0].ks, props[1].ks),
+ MMAX(props[2].ks, props[3].ks));
+ radcoefs_max->kext = MMAX
+ (MMAX(props[0].kext, props[1].kext),
+ MMAX(props[2].kext, props[3].kext));
+
+ return RES_OK;
+}
+
+/* Compute the radiative coefficients of the semi transparent medium as described
+ * in "Effects of multiple scattering on radiative properties of soot fractal
+ * aggregates" J. Yon et al., Journal of Quantitative Spectroscopy and
+ * Radiative Transfer Vol 133, pp. 374-381, 2014 */
+static INLINE void
+radcoefs_compute
+ (struct radcoefs* radcoefs,
+ const struct radcoefs_compute_args* args)
+{
+ /* Input variables */
+ const double lambda = args->lambda; /* [nm] */
+ const double n = args->n;
+ const double kappa = args->kappa;
+ const double Fv = args->soot_volumic_fraction; /* [nm^3(soot)/nm] */
+ const double Np = args->soot_primary_particles_count;
+ const double Dp = args->soot_primary_particles_diameter; /* [nm] */
+ const double kf = args->gyration_radius_prefactor;
+ const double Df = args->fractal_dimension;
+
+ if(Np == 0 || Fv == 0) {
+ radcoefs->ka = 0;
+ radcoefs->ks = 0;
+ radcoefs->kext = 0;
+ } else {
+ /* Precomputed values */
+ const double n2 = n*n;
+ const double kappa2 = kappa*kappa;
+ const double four_n2_kappa2 = 4*n2*kappa2;
+ const double xp = (PI*Dp) / lambda;
+ const double xp3 = xp*xp*xp;
+ const double k = (2*PI) / lambda;
+ const double k2 = k*k;
+
+ /* E(m), m = n + i*kappa */
+ const double tmp0 = (n2 - kappa2 + 2);
+ const double denom = tmp0*tmp0 + four_n2_kappa2;
+ const double Em = (6*n*kappa) / denom;
+
+ /* F(m), m = n + i*kappa */
+ const double tmp1 = ((n2-kappa2 - 1)*(n2-kappa2+2) + four_n2_kappa2)/denom;
+ const double Fm = tmp1*tmp1 + Em*Em;
+
+ /* Gyration radius */
+ const double Rg = 0.5 * Dp * pow(Np/kf, 1.0/Df); /* [nm] */
+ const double Rg2 = Rg*Rg;
+ const double g = pow(1 + 4*k2*Rg2/(3*Df), -Df*0.5);
+
+ /* Cross sections, [nm^3/aggrefate] */
+ const double Cabs = Np * (4*PI*xp3)/(k*k) * Em;
+ const double Csca = Np*Np * (8*PI*xp3*xp3) / (3*k2) * Fm * g;
+
+ /* Radiative coefficients */
+ const double Va = (Np*PI*Dp*Dp*Dp) / 6; /* [nm^3] */
+ const double rho = Fv / Va; /* [#aggregates / nm^3] */
+ const double tmp2 = 1.e-9 * rho;
+
+ ASSERT(denom != 0 && Df != 0 && Dp != 0);
+ radcoefs->ka = tmp2 * Cabs;
+ radcoefs->ks = tmp2 * Csca;
+ radcoefs->kext = radcoefs->ka + radcoefs->ks;
+ }
+}
+
+#endif /* ATRSTM_RADCOEFS_H */
diff --git a/src/atrstm_setup_octrees.c b/src/atrstm_setup_octrees.c
@@ -18,7 +18,7 @@
#include "atrstm_c.h"
#include "atrstm_cache.h"
#include "atrstm_log.h"
-#include "atrstm_optprops.h"
+#include "atrstm_radcoefs.h"
#include "atrstm_partition.h"
#include "atrstm_setup_octrees.h"
#include "atrstm_svx.h"
@@ -86,24 +86,24 @@ register_primitive
}
static void
-voxel_commit_optprops_range
+voxel_commit_radcoefs_range
(float* vox,
const enum atrstm_component cpnt,
- const struct optprops* optprops_min,
- const struct optprops* optprops_max)
+ const struct radcoefs* radcoefs_min,
+ const struct radcoefs* radcoefs_max)
{
double ka[2], ks[2], kext[2];
float vox_ka[2], vox_ks[2], vox_kext[2];
ASSERT(vox);
- ASSERT(optprops_min);
- ASSERT(optprops_max);
+ ASSERT(radcoefs_min);
+ ASSERT(radcoefs_max);
- ka[0] = optprops_min->ka;
- ka[1] = optprops_max->ka;
- ks[0] = optprops_min->ks;
- ks[1] = optprops_max->ks;
- kext[0] = optprops_min->kext;
- kext[1] = optprops_max->kext;
+ ka[0] = radcoefs_min->ka;
+ ka[1] = radcoefs_max->ka;
+ ks[0] = radcoefs_min->ks;
+ ks[1] = radcoefs_max->ks;
+ kext[0] = radcoefs_min->kext;
+ kext[1] = radcoefs_max->kext;
/* Ensure that the single precision bounds include their double precision
* version. */
@@ -115,12 +115,12 @@ voxel_commit_optprops_range
if(kext[1] != (float)kext[1]) kext[1] = nextafterf((float)kext[1], FLT_MAX);
/* Fetch the range of the voxel optical properties */
- vox_ka[0] = vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN);
- vox_ka[1] = vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX);
- vox_ks[0] = vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN);
- vox_ks[1] = vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX);
- vox_kext[0] = vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN);
- vox_kext[1] = vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX);
+ vox_ka[0] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_MIN);
+ vox_ka[1] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_MAX);
+ vox_ks[0] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_MIN);
+ vox_ks[1] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_MAX);
+ vox_kext[0] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MIN);
+ vox_kext[1] = vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MAX);
/* Update the range of the voxel optical properties */
vox_ka[0] = MMIN(vox_ka[0], (float)ka[0]);
@@ -131,12 +131,12 @@ voxel_commit_optprops_range
vox_kext[1] = MMAX(vox_kext[1], (float)kext[1]);
/* Commit the upd to the voxel */
- vox_set(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN, vox_ka[0]);
- vox_set(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX, vox_ka[1]);
- vox_set(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN, vox_ks[0]);
- vox_set(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX, vox_ks[1]);
- vox_set(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN, vox_kext[0]);
- vox_set(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX, vox_kext[1]);
+ vox_set(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_MIN, vox_ka[0]);
+ vox_set(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_MAX, vox_ka[1]);
+ vox_set(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_MIN, vox_ks[0]);
+ vox_set(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_MAX, vox_ks[1]);
+ vox_set(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MIN, vox_kext[0]);
+ vox_set(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MAX, vox_kext[1]);
}
static res_T
@@ -223,21 +223,21 @@ voxelize_partition
intersect = suvm_polyhedron_intersect_aabb(&poly, vxl_low, vxl_upp);
if(intersect != SUVM_INTERSECT_NONE) {
- struct optprops optprops_min;
- struct optprops optprops_max;
+ struct radcoefs radcoefs_min;
+ struct radcoefs radcoefs_max;
float* vox = part_get_voxel(part, mcode[2]);
/* Currently, gas is not handled */
- optprops_min = OPTPROPS_NULL;
- optprops_max = OPTPROPS_NULL;
- voxel_commit_optprops_range
- (vox, ATRSTM_CPNT_GAS, &optprops_min, &optprops_max);
+ radcoefs_min = RADCOEFS_NULL;
+ radcoefs_max = RADCOEFS_NULL;
+ voxel_commit_radcoefs_range
+ (vox, ATRSTM_CPNT_GAS, &radcoefs_min, &radcoefs_max);
- res = primitive_compute_optprops_range
- (atrstm, refract_id, &prim, &optprops_min, &optprops_max);
+ res = primitive_compute_radcoefs_range
+ (atrstm, refract_id, &prim, &radcoefs_min, &radcoefs_max);
if(UNLIKELY(res != RES_OK)) goto error;
- voxel_commit_optprops_range
- (vox, ATRSTM_CPNT_SOOT, &optprops_min, &optprops_max);
+ voxel_commit_radcoefs_range
+ (vox, ATRSTM_CPNT_SOOT, &radcoefs_min, &radcoefs_max);
}
}
}
@@ -429,19 +429,19 @@ voxels_merge_component
FOR_EACH(ivox, 0, nvoxs) {
const float* vox = voxels[ivox];
- ka[0] = MMIN(ka[0], vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN));
- ka[1] = MMAX(ka[1], vox_get(vox, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX));
- ks[0] = MMIN(ks[0], vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN));
- ks[1] = MMAX(ks[1], vox_get(vox, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX));
- kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN));
- kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX));
+ ka[0] = MMIN(ka[0], vox_get(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_MIN));
+ ka[1] = MMAX(ka[1], vox_get(vox, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_MAX));
+ ks[0] = MMIN(ks[0], vox_get(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_MIN));
+ ks[1] = MMAX(ks[1], vox_get(vox, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_MAX));
+ kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MIN));
+ kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MAX));
}
- vox_set(dst, cpnt, ATRSTM_Ka, ATRSTM_SVX_MIN, ka[0]);
- vox_set(dst, cpnt, ATRSTM_Ka, ATRSTM_SVX_MAX, ka[1]);
- vox_set(dst, cpnt, ATRSTM_Ks, ATRSTM_SVX_MIN, ks[0]);
- vox_set(dst, cpnt, ATRSTM_Ks, ATRSTM_SVX_MAX, ks[1]);
- vox_set(dst, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN, kext[0]);
- vox_set(dst, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX, kext[1]);
+ vox_set(dst, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_MIN, ka[0]);
+ vox_set(dst, cpnt, ATRSTM_RADCOEF_Ka, ATRSTM_SVX_MAX, ka[1]);
+ vox_set(dst, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_MIN, ks[0]);
+ vox_set(dst, cpnt, ATRSTM_RADCOEF_Ks, ATRSTM_SVX_MAX, ks[1]);
+ vox_set(dst, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MIN, kext[0]);
+ vox_set(dst, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MAX, kext[1]);
}
static void
@@ -472,8 +472,8 @@ voxels_component_compute_kext_range
/* Compute the Kext range of the submitted voxels */
FOR_EACH(ivox, 0, nvoxs) {
const float* vox = voxels[ivox].data;
- kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MIN));
- kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_Kext, ATRSTM_SVX_MAX));
+ kext[0] = MMIN(kext[0], vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MIN));
+ kext[1] = MMAX(kext[1], vox_get(vox, cpnt, ATRSTM_RADCOEF_Kext, ATRSTM_SVX_MAX));
}
}
diff --git a/src/atrstm_svx.h b/src/atrstm_svx.h
@@ -25,21 +25,21 @@
* linearly as N single precision floating point data, with N computed as
* bellow:
*
- * N = ATRSTM_PROPS_COUNT__ #optical properties per voxel
+ * N = ATRSTM_RADCOEFS_COUNT__ #optical properties per voxel
* * ATRSTM_SVX_OPS_COUNT__ #supported operations on each properties
* * ATRSTM_CPNTS_COUNT__; #components on which properties are defined
*
* In a given voxel, the index `id' in [0, N-1] corresponding to the optical
- * property `enum atrstm_property P' of the component `enum atrstm_component
+ * property `enum atrstm_radcoef P' of the component `enum atrstm_component
* C' according to the operation `enum atrstm_svx_op O' is
* then computed as bellow:
*
* id = C * NFLOATS_PER_CPNT + P * HTSKY_SVX_OPS_COUNT__ + O;
- * NFLOATS_PER_CPNT = ATRSTM_SVX_OPS_COUNT__ * ATRSTM_PROPS_COUNT__;
+ * NFLOATS_PER_CPNT = ATRSTM_SVX_OPS_COUNT__ * ATRSTM_RADCOEFS_COUNT__;
*/
/* Constant defining the number of floating point data per component */
-#define NFLOATS_PER_CPNT (ATRSTM_SVX_OPS_COUNT__ * ATRSTM_PROPS_COUNT__)
+#define NFLOATS_PER_CPNT (ATRSTM_SVX_OPS_COUNT__ * ATRSTM_RADCOEFS_COUNT__)
/* Constant defining the overall number of floating point data of a voxel */
#define NFLOATS_PER_VOXEL (NFLOATS_PER_CPNT * ATRSTM_CPNTS_COUNT__)
@@ -48,12 +48,12 @@ static FINLINE float
vox_get
(const float* vox,
const enum atrstm_component cpnt,
- const enum atrstm_property prop,
+ const enum atrstm_radcoef prop,
const enum atrstm_svx_op op)
{
ASSERT(vox);
ASSERT((unsigned)cpnt < ATRSTM_CPNTS_COUNT__);
- ASSERT((unsigned)prop < ATRSTM_PROPS_COUNT__);
+ ASSERT((unsigned)prop < ATRSTM_RADCOEFS_COUNT__);
ASSERT((unsigned)op < ATRSTM_SVX_OPS_COUNT__);
return vox[cpnt*NFLOATS_PER_CPNT + prop*ATRSTM_SVX_OPS_COUNT__ + op];
}
@@ -62,13 +62,13 @@ static FINLINE void
vox_set
(float* vox,
const enum atrstm_component cpnt,
- const enum atrstm_property prop,
+ const enum atrstm_radcoef prop,
const enum atrstm_svx_op op,
const float val)
{
ASSERT(vox);
ASSERT((unsigned)cpnt < ATRSTM_CPNTS_COUNT__);
- ASSERT((unsigned)prop < ATRSTM_PROPS_COUNT__);
+ ASSERT((unsigned)prop < ATRSTM_RADCOEFS_COUNT__);
ASSERT((unsigned)op < ATRSTM_SVX_OPS_COUNT__);
vox[cpnt*NFLOATS_PER_CPNT + prop*ATRSTM_SVX_OPS_COUNT__ + op] = val;
}
@@ -81,7 +81,7 @@ vox_clear(float* vox)
FOR_EACH(cpnt, 0, ATRSTM_CPNTS_COUNT__) {
int prop;
- FOR_EACH(prop, 0, ATRSTM_PROPS_COUNT__) {
+ FOR_EACH(prop, 0, ATRSTM_RADCOEFS_COUNT__) {
vox_set(vox, cpnt, prop, ATRSTM_SVX_MIN, FLT_MAX);
vox_set(vox, cpnt, prop, ATRSTM_SVX_MAX,-FLT_MAX);
}
diff --git a/src/test_atrstm_optprops.c b/src/test_atrstm_optprops.c
@@ -1,43 +0,0 @@
-/* Copyright (C) 2020 CNRS
- *
- * 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/>. */
-
-#include "atrstm_optprops.h"
-
-int
-main(int argc, char** argv)
-{
- const double ka_ref = 5.7382401729092799E-019;
- const double ks_ref = 7.2169062018378995E-024;
-
- struct optprops optprops;
- struct optprops_compute_args args = OPTPROPS_COMPUTE_ARGS_NULL;
- (void)argc, (void)argv;
-
- args.lambda = 633;
- args.n = 1.90;
- args.kappa = 0.55;
- args.gyration_radius_prefactor = 1.70;
- args.fractal_dimension = 1.75;
- args.soot_volumic_fraction = 1.e-7;
- args.soot_primary_particles_count = 100;
- args.soot_primary_particles_diameter = 1;
-
- optprops_compute(&optprops, &args);
-
- printf("ka = %g; ks = %g\n", optprops.ka, optprops.ks);
- CHK(eq_eps(optprops.ka, ka_ref, ka_ref*1.e-6));
- CHK(eq_eps(optprops.ks, ks_ref, ks_ref*1.e-6));
- return 0;
-}
diff --git a/src/test_atrstm_radcoefs.c b/src/test_atrstm_radcoefs.c
@@ -0,0 +1,43 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * 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/>. */
+
+#include "atrstm_radcoefs.h"
+
+int
+main(int argc, char** argv)
+{
+ const double ka_ref = 5.7382401729092799E-019;
+ const double ks_ref = 7.2169062018378995E-024;
+
+ struct radcoefs radcoefs;
+ struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL;
+ (void)argc, (void)argv;
+
+ args.lambda = 633;
+ args.n = 1.90;
+ args.kappa = 0.55;
+ args.gyration_radius_prefactor = 1.70;
+ args.fractal_dimension = 1.75;
+ args.soot_volumic_fraction = 1.e-7;
+ args.soot_primary_particles_count = 100;
+ args.soot_primary_particles_diameter = 1;
+
+ radcoefs_compute(&radcoefs, &args);
+
+ printf("ka = %g; ks = %g\n", radcoefs.ka, radcoefs.ks);
+ CHK(eq_eps(radcoefs.ka, ka_ref, ka_ref*1.e-6));
+ CHK(eq_eps(radcoefs.ks, ks_ref, ks_ref*1.e-6));
+ return 0;
+}
