star-line

commit fccb72868c926dec8b6458e80c328fb8dad9c3c5
parent 2ee47d0477c22658a4130bf1dcd72107de7a7e6f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 17 Mar 2026 09:06:26 +0100

sln-slab: compute the transmittance

Implement the realisation function that computes the (normal)
transmittance in the slab using the Monte Carlo algorithm by Yaniss
Nyffenegger-Péré et al. [1]. It is estimated by sampling the transition
lines, using the provided acceleration structure to sample them based on
their importance.

Add several debug and runtime checks to ensure that the assertions
regarding the acceleration structure are properly verified. And the very
first tests reveal issues in how the tree is constructed, such that it
is no longer guaranteed that a line is sampled with a valid
probability - that is, less than 1. And that is precisely the main
objective of this tool: to highlight such problems, in this case,
certainly a problem with how the lines are meshed.

[1] Y. Nyffenegger-Péré et al., “Spectrally refined unbiased monte
    carlo estimate of the earth's global radiative cooling”, Proceedings
    of the National Academy of Sciences, 5, 121, e2315492121, 2024.

Diffstat:
M
src/sln_slab.c
 | 
124
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----

1 file changed, 119 insertions(+), 5 deletions(-)
diff --git a/src/sln_slab.c b/src/sln_slab.c
@@ -297,16 +297,130 @@ error:
   goto exit;
 }
 
+static const struct sln_node*
+sample_line
+  (const struct cmd* cmd,
+   struct ssp_rng* rng,
+   const double nu/*[cm^-1]*/,
+   double* out_proba)
+{
+  const struct sln_node* node = NULL;
+  double proba = 1; /* Proba to sample the line */
+  size_t depth = 0;
+  ASSERT(cmd && rng);
+
+  node = sln_tree_get_root(cmd->tree);
+
+  for(depth=0; ;++depth) {
+    const struct sln_node* node0 = NULL;
+    const struct sln_node* node1 = NULL;
+    struct sln_mesh mesh0 = SLN_MESH_NULL;
+    struct sln_mesh mesh1 = SLN_MESH_NULL;
+    double proba0 = 0;
+    double ka_max0 = 0;
+    double ka_max1 = 0;
+    double r = 0;
+
+    if(sln_node_is_leaf(node)) break;
+
+    node0 = sln_node_get_child(node, 0);
+    node1 = sln_node_get_child(node, 1);
+    SLN(node_get_mesh(cmd->tree, node0, &mesh0));
+    SLN(node_get_mesh(cmd->tree, node1, &mesh1));
+
+    ka_max0 = sln_mesh_eval(&mesh0, nu);
+    ka_max1 = sln_mesh_eval(&mesh1, nu);
+    /* TODO handle ka_max[0,1] == 0 */
+#ifndef NDEBUG
+    {
+      struct sln_mesh mesh = SLN_MESH_NULL;
+      double ka_max = 0;
+      SLN(node_get_mesh(cmd->tree, node, &mesh));
+      ka_max = sln_mesh_eval(&mesh, nu);
+      ASSERT(ka_max >= ka_max0 + ka_max1);
+    }
+#endif
+
+    /* Compute the probability to sample the 1st node */
+    proba0 = ka_max0 / (ka_max0 + ka_max1);
+
+    /* Sample the child node */
+    r = ssp_rng_canonical(rng);
+    if(r < proba0) {
+      node = node0;
+      proba *= proba0;
+    } else {
+      node = node1;
+      proba *= (1-proba0);
+    }
+  }
+
+  *out_proba = proba;
+  return node;
+}
+
 static double
 realisation(const struct cmd* cmd, struct ssp_rng* rng)
 {
+  const struct sln_node* node = NULL;
+  struct sln_mesh mesh = SLN_MESH_NULL;
+  double ka_max = 0;
+  double dst_remain = 0;
+  double dst = 0; /* Length */
+  double nu = 0;
+  size_t ncollisions = 0; /* For debug */
   ASSERT(cmd && rng);
-  (void)cmd; /* Avoid unused variable warning */
 
-  if(ssp_rng_canonical(rng) < 0.5) {
-    return 2.0;
-  } else {
-    return 8.0;
+  dst_remain = cmd->args.thickness;
+
+  /* Uniformly sample the spectral dimension */
+  nu = ssp_rng_uniform_double
+    (rng, cmd->args.spectral_range[0], cmd->args.spectral_range[1]);
+
+  /* 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));
+  ka_max = sln_mesh_eval(&mesh, nu);
+
+  for(ncollisions=0; ; ++ncollisions) {
+    /* Sample a traversal distance */
+    dst = ssp_ran_exp(rng, ka_max);
+
+    if(dst > dst_remain) {
+      return cmd->args.spectral_range[1] - cmd->args.spectral_range[0];
+    } else {
+      double proba_abs = 0;
+      double line_proba = 0;
+      double line_ha = 0; /* Value of the line */
+      double r = 0; /* Random number */
+
+      /* Importance sampling of a line */
+      node = sample_line(cmd, rng, nu, &line_proba);
+
+      /* Check whether the sampled node corresponds to a sampled line */
+      #ifndef NDEBUG
+      {
+        struct sln_node_desc desc = SLN_NODE_DESC_NULL;
+        SLN(node_get_desc(node, &desc));
+        ASSERT(desc.nlines == 1);
+      }
+      #endif
+
+      /* Get the line value */
+      line_ha = sln_node_eval(cmd->tree, node, nu);
+
+      proba_abs = line_ha / (line_proba*ka_max);
+      if(proba_abs > 1) {
+        VFATAL("Invalid probability %g\n", ARG1(proba_abs));
+      }
+
+      r = ssp_rng_canonical(rng);
+      if(r < proba_abs) {
+        return 0;
+      } else { /* Null collision */
+        dst_remain -= dst;
+      }
+    }
   }
 }