Remove RNG lib and use gsl for all random number generation.

The GAUSS distribution is now using a Choleski decomposition 
to correlate. One can use either
CORRX,CORRY, CORRZ,R16,R26, R36, R46
or the correlation matrix R a la TRANSPORT.

src/Distribution/CMakeLists.txt

 set (_SRCS
   Distribution.cpp
-  ranlib.cc
   halton1d_sequence.cpp
   LaserProfile.cpp
   )

src/Distribution/Distribution.h

@@ -25,7 +25,7 @@
 
 #include "AbstractObjects/Definition.h"
 #include "Algorithms/PartData.h"
-#include "ranlib.h"
+
 #include "Algorithms/Vektor.h"
 
 #include "Ippl.h"
@@ -36,8 +36,6 @@
 #include <gsl/gsl_histogram.h>
 #include <gsl/gsl_qrng.h>
 
-#define RANLIBX
-
 class Beam;
 class PartBunch;
 class PartBins;
@@ -215,6 +213,7 @@ private:
     void GenerateFlattopT(size_t numberOfParticles);
     void GenerateFlattopZ(size_t numberOfParticles);
     void GenerateGaussZ(size_t numberOfParticles);
+    void GenerateGaussZChol(size_t numberOfParticles);
     void GenerateLongFlattopT(size_t numberOfParticles);
     void GenerateTransverseGauss(size_t numberOfParticles);
     void InitializeBeam(PartBunch &beam);
@@ -287,7 +286,7 @@ private:
     gsl_histogram *energyBinHist_m; /// GSL histogram used to define energy bin
                                     /// structure.
 
-    RANLIB_class *randGenEmit_m;    /// Random number generator for thermal emission
+    gsl_rng *randGenEmit_m;         /// Random number generator for thermal emission
                                     /// models.
 
     // ASTRA and NONE photo emission model.

src/Structure/BoundaryGeometry.cpp

@@ -109,6 +109,9 @@ BoundaryGeometry::BoundaryGeometry() :
     }
     if (!h5FileName_m.empty ())
         initialize ();
+
+    gsl_rng_env_setup();
+    randGen_m = gsl_rng_alloc(gsl_rng_default);    
 }
 
 BoundaryGeometry::BoundaryGeometry(
@@ -129,6 +132,8 @@ BoundaryGeometry::BoundaryGeometry(
     TRayTrace_m = IpplTimings::getTimer ("Ray tracing");
     Tinward_m = IpplTimings::getTimer ("Check inward");
 
+    gsl_rng_env_setup();
+    randGen_m = gsl_rng_alloc(gsl_rng_default);    
 }
 
 BoundaryGeometry::~BoundaryGeometry() {
@@ -144,6 +149,10 @@ BoundaryGeometry::~BoundaryGeometry() {
        if (TriSePartloss_m)
        delete TriSePartloss_m;
      */
+
+    gsl_rng_free(randGen_m);
+
+
 }
 
 bool BoundaryGeometry::canReplaceBy (Object* object) {

src/Structure/BoundaryGeometry.h

@@ -29,6 +29,8 @@ class ElementBase;
 #include "Distribution/ranlib.h"
 #include "Structure/SecondaryEmissionPhysics.h"
 
+#include <gsl/gsl_rng.h>
+
 extern Inform* gmsg;
 
 namespace BGphysics {
@@ -421,6 +423,8 @@ private:
     double parameterFNVYZe_m;   // zero order fit constant for v(y). Default:\f$0.9632\f$.
     double parameterFNVYSe_m;   // second order fit constant for v(y). Default:\f$1.065\f$.
 
+    gsl_rng *randGen_m;         // 
+  
     IpplTimings::TimerRef TPInside_m;   // timers for profiling
     IpplTimings::TimerRef TPreProc_m;
     IpplTimings::TimerRef TRayTrace_m;
@@ -540,8 +544,10 @@ private:
         }
     }
 
-    /*
-       Determine if a point x is outside, inside or just on the boundary.
+
+    inline bool isInside (Vector_t x) {
+        /*
+          DetBermine if a point x is outside, inside or just on the boundary.
        Return true if point is inside boundary otherwise false.
 
        The basic idea is if a line segment starting from the test point has
@@ -553,17 +559,13 @@ private:
        on the boundary. Makesure the end point of the line
        segment is outside the geometry boundary.
     */
-    bool isInside (Vector_t x) {
+
         Vector_t x0 = x;
         Vector_t x1;
         x1[0] = x0[0];
-        //x1[1] = x0[1];
-        RANLIB_class* rGen = new RANLIB_class (265314159, 4);
-        x1[1] = maxcoords_m[1] * (1.1 + rGen->uniform (0.0, 1.0));
-        x1[2] = maxcoords_m[2] * (1.1 + rGen->uniform (0.0, 1.0));
-        //x1[2] = x0[2];
-        delete rGen;
-	
+        x1[1] = maxcoords_m[1] * (1.1 + gsl_rng_uniform(randGen_m));
+        x1[2] = maxcoords_m[2] * (1.1 + gsl_rng_uniform(randGen_m));
+
         /*
           Random number could avoid some specific situation,
           like line parallel to boundary......
@@ -582,7 +584,6 @@ private:
         }
     }
 
-
     /*
       Recursively get inward triangle normal of all surface triangles.