solstice-solver

commit 2b33ca2b14d7478141e4c097e2267b8563ce0216
parent 9734dc2f5e1308a1090c95d320999bd049888739
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue, 26 Sep 2017 10:41:44 +0200

Fix Pillbox distribution.

Pdf was wrong.

Diffstat:
M
src/ssol_ranst_sun_dir.c
 | 
23
+++++++++++++++--------
M
src/ssol_sun.c
 | 
2
+-

2 files changed, 16 insertions(+), 9 deletions(-)
diff --git a/src/ssol_ranst_sun_dir.c b/src/ssol_ranst_sun_dir.c
@@ -40,7 +40,7 @@ struct ran_buie_state {
 };
 
 struct ran_pillbox_state {
-  double radius;
+  double sin2_aperture;
   double basis[9];
 };
 
@@ -243,10 +243,17 @@ ran_pillbox_get
    double dir[3])
 {
   double pt[3];
+  double phi, sin2_theta, sin_theta, cos_theta;
 
-  ASSERT(ran && ran->state.pillbox.radius > 0);
-  ssp_ran_uniform_disk(rng, ran->state.pillbox.radius, pt);
-  pt[2] = 1;
+  ASSERT(ran && 0 <= ran->state.pillbox.sin2_aperture >= 0
+    && ran->state.pillbox.sin2_aperture <= 1);
+  sin2_theta = ssp_rng_uniform_double(rng, 0, ran->state.pillbox.sin2_aperture);
+  sin_theta = sqrt(sin2_theta);
+  cos_theta = sqrt(1 - sin2_theta);
+  phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
+  pt[0] = cos(phi) * sin_theta;
+  pt[1] = sin(phi) * sin_theta;
+  pt[2] = cos_theta;
   d33_muld3(dir, ran->state.pillbox.basis, pt);
   d3_normalize(dir, dir);
   return dir;
@@ -340,12 +347,12 @@ ranst_sun_dir_pillbox_setup
    const double aperture,
    const double dir[3])
 {
-  double radius;
-  if (!ran || !dir || aperture <= 0 || aperture >= PI )
+  double s;
+  if (!ran || !dir || aperture <= 0 || aperture > PI )
     return RES_BAD_ARG;
-  radius = tan(0.5 * aperture);
   ran->get = ran_pillbox_get;
-  ran->state.pillbox.radius = radius;
+  s = sin(aperture);
+    ran->state.pillbox.sin2_aperture = s * s;
   d33_basis(ran->state.pillbox.basis, dir);
   return RES_OK;
 }
diff --git a/src/ssol_sun.c b/src/ssol_sun.c
@@ -174,7 +174,7 @@ ssol_sun_set_pillbox_aperture(struct ssol_sun* sun, const double angle)
 {
   if(!sun
   || angle <= 0
-  || angle > PI * 0.5
+  || angle > PI
   || sun->type != SUN_PILLBOX)
     return RES_BAD_ARG;
   sun->data.pillbox.aperture = angle;