atrstm

commit e86e8da52a10ecb4239e884e0d90edd1a5c074a2
parent ae8c427fdaac180878646e87017394aeb908195a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 15 Jan 2021 13:40:41 +0100

Fix optprops computation

Handle cases where #primary particles or volumic fraction is null.

Diffstat:
M
src/atrstm_optprops.h
 | 
76
++++++++++++++++++++++++++++++++++++++++++----------------------------------

1 file changed, 42 insertions(+), 34 deletions(-)
diff --git a/src/atrstm_optprops.h b/src/atrstm_optprops.h
@@ -111,40 +111,48 @@ optprops_compute
   const double kf = args->gyration_radius_prefactor;
   const double Df = args->fractal_dimension;
 
-  /* 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 = (PI * Dp*Dp*Dp) / 6; /* [nm^3] */
-  const double rho = Fv / Va; /* [#aggregates / nm^3] */
-  const double tmp2 = 1.e-9 * rho;
-  optprops->ka = tmp2 * Cabs;
-  optprops->ks = tmp2 * Csca;
-  optprops->kext = optprops->ka + optprops->ks;
+  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 = (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 */