source: issm/trunk-jpl/src/c/shared/Random/randomgenerator.h@ 26599

/*!\file: randomgenerator.h
 * \brief prototypes for randomgenerator.h
 */

#ifndef _RANDOMGENERATOR_H_
#define _RANDOMGENERATOR_H_

/*Headers:*/
/*{{{*/
#include <cmath>
/*}}}*/

#undef M_PI
#define M_PI 3.141592653589793238462643

class rnd_uniform_distribution
{

  private:
    int a;
    int c;
    unsigned int m;
    unsigned _seed;
    double lbound;
    double ubound;

  public:

    /*constructors, destructors: */
    rnd_uniform_distribution();
    rnd_uniform_distribution(double a_1,double a_2) : _seed( 0 ), a( 1103515245 ), c( 12345 ), m( 2147483648 ), lbound(a_1), ubound(a_2) {}
    ~rnd_uniform_distribution(){}

    void seed( unsigned int s ) { _seed = s; }
    unsigned int get_seed() { return _seed; }
    double generator(){
      _seed = ( a * _seed + c ) % m;
      return (ubound-lbound)*(double) _seed/ m + lbound;
    }

};

class rnd_normal_distribution
{

  private:
    unsigned _seed;
    double mean;
    double sdev;

  public:

    /*constructors, destructors: */
    rnd_normal_distribution() : _seed( 0 ), mean( 0), sdev(1.0) {}
    rnd_normal_distribution(double m,double s) : _seed( 0 ), mean( m ), sdev(s) {}
    ~rnd_normal_distribution(){}

    void seed( unsigned int s ) { _seed = s; }
    double generator()
    {
      rnd_uniform_distribution unifdistri;
      unifdistri.seed(_seed);

      double u1 = unifdistri.generator();
      double u2 = unifdistri.generator();

      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_