#ifndef generators_IS_INCLUDED
#define generators_IS_INCLUDED

#include <cmath>
#include <iostream>
#include "Rnd.h"

//Included from Gnu Scientific Library:
#include <gsl/gsl_math.h>   //Provides mathematical constants and elemetary
			    //mathematical functions.
#include <gsl/gsl_sf_erf.h> //Provides evaluation of the error function.

using namespace std;


/* List of constructor prototypes and short descriptions of parameters:
 *
 * ----------------------------------------------------------------------
 * Envelope1(unsigned long int s1, unsigned long int s2);
 *
 *    (two seeds to random number generators)
 *
 * ----------------------------------------------------------------------
 * Envelope2(unsigned long int s1, unsigned long int s2
 *           unsigned long int s3);
 *
 *    (three seeds to random number generators)
 *
 * ----------------------------------------------------------------------
 * Envelope3(unsigned long int s1, unsigned long int s2
 *           unsigned long int s3, unsigned long int s4);
 *
 *    (four seeds to random number generators)
 *
 * ----------------------------------------------------------------------
 * Envelope4(unsigned long int s1, unsigned long int s2
 *           unsigned long int s3, unsigned long int s4);
 *
 *    (four seeds to random number generators)
 *
 * ----------------------------------------------------------------------
 * InversionMethod(double _error, unsigned long int s1);
 *
 *    (_error: absolute error for the Newton iteration,
 *     s1:     seed to the random number generator)
 *
 * ----------------------------------------------------------------------
 * BoxEnvelope(double _zLTa, double _zGTa,
 *             unsigned long int s1, unsigned long int s2 );
 *
 *    (_zLTa:  should be set to min(a - alpha1, -alpha2)
 *     _zGTa:  should be set to min(a, alpha2)
 *             where alpha1 = 1 or 2,
 *                   alpha2 = 3,
 *                   a      = V.e/sqrt(2T) (fractional fluid velocity
 *                                          along surface normal vector)
 *     s1, s2: two seeds to random number generators)
 *
 * ----------------------------------------------------------------------
 * ReservoirEnvelope(double _zLTa,
 *                   unsigned long int s1, unsigned long int s2 );
 *
 *    (_zLTa:  should be set to min(a - 1, -alphaR)
 *             where alphaR = 3
 *                   a      = V.e/sqrt(2T) (fractional fluid velocity
 *                                          along surface normal vector) 
 *     s1, s2: two seeds to random number generators)
 *
 */


// List of generator type identyfiers
const unsigned int envelope1         = 0x0;
const unsigned int envelope2         = 0x1;
const unsigned int envelope3         = 0x2;
const unsigned int envelope4         = 0x3;
const unsigned int inversionMethod   = 0x4;
const unsigned int boxEnvelope       = 0x5;
const unsigned int reservoirEnvelope = 0x6;


// Base class for generator structs
struct Generator {
  double a;
  
  virtual void set(double _a) {};
  virtual unsigned int isType() {return 0;}
  double generate() {return 0;}
};


struct Envelope1 : public Generator {

  //Unique type specifier
  unsigned int isType() {return envelope1;}
  
  //Generator-specific variables
  double aa;
  Rnd rnd1, rnd2;

  
  //Constructor
  Envelope1(unsigned long int s1, unsigned long int s2) : rnd1(s1), rnd2(s2) {}

  
  void set(double _a) {
    a  = _a;
    aa = _a*_a;
  }

  
  //The generator
  inline double generate() {
    double result;
    while(true) {
      result = -sqrt(aa - log(rnd1.getUniform()));
      if ((result-a)/result > rnd2.getUniform())
	return result;
    }
  }

}; // end struct Envelope1


struct Envelope2 : public Generator {

  //Unique type specifier
  unsigned int isType() {return envelope2;}

  //Generator-specific variables
  double z_a, beta_a, z_a_2, beta_a_2, aMz_aINV, aMbeta_a, c;
  Rnd rnd1, rnd2, rnd3;

  
  //Constructor
  Envelope2(unsigned long int s1,
	    unsigned long int s2,
	    unsigned long int s3) :
    rnd1(s1), rnd2(s2), rnd3(s3) {
  }


  void set(double _a) {
    a            = _a;
    z_a          = 0.5*(a - sqrt(a*a + 2));
    beta_a       = a - (1 - a)*(a - z_a);
    z_a_2        = z_a*z_a;
    beta_a_2     = beta_a*beta_a;
    double aMz_a = a - z_a;
    aMz_aINV     = 1./aMz_a;
    aMbeta_a     = a - beta_a;
    c            = exp(-beta_a_2)/(exp(-beta_a_2) + 2*aMz_a*aMbeta_a*exp(-z_a_2));
  }

  
  //The generator
  inline double generate() {
    double result;

    while(true) {
      if (c > rnd1.getUniform()) {
	result = -sqrt(beta_a_2 - log(rnd2.getUniform()));
	if ((result-a)/result > rnd3.getUniform())
	  return result;
      }
      else {
	result = beta_a + aMbeta_a*rnd2.getUniform();
	if ((a-result)*aMz_aINV*exp(z_a_2 - result*result) > rnd3.getUniform())
	  return result;
      }
    }

  }

}; // end struct Envelope2


struct Envelope3 : public Generator {

  //Unique type specifier
  unsigned int isType() {return envelope3;}
  
