Index: /issm/trunk-jpl/src/c/shared/Random/randomgenerator.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/randomgenerator.cpp	(revision 26554)
+++ /issm/trunk-jpl/src/c/shared/Random/randomgenerator.cpp	(revision 26554)
@@ -0,0 +1,6 @@
+/*!\file: randomgenerator
+ * \brief random number generating functions
+ */
+
+#include <iostream>
+#include "./randomgenerator.h"
Index: /issm/trunk-jpl/src/c/shared/Random/randomgenerator.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/randomgenerator.h	(revision 26554)
+++ /issm/trunk-jpl/src/c/shared/Random/randomgenerator.h	(revision 26554)
@@ -0,0 +1,64 @@
+/*!\file: randomgenerator.h
+ * \brief prototypes for randomgenerator.h
+ */
+
+#ifndef _RANDOMGENERATOR_H_
+#define _RANDOMGENERATOR_H_
+
+#undef M_PI
+#define M_PI 3.141592653589793238462643
+
+class uniform_distribution_rnd
+{
+
+  private:
+    int a;
+    int c;
+    unsigned int m;
+    unsigned _seed;
+    double a1;
+    double a2;
+
+    int rnd_int() { return( _seed = ( a * _seed + c ) % m ); }
+
+  public:
+    uniform_distribution_rnd() : _seed( 0 ), a( 1103515245 ), c( 12345 ), m( 2147483648 ), a1(0.0), a2(1.0) {}
+    uniform_distribution_rnd(double a_1,double a_2) : _seed( 0 ), a( 1103515245 ), c( 12345 ), m( 2147483648 ), a1(a_1), a2(a_2) {}
+    void seed( unsigned int s ) { _seed = s; }
+    unsigned int get_seed() { return _seed; }
+    double rnd() { return (a2-a1)*(double) rnd_int()/ m + a1; }
+
+};
+
+class normal_distribution_rnd
+{
+
+  private:
+    unsigned _seed;
+    double mean;
+    double sdev;
+
+  public:
+    normal_distribution_rnd() : _seed( 0 ), mean( 0), sdev(1.0) {}
+    normal_distribution_rnd(double m,double s) : _seed( 0 ), mean( m ), sdev(s) {}
+    void seed( unsigned int s ) { _seed = s; }
+    double rnd()
+    {
+      uniform_distribution_rnd unifdistri;
+      unifdistri.seed(_seed);
+
+      double u1 = unifdistri.rnd();
+      double u2 = unifdistri.rnd();
+
+      double R = sqrt(-2*log(u1));
+      double theta = 2*M_PI*u2;
+
+      seed(unifdistri.get_seed());
+
+      return mean + sdev * (R*cos(theta));
+
+    }
+
+};
+
+#endif //ifndef _RANDOMGENERATOR_H_
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp	(revision 26553)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp	(revision 26554)
@@ -8,37 +8,6 @@
 #include "../modules/modules.h"
 #include "../analyses/analyses.h"
+#include "../shared/Random/randomgenerator.h"
 #include <random>
-
-#undef M_PI
-#define M_PI 3.141592653589793238462643
-
-double lcg(long long* seed){/*{{{*/
-
-
-	int N = 5;
-	long long M = 1LL<<32;
-	long long a = 1812433253;
-	long long c = 1;
-
-	for(int i=0;i<N;i++){
-		*seed = a * *seed + c;
-		*seed = *seed%M;
-	}
-
-  return ((double) *seed) / ((double) M);
-
-}
-
-double BoxMuller(long long* pseed){
-
-	double u1 = lcg(pseed);
-	double u2 = lcg(pseed);
-
-	double R = sqrt(-2*log(u1));
-	double theta = 2*M_PI*u2;
-
-	return R*cos(theta);
-
-}
 
 void GaussianVector(Vector<IssmDouble>* ppf,int seed){/*{{{*/
@@ -49,4 +18,5 @@
 
 	/*Define seed*/
+	normal_distribution_rnd distribution;
 	if(seed<0){
 		std::random_device rd;
@@ -59,7 +29,5 @@
 		seed = seed + 783728*my_rank; // change default seed for parallel simulations (by considering an arbitrary shif based on the rank number)
 	}
-
-	long long seedbis;
- 	seedbis  = (long long) seed;
+	distribution.seed(seed);
 
 	int        *local_indices = NULL;
@@ -70,5 +38,5 @@
 	ppf->GetLocalSize(&M);
 	for(int i=0;i<M;i++){
-		rdnumber = BoxMuller(&seedbis);
+		rdnumber = distribution.rnd();
 		ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
 	}