--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/samples/main-random-variable.cc	Thu Apr 16 10:07:48 2009 +0200
@@ -0,0 +1,377 @@
+/* -*-  Mode: C++; c-file-style: "gnu"; indent-tabs-mode:nil; -*- */
+/*
+ * Copyright (c) 2008 Timo Bingmann
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License version 2 as
+ * published by the Free Software Foundation;
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *
+ * Author: Timo Bingmann <timo.bingmann@student.kit.edu>
+ */
+#include "ns3/random-variable.h"
+#include "ns3/gnuplot.h"
+#include <map>
+#include <cmath>
+
+using namespace ns3;
+
+/// Round a double number to the given precision. e.g. dround(0.234, 0.1) = 0.2
+/// and dround(0.257, 0.1) = 0.3
+double dround(double number, double precision)
+{
+  number /= precision;
+  if (number >= 0)
+    number = std::floor(number + 0.5);
+  else
+    number = std::ceil(number - 0.5);
+  number *= precision;
+  return number;
+}
+
+static GnuplotDataset
+Histogramm (RandomVariable rndvar, unsigned int probes, double precision, const std::string& title, bool notcontinous = false)
+{
+  typedef std::map<double, unsigned int> histogramm_maptype;
+  histogramm_maptype histogramm;
+
+  for(unsigned int i = 0; i < probes; ++i)
+    {
+      double val = dround( rndvar.GetValue(), precision );
+
+      ++histogramm[val];
+    }
+
+  Gnuplot2dDataset data;
+  data.SetTitle(title);
+
+  if (notcontinous)
+    {
+      data.SetStyle(Gnuplot2dDataset::IMPULSES);
+    }
+
+  for(histogramm_maptype::const_iterator hi = histogramm.begin();
+      hi != histogramm.end(); ++hi)
+    {
+      data.Add(hi->first, (double)hi->second / (double)probes / precision);
+    }
+
+  return data;
+}
+
+int main (int argc, char *argv[])
+{
+  unsigned int probes = 1000000;
+  double precision = 0.01;
+
+  GnuplotCollection gnuplots("main-random-variables.pdf");
+  gnuplots.SetTerminal("pdf enhanced");
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("UniformVariable");
+    plot.AppendExtra("set yrange [0:]");
+
+    plot.AddDataset( Histogramm(UniformVariable(0.0, 1.0), probes, precision,
+                                "UniformVariable [0.0 .. 1.0)") );
+
+    plot.AddDataset( Gnuplot2dFunction("0.1",
+                                       "0 <= x && x <= 1 ? 1.0 : 0") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("ExponentialVariable");
+    plot.AppendExtra("set xrange [0:8]");
+    plot.AppendExtra("ExpDist(x,l) = 1/l * exp(-1/l * x)");
+
+    plot.AddDataset( Histogramm(ExponentialVariable(0.5), probes, precision,
+                                "ExponentialVariable m=0.5") );
+
+    plot.AddDataset( Gnuplot2dFunction("ExponentialDistribution mean 0.5",
+                                       "ExpDist(x, 0.5)") );
+
+    plot.AddDataset( Histogramm(ExponentialVariable(1.0), probes, precision,
+                                "ExponentialVariable m=1") );
+
+    plot.AddDataset( Gnuplot2dFunction("ExponentialDistribution mean 1.0",
+                                       "ExpDist(x, 1.0)") );
+
+    plot.AddDataset( Histogramm(ExponentialVariable(1.5), probes, precision,
+                                "ExponentialVariable m=1.5") );
+
+    plot.AddDataset( Gnuplot2dFunction("ExponentialDistribution mean 1.5",
+                                       "ExpDist(x, 1.5)") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("ParetoVariable");
+    plot.AppendExtra("set xrange [0:2]");
+
+    plot.AddDataset( Histogramm(ParetoVariable(1.0, 1.5), probes, precision,
+                                "ParetoVariable m=1.0 s=1.5") );
+
+    plot.AddDataset( Histogramm(ParetoVariable(1.0, 2.0), probes, precision,
+                                "ParetoVariable m=1.0 s=2.0") );
+
+    plot.AddDataset( Histogramm(ParetoVariable(1.0, 2.5), probes, precision,
+                                "ParetoVariable m=1.0 s=2.5") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("WeibullVariable");
+    plot.AppendExtra("set xrange [0:3]");
+
+    plot.AddDataset( Histogramm(WeibullVariable(1.0, 1.0), probes, precision,
+                                "WeibullVariable m=1.0 s=1.0") );
+
+    plot.AddDataset( Histogramm(WeibullVariable(1.0, 2.0), probes, precision,
+                                "WeibullVariable m=1.0 s=2.0") );
+
+    plot.AddDataset( Histogramm(WeibullVariable(1.0, 3.0), probes, precision,
+                                "WeibullVariable m=1.0 s=3.0") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("NormalVariable");
+    plot.AppendExtra("set xrange [-3:3]");
+    plot.AppendExtra("NormalDist(x,m,s) = 1 / (s * sqrt(2*pi)) * exp(-1.0 / 2.0 * ((x-m) / s)**2)");
+
+    plot.AddDataset( Histogramm(NormalVariable(0.0, 1.0), probes, precision,
+                                "NormalVariable m=0.0 v=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("NormalDist {/Symbol m}=0.0 {/Symbol s}=1.0",
+                                       "NormalDist(x,0.0,1.0)") );
+
+    plot.AddDataset( Histogramm(NormalVariable(0.0, 2.0), probes, precision,
+                                "NormalVariable m=0.0 v=2.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("NormalDist {/Symbol m}=0.0 {/Symbol s}=sqrt(2.0)",
+                                       "NormalDist(x,0.0,sqrt(2.0))") );
+
+    plot.AddDataset( Histogramm(NormalVariable(0.0, 3.0), probes, precision,
+                                "NormalVariable m=0.0 v=3.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("NormalDist {/Symbol m}=0.0 {/Symbol s}=sqrt(3.0)",
+                                       "NormalDist(x,0.0,sqrt(3.0))") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("EmpiricalVariable");
+    plot.AppendExtra("set xrange [*:*]");
+
+    EmpiricalVariable emp1;
+    emp1.CDF(0.0,  0.0 / 15.0);
+    emp1.CDF(0.2,  1.0 / 15.0);
+    emp1.CDF(0.4,  3.0 / 15.0);
+    emp1.CDF(0.6,  6.0 / 15.0);
+    emp1.CDF(0.8, 10.0 / 15.0);
+    emp1.CDF(1.0, 15.0 / 15.0);
+
+    plot.AddDataset( Histogramm(emp1, probes, precision,
+                                "EmpiricalVariable (Stairs)") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("DeterministicVariable");
+    plot.AppendExtra("set xrange [*:*]");
+
+    double values[] = { 0.0, 0.2, 0.2, 0.4, 0.2, 0.6, 0.8, 0.8, 1.0 };
+    DeterministicVariable det1 (values, sizeof(values) / sizeof(values[0]));
+
+    plot.AddDataset( Histogramm(det1, probes, precision,
+                                "DeterministicVariable", true) );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("LogNormalVariable");
+    plot.AppendExtra("set xrange [0:3]");
+
+    plot.AppendExtra("LogNormalDist(x,m,s) = 1.0/x * NormalDist(log(x), m, s)");
+
+    plot.AddDataset( Histogramm(LogNormalVariable(0.0, 1.0), probes, precision,
+                                "LogNormalVariable m=0.0 s=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("LogNormalDist(x, 0.0, 1.0)",
+                                       "LogNormalDist(x, 0.0, 1.0)") );
+
+    plot.AddDataset( Histogramm(LogNormalVariable(0.0, 0.5), probes, precision,
+                                "LogNormalVariable m=0.0 s=0.5") );
+
+    plot.AddDataset( Histogramm(LogNormalVariable(0.0, 0.25), probes, precision,
+                                "LogNormalVariable m=0.0 s=0.25") );
+
+    plot.AddDataset( Gnuplot2dFunction("LogNormalDist(x, 0.0, 0.25)",
+                                       "LogNormalDist(x, 0.0, 0.25)") );
+
+    plot.AddDataset( Histogramm(LogNormalVariable(0.0, 0.125), probes, precision,
+                                "LogNormalVariable m=0.0 s=0.125") );
+
+    plot.AddDataset( Histogramm(LogNormalVariable(0.0, 2.0), probes, precision,
+                                "LogNormalVariable m=0.0 s=2.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("LogNormalDist(x, 0.0, 2.0)",
+                                       "LogNormalDist(x, 0.0, 2.0)") );
+
+    plot.AddDataset( Histogramm(LogNormalVariable(0.0, 2.5), probes, precision,
+                                "LogNormalVariable m=0.0 s=2.5") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("TriangularVariable");
+    plot.AppendExtra("set xrange [*:*]");
+
+    plot.AddDataset( Histogramm(TriangularVariable(0.0, 1.0, 0.5), probes, precision,
+                                "TriangularVariable [0.0 .. 1.0) m=0.5") );
+
+    plot.AddDataset( Histogramm(TriangularVariable(0.0, 1.0, 0.4), probes, precision,
+                                "TriangularVariable [0.0 .. 1.0) m=0.4") );
+
+    plot.AddDataset( Histogramm(TriangularVariable(0.0, 1.0, 0.65), probes, precision,
+                                "TriangularVariable [0.0 .. 1.0) m=0.65") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("GammaVariable");
+    plot.AppendExtra("set xrange [0:10]");
+    plot.AppendExtra("set yrange [0:1]");
+    plot.AppendExtra("GammaDist(x,a,b) = x**(a-1) * 1/b**a * exp(-x/b) / gamma(a)");
+
+    plot.AppendExtra("set label 1 '{/Symbol g}(x,{/Symbol a},{/Symbol b}) = x^{/Symbol a-1} e^{-x {/Symbol b}^{-1}} ( {/Symbol b}^{/Symbol a} {/Symbol G}({/Symbol a}) )^{-1}' at 0.7, 0.9");
+
+    plot.AddDataset( Histogramm(GammaVariable(1.0, 1.0), probes, precision,
+                                "GammaVariable a=1.0 b=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("{/Symbol g}(x, 1.0, 1.0)",
+                                       "GammaDist(x, 1.0, 1.0)") );
+
+    plot.AddDataset( Histogramm(GammaVariable(1.5, 1.0), probes, precision,
+                                "GammaVariable a=1.5 b=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("{/Symbol g}(x, 1.5, 1.0)",
+                                       "GammaDist(x, 1.5, 1.0)") );
+
+    plot.AddDataset( Histogramm(GammaVariable(2.0, 1.0), probes, precision,
+                                "GammaVariable a=2.0 b=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("{/Symbol g}(x, 2.0, 1.0)",
+                                       "GammaDist(x, 2.0, 1.0)") );
+
+    plot.AddDataset( Histogramm(GammaVariable(4.0, 1.0), probes, precision,
+                                "GammaVariable a=4.0 b=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("{/Symbol g}(x, 4.0, 1.0)",
+                                       "GammaDist(x, 4.0, 1.0)") );
+
+    plot.AddDataset( Histogramm(GammaVariable(2.0, 2.0), probes, precision,
+                                "GammaVariable a=2.0 b=2.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("{/Symbol g}(x, 2.0, 2.0)",
+                                       "GammaDist(x, 2.0, 2.0)") );
+
+    plot.AddDataset( Histogramm(GammaVariable(2.5, 3.0), probes, precision,
+                                "GammaVariable a=2.5 b=3.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("{/Symbol g}(x, 2.5, 3.0)",
+                                       "GammaDist(x, 2.5, 3.0)") );
+
+    plot.AddDataset( Histogramm(GammaVariable(2.5, 4.5), probes, precision,
+                                "GammaVariable a=2.5 b=4.5") );
+
+    plot.AddDataset( Gnuplot2dFunction("{/Symbol g}(x, 2.5, 4.5)",
+                                       "GammaDist(x, 2.5, 4.5)") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  {
+    Gnuplot plot;
+    plot.SetTitle("ErlangVariable");
+    plot.AppendExtra("set xrange [0:10]");
+    plot.AppendExtra("ErlangDist(x,k,l) = x**(k-1) * 1/l**k * exp(-x/l) / (k-1)!");
+
+    plot.AppendExtra("set label 1 'Erlang(x,k,{/Symbol l}) = x^{k-1} e^{-x {/Symbol l}^{-1}} ( {/Symbol l}^k (k-1)! )^{-1}' at 0.7, 0.9");
+
+    plot.AddDataset( Histogramm(ErlangVariable(1, 1.0), probes, precision,
+                                "ErlangVariable k=1 {/Symbol l}=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("Erlang(x, 1, 1.0)",
+                                       "ErlangDist(x, 1, 1.0)") );
+
+    plot.AddDataset( Histogramm(ErlangVariable(2, 1.0), probes, precision,
+                                "ErlangVariable k=2 {/Symbol l}=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("Erlang(x, 2, 1.0)",
+                                       "ErlangDist(x, 2, 1.0)") );
+
+    plot.AddDataset( Histogramm(ErlangVariable(3, 1.0), probes, precision,
+                                "ErlangVariable k=3 {/Symbol l}=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("Erlang(x, 3, 1.0)",
+                                       "ErlangDist(x, 3, 1.0)") );
+
+    plot.AddDataset( Histogramm(ErlangVariable(5, 1.0), probes, precision,
+                                "ErlangVariable k=5 {/Symbol l}=1.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("Erlang(x, 5, 1.0)",
+                                       "ErlangDist(x, 5, 1.0)") );
+
+    plot.AddDataset( Histogramm(ErlangVariable(2, 2.0), probes, precision,
+                                "ErlangVariable k=2 {/Symbol l}=2.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("Erlang(x, 2, 2.0)",
+                                       "ErlangDist(x, 2, 2.0)") );
+
+    plot.AddDataset( Histogramm(ErlangVariable(2, 3.0), probes, precision,
+                                "ErlangVariable k=2 {/Symbol l}=3.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("Erlang(x, 2, 3.0)",
+                                       "ErlangDist(x, 2, 3.0)") );
+
+    plot.AddDataset( Histogramm(ErlangVariable(2, 5.0), probes, precision,
+                                "ErlangVariable k=2 {/Symbol l}=5.0") );
+
+    plot.AddDataset( Gnuplot2dFunction("Erlang(x, 2, 5.0)",
+                                       "ErlangDist(x, 2, 5.0)") );
+
+    gnuplots.AddPlot(plot);
+  }
+
+  gnuplots.GenerateOutput(std::cout);
+
+  return 0;
+}

--- a/samples/wscript	Thu Apr 16 09:25:32 2009 +0200
+++ b/samples/wscript	Thu Apr 16 10:07:48 2009 +0200
@@ -10,6 +10,9 @@
     obj = bld.create_ns3_program('main-simulator')
     obj.source = 'main-simulator.cc'
 
+    obj = bld.create_ns3_program('main-random-variable')
+    obj.source = 'main-random-variable.cc'
+
     obj = bld.create_ns3_program('main-packet-header', ['common', 'simulator'])
     obj.source = 'main-packet-header.cc'
 

--- a/src/core/random-variable.cc	Thu Apr 16 09:25:32 2009 +0200
+++ b/src/core/random-variable.cc	Thu Apr 16 10:07:48 2009 +0200
@@ -1220,6 +1220,230 @@
 
 //-----------------------------------------------------------------------------
 //-----------------------------------------------------------------------------
+// GammaVariableImpl
+class GammaVariableImpl : public RandomVariableBase
+{
+public:
+  /**
+   * \param alpha alpha parameter of the gamma distribution
+   * \param beta beta parameter of the gamma distribution
+   */
+  GammaVariableImpl (double alpha, double beta);
+
+  /**
+   * \return A random value from this distribution
+   */
+  virtual double GetValue ();
+
+  /**
+   * \return A random value from the gamma distribution with parameters alpha
+   * and beta
+   */
+  double GetValue(double alpha, double beta);
+
+  virtual RandomVariableBase* Copy(void) const;
+
+private:
+  double m_alpha;
+  double m_beta;
+  NormalVariable m_normal;
+};
+
+
+RandomVariableBase* GammaVariableImpl::Copy () const
+{
+  return new GammaVariableImpl (m_alpha, m_beta);
+}
+
+GammaVariableImpl::GammaVariableImpl (double alpha, double beta)
+  : m_alpha(alpha), m_beta(beta) 
+{
+}
+
+double
+GammaVariableImpl::GetValue ()
+{
+  return GetValue(m_alpha, m_beta);
+}
+
+/*
+  The code for the following generator functions was adapted from ns-2
+  tools/ranvar.cc
+
+  Originally the algorithm was devised by Marsaglia in 2000:
+  G. Marsaglia, W. W. Tsang: A simple method for gereating Gamma variables
+  ACM Transactions on mathematical software, Vol. 26, No. 3, Sept. 2000
+
+  The Gamma distribution density function has the form 
+
+	                     x^(alpha-1) * exp(-x/beta)
+	p(x; alpha, beta) = ----------------------------
+                             beta^alpha * Gamma(alpha)
+
+  for x > 0.
+*/
+double
+GammaVariableImpl::GetValue (double alpha, double beta)
+{
+  if(!m_generator)
+    {
+      m_generator = new RngStream();
+    }
+
+  if (alpha < 1)
+    {
+      double u = m_generator->RandU01 ();
+      return GetValue(1.0 + alpha, beta) * pow (u, 1.0 / alpha);
+    }
+	
+  double x, v, u;
+  double d = alpha - 1.0 / 3.0;
+  double c = (1.0 / 3.0) / sqrt (d);
+
+  while (1)
+    {
+      do
+        {
+          x = m_normal.GetValue ();
+          v = 1.0 + c * x;
+        } while (v <= 0);
+
+      v = v * v * v;
+      u = m_generator->RandU01 ();
+      if (u < 1 - 0.0331 * x * x * x * x)
+        {
+          break;
+        }
+      if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
+        {
+          break;
+        }
+    }
+
+  return beta * d * v;
+}
+
+GammaVariable::GammaVariable ()
+  : RandomVariable (GammaVariableImpl (1.0, 1.0))
+{
+}
+
+GammaVariable::GammaVariable (double alpha, double beta)
+  : RandomVariable (GammaVariableImpl (alpha, beta))
+{
+}
+
+double GammaVariable::GetValue(void) const
+{
+  return this->RandomVariable::GetValue ();
+}
+
+double GammaVariable::GetValue(double alpha, double beta) const
+{
+  return ((GammaVariableImpl*)Peek())->GetValue(alpha, beta);
+}
+
+//-----------------------------------------------------------------------------
+//-----------------------------------------------------------------------------
+// ErlangVariableImpl
+
+class ErlangVariableImpl : public RandomVariableBase
+{
+public:
+  /**
+   * \param k k parameter of the Erlang distribution
+   * \param lambda lambda parameter of the Erlang distribution
+   */
+  ErlangVariableImpl (unsigned int k, double lambda);
+
+  /**
+   * \return A random value from this distribution
+   */
+  virtual double GetValue ();
+
+  /**
+   * \return A random value from the Erlang distribution with parameters k and
+   * lambda.
+   */
+  double GetValue(unsigned int k, double lambda);
+
+  virtual RandomVariableBase* Copy(void) const;
+
+private:
+  unsigned int m_k;
+  double m_lambda;
+};
+
+
+RandomVariableBase* ErlangVariableImpl::Copy () const
+{
+  return new ErlangVariableImpl (m_k, m_lambda);
+}
+
+ErlangVariableImpl::ErlangVariableImpl (unsigned int k, double lambda)
+  : m_k(k), m_lambda(lambda)
+{
+}
+
+double
+ErlangVariableImpl::GetValue ()
+{
+  return GetValue(m_k, m_lambda);
+}
+
+/*
+  The code for the following generator functions was adapted from ns-2
+  tools/ranvar.cc
+
+  The Erlang distribution density function has the form 
+
+	                   x^(k-1) * exp(-x/lambda)
+	p(x; k, lambda) = ---------------------------
+                             lambda^k * (k-1)!
+
+  for x > 0.
+*/
+double
+ErlangVariableImpl::GetValue (unsigned int k, double lambda)
+{
+  if(!m_generator)
+    {
+      m_generator = new RngStream();
+    }
+
+  ExponentialVariable exponential(lambda);
+
+  double result = 0;
+  for (unsigned int i = 0; i < k; ++i)
+    {
+      result += exponential.GetValue();
+    }
+
+  return result;
+}
+
+ErlangVariable::ErlangVariable ()
+  : RandomVariable (ErlangVariableImpl (1, 1.0))
+{
+}
+
+ErlangVariable::ErlangVariable (unsigned int k, double lambda)
+  : RandomVariable (ErlangVariableImpl (k, lambda))
+{
+}
+
+double ErlangVariable::GetValue(void) const
+{
+  return this->RandomVariable::GetValue ();
+}
+
+double ErlangVariable::GetValue(unsigned int k, double lambda) const
+{
+  return ((ErlangVariableImpl*)Peek())->GetValue(k, lambda);
+}
+
+//-----------------------------------------------------------------------------
+//-----------------------------------------------------------------------------
 // TriangularVariableImpl methods
 class TriangularVariableImpl : public RandomVariableBase {
 public:

--- a/src/core/random-variable.h	Thu Apr 16 09:25:32 2009 +0200
+++ b/src/core/random-variable.h	Thu Apr 16 10:07:48 2009 +0200
@@ -36,7 +36,7 @@
  *
  */
 
-namespace ns3{
+namespace ns3 {
 
 class RandomVariableBase;
 
@@ -422,7 +422,7 @@
    * value and a shape (alpha) parameter of 1.5.
    * \param m mean value of the distribution
    */
-   WeibullVariable(double m) ;
+  WeibullVariable(double m) ;
 
   /**
    * Constructs a weibull random variable with the specified mean
@@ -560,7 +560,6 @@
   explicit DeterministicVariable(double* d, uint32_t c);
 };
 
-
 /**
  * \brief Log-normal Distributed random var
  * \ingroup randomvariable
@@ -594,6 +593,108 @@
 };
 
 /**
+ * \brief Gamma Distributed Random Variable
+ * \ingroup randomvariable
+ *
+ * GammaVariable defines a random variable with gamma distribution. 
+ *
+ * This class supports the creation of objects that return random numbers
+ * from a fixed gamma distribution. It also supports the generation of 
+ * single random numbers from various gamma distributions.
+ *
+ * The probability density function is defined over the interval [0,+inf) as:
+ * \f$ x^{\alpha-1} \frac{e^{-\frac{x}{\beta}}}{\beta^\alpha \Gamma(\alpha)}\f$
+ * where \f$ mean = \alpha\beta \f$ and 
+ * \f$ variance = \alpha \beta^2\f$
+ */
+class GammaVariable : public RandomVariable 
+{
+public:
+  /**
+   * Constructs a gamma random variable with alpha = 1.0 and beta = 1.0
+   */
+  GammaVariable ();
+
+  /**
+   * \param alpha alpha parameter of the gamma distribution
+   * \param beta beta parameter of the gamma distribution
+   */
+  GammaVariable (double alpha, double beta);
+
+  /**
+   * \brief call RandomVariable::GetValue
+   * \return A floating point random value 
+   *
+   * Note: we have to re-implement this method here because the method is
+   * overloaded below for the two-argument variant and the c++ name resolution
+   * rules don't work well with overloads split between parent and child
+   * classes.
+   */
+  double GetValue (void) const;
+  
+  /**
+   * \brief Returns a gamma random distributed double with parameters alpha and beta.
+   * \param alpha alpha parameter of the gamma distribution
+   * \param beta beta parameter of the gamma distribution
+   * \return A floating point random value
+   */
+  double GetValue(double alpha, double beta) const;
+};
+
+/**
+ * \brief Erlang Distributed Random Variable
+ * \ingroup randomvariable
+ *
+ * ErlangVariable defines a random variable with Erlang distribution.
+ *
+ * The Erlang distribution is a special case of the Gamma distribution where k
+ * (= alpha) is a non-negative integer. Erlang distributed variables can be
+ * generated using a much faster algorithm than gamma variables.
+ *
+ * This class supports the creation of objects that return random numbers from
+ * a fixed Erlang distribution. It also supports the generation of single
+ * random numbers from various Erlang distributions.
+ *
+ * The probability density function is defined over the interval [0,+inf) as:
+ * \f$ \frac{x^{k-1} e^{-\frac{x}{\lambda}}}{\lambda^k (k-1)!}\f$
+ * where \f$ mean = k \lambda \f$ and 
+ * \f$ variance = k \lambda^2\f$
+ */
+class ErlangVariable : public RandomVariable 
+{
+public:
+  /**
+   * Constructs an Erlang random variable with k = 1 and lambda = 1.0
+   */
+  ErlangVariable ();
+
+  /**
+   * \param k k parameter of the Erlang distribution. Must be a non-negative integer.
+   * \param lambda lambda parameter of the Erlang distribution
+   */
+  ErlangVariable (unsigned int k, double lambda);
+
+  /**
+   * \brief call RandomVariable::GetValue
+   * \return A floating point random value 
+   *
+   * Note: we have to re-implement this method here because the method is
+   * overloaded below for the two-argument variant and the c++ name resolution
+   * rules don't work well with overloads split between parent and child
+   * classes.
+   */
+  double GetValue (void) const;
+  
+  /**
+   * \brief Returns an Erlang random distributed double with parameters k and lambda.
+   * \param k k parameter of the Erlang distribution. Must be a non-negative integer.
+   * \param lambda lambda parameter of the Erlang distribution
+   * \return A floating point random value
+   */
+  double GetValue(unsigned int k, double lambda) const;
+};
+
+/**
  * \brief Triangularly Distributed random var
  * \ingroup randomvariable
  *