From f833d348e0337947c3852b591f4af98b5e2c7cdc Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 17 Dec 2021 08:01:15 +0100
Subject: [PATCH] testing random sampling

---
 PRIVATE      |  2 +-
 src/math.f90 | 28 ++++++++++++++++++++++++----
 2 files changed, 25 insertions(+), 5 deletions(-)

diff --git a/PRIVATE b/PRIVATE
index 2ad27552c..96c32ba42 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 2ad27552c43316735b6ef425737fe3c8a5231598
+Subproject commit 96c32ba4237a51eaad92cd139e1a716ee5b32493
diff --git a/src/math.f90 b/src/math.f90
index 3303f3b4c..2d8f8fb3b 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -895,7 +895,7 @@ pure function math_33toVoigt6_stress(sigma) result(sigma_tilde)
 
 
   sigma_tilde = [sigma(1,1), sigma(2,2), sigma(3,3), &
-                 sigma(3,2), sigma(3,1), sigma(1,2)] 
+                 sigma(3,2), sigma(3,1), sigma(1,2)]
 
 end function math_33toVoigt6_stress
 
@@ -910,7 +910,7 @@ pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)
 
 
   epsilon_tilde = [          epsilon(1,1),           epsilon(2,2),           epsilon(3,3), &
-                   2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)] 
+                   2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)]
 
 end function math_33toVoigt6_strain
 
@@ -970,11 +970,11 @@ impure elemental subroutine math_normal(x,mu,sigma)
 
   real(pReal), intent(out) :: x
   real(pReal), intent(in), optional :: mu, sigma
-  
+
   real(pReal) :: sigma_, mu_
   real(pReal), dimension(2) :: rnd
 
-  
+
   if (present(mu)) then
     mu_ = mu
   else
@@ -1431,6 +1431,26 @@ subroutine selfTest
   if (dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) &
     error stop 'math_LeviCivita'
 
+  normal_distribution: block
+    real(pReal), dimension(500000) :: r
+    real(pReal) :: mu, sigma
+
+    call random_number(mu)
+    call random_number(sigma)
+
+    sigma = 1.0_pReal + sigma*5.0_pReal
+    mu = (mu-0.5_pReal)*10_pReal
+
+    call math_normal(r,mu,sigma)
+
+    if (abs(mu -sum(r)/real(size(r),pReal))>5.0e-2_pReal) &
+      error stop 'math_normal(mu)'
+
+    mu = sum(r)/real(size(r),pReal)
+    if (abs(sigma**2 -1.0_pReal/real(size(r)-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
+      error stop 'math_normal(sigma)'
+  end block normal_distribution
+
 end subroutine selfTest
 
 end module math