htrdr

commit bec0fdc8f1f2f180f4ef6dc51bbdd494818ab8d8
parent 62babff860e9dd45960eeb023ec255d8e10ec3d9
Author: Arfaux Anthony <anthony.arfaux@obspm.fr>
Date:   Mon, 19 May 2025 14:25:12 +0200

fix LW NaN issue

Modify the behavior in lw radiative budget mode to avoid the apparition
of infinity and NaN related to the absence of direct contribution.

Diffstat:
M
src/planets/htrdr_planets_solve_volrad_budget.c
 | 
20
+++++++++++++-------

1 file changed, 13 insertions(+), 7 deletions(-)
diff --git a/src/planets/htrdr_planets_solve_volrad_budget.c b/src/planets/htrdr_planets_solve_volrad_budget.c
@@ -230,6 +230,11 @@ realisation
   position_sampling(args, volrad_mesh_desc, pos);
   ssp_ran_sphere_uniform(args->rng, dir, &dir_pdf);
 
+  S = get_source(cmd, pos, wlen); /* [W/m^2/sr/m] */
+
+  ka = get_ka(cmd, pos, iband, iquad);
+  wlen_pdf_m = wlen_pdf_nm * 1.e9; /* Transform pdf from nm^-1 to m^-1 */
+
   /* Compute the radiance in W/m^2/sr/m */
   d3_set(rad_args.path_org, pos);
   rad_args.rng = args->rng;
@@ -243,9 +248,12 @@ realisation
      * path for the sampled direction and position: the radiance is considered
      * as purely diffuse. */
     d3_set(rad_args.path_dir, dir);
-    L_direct  = 0.0;
     L_diffuse = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
 
+    /* Calculate the weights [W/m^3] */
+    weights[DIRECT]  = 0.0;
+    weights[DIFFUSE] = ka * (L_diffuse - S) / (wlen_pdf_m * dir_pdf);
+
   } else {
     /* In the so-called shortwave region (actually, the radiation due the
      * external source) is decomposed in its direct and diffuse components */
@@ -260,16 +268,14 @@ realisation
     d3_set(rad_args.path_dir, dir);
     rad_args.component = PLANETS_RADIANCE_CPNT_DIFFUSE;
     L_diffuse = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */
-  }
 
-  S = get_source(cmd, pos, wlen); /* [W/m^2/sr/m] */
 
-  ka = get_ka(cmd, pos, iband, iquad);
-  wlen_pdf_m = wlen_pdf_nm * 1.e9; /* Transform pdf from nm^-1 to m^-1 */
+    /* Calculate the weights [W/m^3] */
+    weights[DIRECT]  = ka * (L_direct  - S) / (wlen_pdf_m * dir_src_pdf);
+    weights[DIFFUSE] = ka * (L_diffuse - S) / (wlen_pdf_m * dir_pdf);
+  }
 
   /* Calculate the weights [W/m^3] */
-  weights[DIRECT]  = ka * (L_direct  - S) / (wlen_pdf_m * dir_src_pdf);
-  weights[DIFFUSE] = ka * (L_diffuse - S) / (wlen_pdf_m * dir_pdf);
   weights[TOTAL]   = weights[DIRECT] + weights[DIFFUSE];
 }