commit 576f8d8e1fa7cf57c2c2af93b70e62459b69c7ef
parent 53207b075ac1f7b5c5db5210048a7376cedef533
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 23 Mar 2026 16:58:04 +0100
sln-slab: add the estimate of emissivity
The emissivity is computed in the same realisation function, using the
same path space. The realisation function simply outputs 2 weights
instead of 1.
This new computation requires the evaluation of the monochromatic
Planck's function. The Star-Blackbody library, which provides this
function, is therefore added as a dependency of the sln-slab utility,
and therefore of the Star-Line project. However, the monochromatic
Planck's function is very simple and could therefore be copied/pasted
into the sln-slab source code, if adding a new dependency was a problem.
But since this dependency does not seem to be a problem, it is better to
leave Planck's function in one place.
Diffstat:
5 files changed, 103 insertions(+), 40 deletions(-)
diff --git a/Makefile b/Makefile
@@ -69,6 +69,7 @@ libsln.o: $(OBJ)
.config: config.mk
$(PKG_CONFIG) --atleast-version $(LBLU_VERSION) lblu
$(PKG_CONFIG) --atleast-version $(RSYS_VERSION) rsys
+ $(PKG_CONFIG) --atleast-version $(SBB_VERSION) sbb
$(PKG_CONFIG) --atleast-version $(SHTR_VERSION) shtr
$(PKG_CONFIG) --atleast-version $(SSP_VERSION) star-sp
echo "config done" > $@
@@ -94,8 +95,8 @@ LIBS_UTIL = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys sln-local shtr)
CFLAGS_UTIL = -std=c89 $(CFLAGS_EXE) $(INCS_UTIL)
LDFLAGS_UTIL = $(LDFLAGS_EXE) $(LIBS_UTIL)
-INCS_SLAB = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys sln-local shtr star-sp)
-LIBS_SLAB = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys sln-local shtr star-sp) -lm
+INCS_SLAB = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --cflags rsys sln-local sbb shtr star-sp)
+LIBS_SLAB = $$($(PKG_CONFIG_LOCAL) $(PCFLAGS) --libs rsys sln-local sbb shtr star-sp) -lm
CFLAGS_SLAB = -std=c89 $(CFLAGS_EXE) $(INCS_SLAB) -fopenmp
LDFLAGS_SLAB = $(LDFLAGS_EXE) $(LIBS_SLAB)
diff --git a/README.md b/README.md
@@ -10,6 +10,7 @@ Build an acceleration structure to speed up line importance sampling.
- Line By Line Utils
- RSys
- Star HITRAN
+- Star BlackBody
## Installation
diff --git a/config.mk b/config.mk
@@ -31,7 +31,10 @@ PCFLAGS = $(PCFLAGS_$(LIB_TYPE))
LBLU_VERSION = 0.0.1
RSYS_VERSION = 0.14
SHTR_VERSION = 0
-SSP_VERSION = 0.15 # Dependency for the sln-slab util
+
+# Dependencies for the sln-slab util
+SSP_VERSION = 0.15
+SBB_VERSION = 0.0
INCS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags lblu rsys shtr)
LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs lblu rsys shtr) -lm
diff --git a/doc/sln-slab.1 b/doc/sln-slab.1
@@ -15,7 +15,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/>.
-.Dd March 19, 2026
+.Dd March 23, 2026
.Dt SLN-SLAB 1
.Os
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@@ -36,10 +36,10 @@
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
.Sh DESCRIPTION
.Nm
-calculates the transmittance in a one-dimensional homogeneous slab of
-arbitrary thickness using a Monte Carlo algorithm that samples the
-spectral lines that make up the gas mixture.
-This computation is accelerated by sampling the lines based on the
+calculates the transmittance and the emissivity in a one-dimensional
+homogeneous slab of arbitrary thickness using a Monte Carlo algorithm
+that samples the spectral lines that make up the gas mixture.
+These computations are accelerated by sampling the lines based on the
magnitude of their contribution to the mixture’s spectrum, so that few
Monte Carlo runs are required to estimate the transmissivity with a high
degree of confidence.
@@ -54,20 +54,21 @@ More than just a numerical simulation tool,
.Nm
is primarily designed to validate the aforementioned acceleration
structure in relation to its intended use, namely radiative transfer
-calculations.
+computations.
Thus, not only could an error be returned in the event of a problem with
the structure or its use, but the computed value can also contribute to
-this validation through its comparison with the result of a calculation
+this validation through its comparison with the result of a computation
of the same quantity performed by another radiative transfer code.
.Pp
The output of
.Nm
-displays the estimated transmittance, its standard deviation, and the
-number of Monte Carlo realisations rejected due to issues encountered
-during the computation, such as numerical uncertainty.
-The displayed values are presented as follows:
+displays the estimated transmittance and emissivity, their standard
+deviation, and the number of Monte Carlo realisations rejected due to
+issues encountered during the computation, such as numerical
+uncertainty.
+Each estimate is displayed on a line formatted as follows:
.Bd -literal -offset Ds
-"%e %e %lu\en", transmittance, std_err, rejects_count
+"%-16s: %e +/- %e; %lu\en", name, estimate, std_err, rejects_count
.Ed
.Pp
The options are as follows:
@@ -132,8 +133,8 @@ The maximum is 3.
.Ex -std
.\""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
.Sh EXAMPLES
-Estimate the transmittance between 100 and 2500 cm^-1 for a slab 2
-meters thick.
+Estimate the transmittance and emissivity between 100 and 2500 cm^-1 for
+a slab 2 meters thick.
The slab consists of a homogeneous gas mixture of H2O, CO2 and CO
molecules.
The thermodynamic properties of the mixture, such as its pressure,
diff --git a/src/sln_slab.c b/src/sln_slab.c
@@ -21,6 +21,7 @@
#include "sln.h"
#include <star/shtr.h>
+#include <star/sbb.h>
#include <star/ssp.h>
#include <rsys/cstr.h>
@@ -30,6 +31,14 @@
#include <unistd.h> /* getopt */
+enum estimate {
+ TRANSMISSIVITY,
+ EMISSIVITY,
+ ESTIMATE_COUNT__
+};
+
+#define WAVENUMBER_TO_WAVELENGTH(Nu/* [cm^-1] */) (1.e-2/(Nu))/*[m]*/
+
struct args {
const char* tree; /* Acceleration structure */
const char* molparams;
@@ -65,7 +74,6 @@ struct accum {
size_t count;
};
#define ACCUM_NULL__ {0}
-static const struct accum ACCUM_NULL = ACCUM_NULL__;
/*******************************************************************************
* Helper functions
@@ -405,10 +413,38 @@ check_proba(const struct cmd* cmd, const double proba)
return RES_BAD_ARG;
}
+static FINLINE double /* [W/m^2/sr/cm^-1] */
+planck
+ (const double nu/* [cm^-1] */,
+ const double T/* [K] */)
+{
+ const double lambda = WAVENUMBER_TO_WAVELENGTH(nu); /* [m] */
+ const double planck1 = sbb_planck_monochromatic(lambda, T); /* [W/m^2/sr/m] */
+ const double planck2 = planck1 * (1.0e-2/(nu*nu)); /* [W/m^2/sr/cm^-1] */
+ ASSERT(T >= 0);
+ return planck2;
+}
+
+static INLINE const char*
+estimate_cstr(const enum estimate estimate)
+{
+ const char* cstr = NULL;
+ switch(estimate) {
+ case TRANSMISSIVITY: cstr="transmissivity"; break;
+ case EMISSIVITY: cstr="emissivity"; break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ return cstr;
+}
+
static res_T
-realisation(const struct cmd* cmd, struct ssp_rng* rng, double* out_weight)
+realisation
+ (const struct cmd* cmd,
+ struct ssp_rng* rng,
+ double out_weights[ESTIMATE_COUNT__])
{
/* Acceleration structure */
+ struct sln_tree_desc tree_desc = SLN_TREE_DESC_NULL;
const struct sln_node* node = NULL;
struct sln_mesh mesh = SLN_MESH_NULL;
@@ -418,11 +454,16 @@ realisation(const struct cmd* cmd, struct ssp_rng* rng, double* out_weight)
size_t ncollisions = 0; /* Number of null collisions */
/* Miscellaneous */
- double nu = 0; /* Sampled wavenumber */
- double w = 0; /* Monte Carlo weight */
+ double w[ESTIMATE_COUNT__] = {0, 0}; /* Monte Carlo weight */
+ double nu = 0; /* Sampled wavenumber [cm^-1] */
+ double nu_range = 0; /* Spectral range [cm^-1] */
+ int i = 0;
res_T res = RES_OK;
- ASSERT(cmd && rng && out_weight); /* Check pre-conditions */
+ ASSERT(cmd && rng && out_weights); /* Check pre-conditions */
+
+ /* Precompute the spectral range */
+ nu_range = cmd->args.spectral_range[1] - cmd->args.spectral_range[0];
/* Initialize the total distance to traverse with the thickness of the slab */
dst_remain = cmd->args.thickness;
@@ -431,6 +472,8 @@ realisation(const struct cmd* cmd, struct ssp_rng* rng, double* out_weight)
nu = ssp_rng_uniform_double
(rng, cmd->args.spectral_range[0], cmd->args.spectral_range[1]);
+ SLN(tree_get_desc(cmd->tree, &tree_desc));
+
/* Retrieve the ka_max of the spectrum at the sampled nu */
node = sln_tree_get_root(cmd->tree);
SLN(node_get_mesh(cmd->tree, node, &mesh));
@@ -445,7 +488,11 @@ realisation(const struct cmd* cmd, struct ssp_rng* rng, double* out_weight)
dst = ssp_ran_exp(rng, ka_max); /* Sample a traversal distance */
- if(dst > dst_remain) { w=1.0; break; } /* No absorption in the slab */
+ if(dst > dst_remain) { /* No absorption in the slab */
+ w[TRANSMISSIVITY] = 1.0;
+ w[EMISSIVITY] = 0.0;
+ break;
+ }
/* Importance sampling of a line */
node = sample_line(cmd, rng, nu, &line_proba);
@@ -458,16 +505,20 @@ realisation(const struct cmd* cmd, struct ssp_rng* rng, double* out_weight)
if((res = check_proba(cmd, proba_abs)) != RES_OK) goto error;
r = ssp_rng_canonical(rng);
- if(r < proba_abs) { w=0; break; } /* An absorption occurs */
+ if(r < proba_abs) { /* An absorption occurs */
+ w[TRANSMISSIVITY] = 0.0;
+ w[EMISSIVITY] = planck(nu, tree_desc.temperature)*nu_range; /*[W/m^2/sr]*/
+ break;
+ }
dst_remain -= dst; /* This was a null transition. Go on */
}
exit:
- *out_weight = w;
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) out_weights[i] = w[i];
return res;
error:
- w = NaN;
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) w[i] = NaN;
goto exit;
}
@@ -478,10 +529,7 @@ cmd_run(const struct cmd* cmd)
struct ssp_rng** rngs = NULL;
/* Monte Carlo */
- struct accum accum = ACCUM_NULL;
- double E = 0; /* Expected value */
- double V = 0; /* Variance */
- double SE = 0; /* Standard Error */
+ struct accum accum[ESTIMATE_COUNT__] = {0};
int64_t i = 0; /* Index of the realisation */
size_t nrejects = 0; /* Number of rejected realisations */
@@ -506,20 +554,23 @@ cmd_run(const struct cmd* cmd)
#pragma omp parallel for schedule(static)
for(i = 0; i < (int64_t)nrealisations; ++i) {
+ double w[ESTIMATE_COUNT__] = {0}; /* Monte Carlo weights */
const int ithread = omp_get_thread_num();
int pcent = 0;
- double w = 0; /* Monte Carlo weight */
res_T res_realisation = RES_OK;
- res_realisation = realisation(cmd, rngs[ithread], &w);
+ res_realisation = realisation(cmd, rngs[ithread], w);
#pragma omp critical
{
/* Update the Monte Carlo accumulator */
if(res_realisation == RES_OK) {
- accum.sum += w;
- accum.sum2 += w*w;
- accum.count += 1;
+ int iestim = 0;
+ FOR_EACH(iestim, 0, ESTIMATE_COUNT__) {
+ accum[iestim].sum += w[iestim];
+ accum[iestim].sum2 += w[iestim]*w[iestim];
+ accum[iestim].count += 1;
+ }
}
if(cmd->args.verbose >= 3) {
@@ -536,12 +587,18 @@ cmd_run(const struct cmd* cmd)
#undef PROGRESS_MSG
- E = accum.sum / (double)accum.count;
- V = accum.sum2 / (double)accum.count - E*E;
- SE = sqrt(V/(double)accum.count);
- nrejects = nrealisations - accum.count;
+ nrejects = nrealisations - accum[0].count;
+
+ FOR_EACH(i, 0, ESTIMATE_COUNT__) {
+ const double E = accum[i].sum / (double)accum[i].count;
+ const double V = accum[i].sum2 / (double)accum[i].count - E*E;
+ const double SE = sqrt(V/(double)accum[i].count);
- printf("%e %e %lu\n", E, SE, (unsigned long)nrejects);
+ /* Assume that the number of realisations is the same for all estimates */
+ ASSERT(accum[i].count == accum[0].count);
+
+ printf("%-16s: %e +/- %e; %lu\n", estimate_cstr(i), E, SE, (unsigned long)nrejects);
+ }
exit:
delete_per_thread_rngs(cmd, rngs);