  //Generator-specific variables
  double c1, c2;
  Rnd rnd1, rnd2, rnd3, rnd4;

  
  //Constructor
  Envelope3(unsigned long int s1,
	    unsigned long int s2,
	    unsigned long int s3,
	    unsigned long int s4 ) :
    rnd1(s1), rnd2(s2), rnd3(s3), rnd4(s4) {}

  
  void set(double _a) {
    a  = _a;
    c1 = a*M_SQRTPI/(a*M_SQRTPI + 1 + a*a);
    c2 = (a*M_SQRTPI + 1)/(a*M_SQRTPI + 1 + a*a);
    /* gsl-specific content:
     * const double M_SQRTPI; (square root of pi)
     */
  }

  
  //The generator
  inline double generate() {
    double res;
    while(true) {
      res = rnd1.getUniform();
      if      (c1 > res)
	return -fabs(rnd2.getGaussian(M_SQRT1_2));
        /* gsl-specific content:
	 * const double M_SQRT1_2; (square root of 0.5)
	 */

      else if (c2 > res)
	return -sqrt(-log(rnd3.getUniform()));
      else {
	res = (1 - sqrt(rnd3.getUniform()))*a;
	if (exp(-res*res) > rnd4.getUniform())
	  return res;
      }
    }
  }

}; // end struct Envelope3


struct Envelope4 : public Generator {

  //Unique type specifier
  unsigned int isType() {return envelope4;}
  
  //Generator-specific variables
  double aINV, c;
  Rnd rnd1, rnd2, rnd3, rnd4;

  
  //Constructor
  Envelope4(unsigned long int s1,
	    unsigned long int s2,
	    unsigned long int s3,
	    unsigned long int s4 ) :
    rnd1(s1), rnd2(s2), rnd3(s3), rnd4(s4) {}


  void set(double _a) {
    a    = _a;
    aINV = 1./a;
    c    = 1/(2*a*M_SQRTPI + 1);
    /* gsl-specific content:
     * const double M_SQRTPI; (square root of pi)
     */
  }

  
  //The generator
  inline double generate() {
    double result;
    while(true) {
      if (c > rnd1.getUniform())
	result =  -sqrt(-log(rnd2.getUniform()));
      else
	result = rnd3.getGaussian(M_SQRT1_2);
        /* gsl-specific content:
	 * const double M_SQRT1_2; (square root of 0.5)
	 */

      if ((a - result)*aINV > rnd4.getUniform())
	return result;
    }
  }

}; // end struct Envelope4


struct InversionMethod : public Generator {

  //Unique type specifier
  unsigned int isType() {return inversionMethod;}
  
  //Generator-specific variables
  double error, m_a, m_aINV, c, z0;
  Rnd rnd1;

  
  //Constructor
  InversionMethod(double _error, unsigned long int s1) :
    error(_error), rnd1(s1) {}


  void set(double _a) {
    a = _a;
    m_a    = exp(-a*a) + a*M_SQRTPI*(1 + gsl_sf_erf(a));
    m_aINV = 1./m_a;
    c      = a*M_SQRTPI;
    z0     = 0.5*(a - sqrt(a*a +2));
    /* gsl-specific content:
     * const double M_SQRTPI; (square root of pi)
     */
  }

  
  //The generator
  inline double generate() {
    double u = rnd1.getUniform();

    //Initial guess
    double z = z0;
    double errorCurrent;

    while ( fabs(errorCurrent = (exp(-z*z) + c*(1 + gsl_sf_erf(z)))*m_aINV - u) >= error)
      z -= errorCurrent*m_a/(2*(a - z)*exp(-z*z));
    /* gsl-specific content:
     * double gsl_sf_erf(double); (the error function)
     */

    return z;

  }

}; // end struct InversionMethod


struct BoxEnvelope : public Generator {

  //Unique type specifier
  unsigned int isType() {return boxEnvelope;}
  
  //Generator-specific variables
  double zLTa, zGTa, aMz_aINV, z_a_2, c;
  Rnd rnd1, rnd2;

  
  //Constructor
  BoxEnvelope(double _zLTa, double _zGTa,
	      unsigned long int s1,
	      unsigned long int s2 ) :
    zLTa(_zLTa), zGTa(_zGTa), rnd1(s1), rnd2(s2) {}


  void set(double _a) {
    double z_a  = 0.5*(a - sqrt(a*a +2));
    aMz_aINV    = 1./(a - z_a);
    z_a_2       = z_a*z_a;
    c           = zGTa - zLTa;
  }

  
  //The generator
  inline double generate() {
    double z;

    do
      z = zLTa + c*rnd1.getUniform();
    while ((a - z)*aMz_aINV*exp(z_a_2 - z*z) <= rnd2.getUniform());

    return z;

  }

}; // end struct BoxEnvelope4


struct ReservoirEnvelope : public Generator {

  //Unique type specifier
  unsigned int isType() {return reservoirEnvelope;}
  
  //Generator-specific variables
  double zLTa, c;
  Rnd rnd1, rnd2;

  
  //Constructor
  ReservoirEnvelope(double _zLTa,
		    unsigned long int s1,
		    unsigned long int s2 ) :
    zLTa(_zLTa), rnd1(s1), rnd2(s2) {}


  void set(double _a) {
    a = _a;
    c = 1./(a - zLTa);
  }

  
  //The generator
  inline double generate() {
    double result;
    
    do
      if (a <= 0)
	result = -fabs(rnd1.getGaussian(M_SQRT1_2));
      else
	result = rnd1.getGaussian(M_SQRT1_2);
    while ( (a-result)*c <= rnd2.getUniform() );
    /* gsl-specific content:
     * const double M_SQRT1_2; (square root of 0.5)
     */

    return result;

  }

}; // end struct ReservoirEnvelope

#endif