atrstm

commit 3b694e0b921857e2149d129a8619d0cd19894064
parent cc23f79b77a6aec67d501013396c3ce941d6c975
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 21 Jan 2021 12:15:22 +0100

Fix and update the API of the radcoefs_compute function

Diffstat:
M
src/atrstm_radcoefs.h
 | 
68
+++++++++++++++++++++++++++++++++++++++++++-------------------------
M
src/test_atrstm_radcoefs.c
 | 
25
+++++++++++++++++++++++--

2 files changed, 66 insertions(+), 27 deletions(-)
diff --git a/src/atrstm_radcoefs.h b/src/atrstm_radcoefs.h
@@ -16,10 +16,12 @@
 #ifndef ATRSTM_RADCOEFS_H
 #define ATRSTM_RADCOEFS_H
 
+#include "atrstm.h"
+
 #include <rsys/math.h>
 #include <rsys/rsys.h>
 
-/* Radiative coefficients */
+/* Radiative coefficients in m^-1 */
 struct radcoefs {
   double ka; /* Absorption coefficient */
   double ks; /* Scattering coefficient */
@@ -37,6 +39,7 @@ struct radcoefs_compute_args {
   double soot_volumic_fraction; /* In m^3(soot) / m^3 */
   double soot_primary_particles_count;
   double soot_primary_particles_diameter; /* In nm */
+  int radcoefs_mask; /* Combination of atrstm_radcoef_flag */
 };
 #define RADCOEFS_COMPUTE_ARGS_NULL__ {0}
 static const struct radcoefs_compute_args RADCOEFS_COMPUTE_ARGS_NULL =
@@ -118,12 +121,16 @@ radcoefs_compute
   const double Dp = args->soot_primary_particles_diameter; /* [nm] */
   const double kf = args->gyration_radius_prefactor;
   const double Df = args->fractal_dimension;
+  const int ka = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) != 0;
+  const int ks = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) != 0;
+  const int kext = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) != 0;
+
+  /* Clean up radiative coefficients */
+  radcoefs->ka = 0;
+  radcoefs->ks = 0;
+  radcoefs->kext = 0;
 
-  if(Np == 0 || Fv == 0) {
-    radcoefs->ka = 0;
-    radcoefs->ks = 0;
-    radcoefs->kext = 0;
-  } else {
+  if(Np != 0 && Fv != 0 && (args->radcoefs_mask & ATRSTM_RADCOEFS_MASK_ALL)) {
     /* Precomputed values */
     const double n2 = n*n;
     const double kappa2 = kappa*kappa;
@@ -133,33 +140,44 @@ radcoefs_compute
     const double k = (2*PI) / lambda;
     const double k2 = k*k;
 
+    const double Va = (Np*PI*Dp*Dp*Dp) / 6; /* [nm^3] */
+    const double rho = Fv / Va; /* [#aggregates / nm^3] */
+    const double tmp0 = 1.e9 * rho;
+
     /* E(m), m = n + i*kappa */
-    const double tmp0 = (n2 - kappa2 + 2);
-    const double denom = tmp0*tmp0 + four_n2_kappa2;
+    const double tmp1 = (n2 - kappa2 + 2);
+    const double denom = tmp1*tmp1 + 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;
+    if(ka || kext) {
+      /* Absorption cross section, [nm^3/aggrefate] */
+      const double Cabs = Np * (4*PI*xp3)/(k2) * 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);
+      /* Absorption coefficient, [m^-1] */
+      radcoefs->ka = tmp0 * Cabs;
+    }
 
-    /* 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;
+    if(ks || kext) {
+      /* F(m), m = n + i*kappa */
+      const double tmp2 = ((n2-kappa2 - 1)*(n2-kappa2+2) + four_n2_kappa2)/denom;
+      const double Fm = tmp2*tmp2 + Em*Em;
 
-    /* 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;
+      /* 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);
+
+      /* Scattering cross section, [nm^3/aggrefate] */
+      const double Csca = Np*Np * (8*PI*xp3*xp3) / (3*k2) * Fm * g;
+
+      /* Scattering coefficient, [m^-1] */
+      radcoefs->ks = tmp0 * Csca;
+    }
 
+    if(kext) {
+      radcoefs->kext = radcoefs->ka + radcoefs->ks;
+    }
     ASSERT(denom != 0 && Df != 0 && Dp != 0);
-    radcoefs->ka = tmp2 * Cabs;
-    radcoefs->ks = tmp2 * Csca;
-    radcoefs->kext = radcoefs->ka + radcoefs->ks;
   }
 }
 
diff --git a/src/test_atrstm_radcoefs.c b/src/test_atrstm_radcoefs.c
@@ -18,8 +18,8 @@
 int
 main(int argc, char** argv)
 {
-  const double ka_ref = 5.7382401729092799E-019;
-  const double ks_ref = 7.2169062018378995E-024;
+  const double ka_ref = 5.7382401729092799E-1;
+  const double ks_ref = 7.2169062018378995E-6;
 
   struct radcoefs radcoefs;
   struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL;
@@ -33,11 +33,32 @@ main(int argc, char** argv)
   args.soot_volumic_fraction = 1.e-7;
   args.soot_primary_particles_count = 100;
   args.soot_primary_particles_diameter = 1;
+  args.radcoefs_mask = ATRSTM_RADCOEFS_MASK_ALL;
 
   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));
+  CHK(eq_eps(radcoefs.kext, ka_ref+ks_ref, (ka_ref+ks_ref)*1e-6));
+
+  args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ka;
+  radcoefs_compute(&radcoefs, &args);
+  CHK(eq_eps(radcoefs.ka, ka_ref, ka_ref*1.e-6));
+  CHK(radcoefs.ks == 0);
+  CHK(radcoefs.kext == 0);
+
+  args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks;
+  radcoefs_compute(&radcoefs, &args);
+  CHK(radcoefs.ka == 0);
+  CHK(eq_eps(radcoefs.ks, ks_ref, ks_ref*1.e-6));
+  CHK(radcoefs.kext == 0);
+
+  args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Kext;
+  /* Note that actually even though Ka and Ks are note required they are
+   * internally computed and thus are returned to the caller */
+  radcoefs_compute(&radcoefs, &args);
+  CHK(eq_eps(radcoefs.kext, ka_ref+ks_ref, (ka_ref+ks_ref)*1e-6));
+
   return 0;
 }