solstice-solver

commit aef7d5b3ef69e8b181fbfe7fdbe94cb8c4fec549
parent f5c7ea7570a9420e8932768c4e6c6e5f2a76ea4d
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue, 19 Dec 2017 16:09:28 +0100

Fix gaussian sunshape

Diffstat:
M
src/ssol_ranst_sun_dir.c
 | 
12
++++++++----

1 file changed, 8 insertions(+), 4 deletions(-)
diff --git a/src/ssol_ranst_sun_dir.c b/src/ssol_ranst_sun_dir.c
@@ -275,15 +275,19 @@ ran_gaussian_get
    double dir[3])
 {
   double pt[3];
-  double phi, theta, sin_theta;
+  double phi, cos_theta, sin_theta;
 
+  /* The following is not truly a gaussian sunshape,
+   * but an accurate enough approximation for small angles */
   ASSERT(ran && 0 <= ran->state.gaussian.std_dev);
-  theta = ssp_ran_gaussian(rng, 0, ran->state.gaussian.std_dev);
-  sin_theta = sin(theta);
+  sin_theta
+    = ran->state.gaussian.std_dev * sqrt(-2 * log(ssp_rng_canonical(rng)));
+  sin_theta = MMIN(1, sin_theta); /* macro: don't merge with previous line! */
+  cos_theta = sin2cos(sin_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);
+  pt[2] = cos_theta;
   d33_muld3(dir, ran->state.gaussian.basis, pt);
   d3_normalize(dir, dir);
   return dir;