From a869754bc425db63fc23d5f3c4aafa175ee2bad4 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 8 Apr 2015 13:08:08 +0100
Subject: [PATCH] StochasticDispersionRAS: Corrected spherical distribution of
 UTurb direction Resolves bug-report
 http://www.openfoam.org/mantisbt/view.php?id=1650

---
 .../StochasticDispersionRAS.C                 | 36 +++++++++++--------
 .../StochasticDispersionRAS.C                 | 14 +++-----
 2 files changed, 26 insertions(+), 24 deletions(-)

diff --git a/src/lagrangian/Turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C b/src/lagrangian/Turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
index badebbc1..8da932eb 100644
--- a/src/lagrangian/Turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
+++ b/src/lagrangian/Turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "StochasticDispersionRAS.H"
+#include "constants.H"
+
+using namespace Foam::constant::mathematical;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -89,29 +92,34 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
 
         if (tTurb > tTurbLoc)
         {
-            tTurb = 0.0;
+            tTurb = 0;
 
-            scalar sigma = sqrt(2.0*k/3.0);
-            vector dir = 2.0*rnd.sample01<vector>() - vector::one;
-            dir /= mag(dir) + SMALL;
+            const scalar sigma = sqrt(2*k/3.0);
+
+            // Calculate a random direction dir distributed uniformly
+            // in spherical coordinates
+
+            const scalar theta = rnd.sample01<scalar>()*twoPi;
+            const scalar u = 2*rnd.sample01<scalar>() - 1;
+
+            const scalar a = sqrt(1 - sqr(u));
+            const vector dir(a*cos(theta), a*sin(theta), u);
 
             // Numerical Recipes... Ch. 7. Random Numbers...
-            scalar x1 = 0.0;
-            scalar x2 = 0.0;
-            scalar rsq = 10.0;
-            while ((rsq > 1.0) || (rsq == 0.0))
+            scalar x1 = 0;
+            scalar x2 = 0;
+            scalar rsq = 10;
+            while ((rsq > 1) || (rsq == 0))
             {
-                x1 = 2.0*rnd.sample01<scalar>() - 1.0;
-                x2 = 2.0*rnd.sample01<scalar>() - 1.0;
+                x1 = 2*rnd.sample01<scalar>() - 1;
+                x2 = 2*rnd.sample01<scalar>() - 1;
                 rsq = x1*x1 + x2*x2;
             }
 
-            scalar fac = sqrt(-2.0*log(rsq)/rsq);
-
+            scalar fac = sqrt(-2*log(rsq)/rsq);
             fac *= mag(x1);
 
             UTurb = sigma*fac*dir;
-
         }
     }
     else
diff --git a/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C b/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
index 165812f1..8da932eb 100644
--- a/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
+++ b/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
@@ -99,17 +99,11 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
             // Calculate a random direction dir distributed uniformly
             // in spherical coordinates
 
-            const scalar theta = rnd.sample01<scalar>()*pi;
-            const scalar phi = rnd.sample01<scalar>()*twoPi;
+            const scalar theta = rnd.sample01<scalar>()*twoPi;
+            const scalar u = 2*rnd.sample01<scalar>() - 1;
 
-            // Optimising compilers will use the sincos function
-            const scalar sinTheta = sin(theta);
-            const scalar cosTheta = cos(theta);
-
-            const scalar sinPhi = sin(phi);
-            const scalar cosPhi = cos(phi);
-
-            const vector dir(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta);
+            const scalar a = sqrt(1 - sqr(u));
+            const vector dir(a*cos(theta), a*sin(theta), u);
 
             // Numerical Recipes... Ch. 7. Random Numbers...
             scalar x1 = 0;