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

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

#ifndef _RANDOMGENERATOR_H_
#define _RANDOMGENERATOR_H_

#undef M_PI
#define M_PI 3.141592653589793238462643

class rnd_uniform_distribution
{

  private:
    int a;
    int c;
    unsigned int m;
    unsigned _seed;
    double a1;
    double a2;

    int drnd() { return( _seed = ( a * _seed + c ) % m ); }

  public:

    /*constructors, destructors: */
    rnd_uniform_distribution() : _seed( 0 ), a( 1103515245 ), c( 12345 ), m( 2147483648 ), a1(0.0), a2(1.0) {}
    rnd_uniform_distribution(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 generator() { return (a2-a1)*(double) drnd()/ m + a1; }

};

class rnd_normal_distribution
{

  private:
    unsigned _seed;
    double mean;
    double sdev;

  public:
    rnd_normal_distribution() : _seed( 0 ), mean( 0), sdev(1.0) {}
    rnd_normal_distribution(double m,double s) : _seed( 0 ), mean( m ), sdev(s) {}
    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_