#ifndef MaxwellInflowGen_IS_INCLUDED
#define MaxwellInflowGen_IS_INCLUDED

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

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

using namespace std;


class MaxwellInflowGen {

 private:
  Generator* generator;
  unsigned int type;

  Rnd rnd1, rnd2;
  double c;

  
  //Physicals
  double rho, T;        //Density, temperature
  double V[3];          //Fluid velocity
  double phi, theta;    //Inwards normal vector orientation (rotation, azimuth)
  double normal[3];     //Inwards normal vector
  double orthoBasis[9]; //orthogonal matrix
  double a;
  double prefactor;     //Normalizing prefactor for the distribution

 public:

  //Constructor
  MaxwellInflowGen(double _rho,
		   double _T,
		   double _V[],
		   double _phi,
		   double _theta,
		   unsigned long int s1,
		   unsigned long int s2,
		   Generator* _generator ) :
    rho(_rho), T(_T), phi(_phi), theta(_theta), rnd1(s1), rnd2(s2), generator(_generator) {

    type = _generator->isType();

    V[0] = _V[0];
    V[1] = _V[1];
    V[2] = _V[2];

    c = sqrt(2*T);

    int i = 0;
    orthoBasis[i++] = cos(phi)*cos(theta);
    orthoBasis[i++] = -sin(phi);
    orthoBasis[i++] = cos(phi)*sin(theta);
    orthoBasis[i++] = sin(phi)*cos(theta);
    orthoBasis[i++] = cos(phi);
    orthoBasis[i++] = sin(phi)*sin(theta);
    orthoBasis[i++] = -sin(theta);
    orthoBasis[i++] = 0;
    orthoBasis[i++] = cos(theta);

    normal[0] = cos(phi)*sin(theta);
    normal[1] = sin(phi)*sin(theta);
    normal[2] = cos(theta);

    a = (normal[0]*V[0] +
	 normal[1]*V[1] +
	 normal[2]*V[2]  )/c;

    generator->set(a);

    prefactor = 1./(2*M_PI*(exp(-a*a) + a*M_SQRTPI*(1 + gsl_sf_erf(a)))*T*T);
    /* gsl-specific content:
     * const double M_PI; (pi)
     * const double M_SQRTPI; (square root of pi) 
     * double gsl_sf_erf(double); (the error function)
     */

  } // end constructor

  
  void generateVelocity(double result[]) {
    double tmp[3];

    //0-direction velocity distributed normally with mean 0, variance 0.5
    tmp[0] = rnd1.getGaussian(M_SQRT1_2);
    /* gsl-specific content:
     * const double M_SQRT1_2; (square root of 0.5)
     */

    //1-direction velocity distributed normally with mean 0, variance 0.5
    tmp[1] = rnd2.getGaussian(M_SQRT1_2);
    /* gsl-specific content:
     * const double M_SQRT1_2; (square root of 0.5)
     */

    //2-direction velocity
    switch (type) {
    case envelope1:
      tmp[2] = - ((Envelope1*) generator)->generate();
      break;
    case envelope2:
      tmp[2] = - ((Envelope2*) generator)->generate();
      break;
    case envelope3:
      tmp[2] = - ((Envelope3*) generator)->generate();
      break;
    case envelope4:
      tmp[2] = - ((Envelope4*) generator)->generate();
      break;
    case inversionMethod:
      tmp[2] = - ((InversionMethod*) generator)->generate();
      break;
    case boxEnvelope:
      tmp[2] = - ((BoxEnvelope*) generator)->generate();
      break;
    case reservoirEnvelope:
      tmp[2] = - ((ReservoirEnvelope*) generator)->generate();
      break;
    }

    //Transform all 3 velocity components back into global frame
    result[0] = V[0] + c*(orthoBasis[0]*tmp[0] +
			  orthoBasis[1]*tmp[1] +
			  orthoBasis[2]*tmp[2]  );
    result[1] = V[1] + c*(orthoBasis[3]*tmp[0] +
			  orthoBasis[4]*tmp[1] +
			  orthoBasis[5]*tmp[2]  );
    result[2] = V[2] + c*(orthoBasis[6]*tmp[0] +
			  orthoBasis[7]*tmp[1] +
			  orthoBasis[8]*tmp[2]  );

  } // end generateVelocity()

  
  /* The following function uses the stl (standard template library) class
   * "string", and returns a string-object containing essential information
   * stored in its MaxwellInflowGen-object. The whole function is not
   * essential for the generator and may be commented out entirely.
   */
  string summary() {
      string str;
      char   buffer[100];

      snprintf(buffer, 100, "Normal vector:\n\t[%10.5f %10.5f %10.5f]\n\n",
	       normal[0], normal[1], normal[2]);
      str += buffer;

      snprintf(buffer, 100, "Orthogonal matrix:\n");
      str += buffer;
      snprintf(buffer, 100, "\t[%10.5f %10.5f %10.5f\n",
	       orthoBasis[0], orthoBasis[1], orthoBasis[2]);
      str += buffer;
      snprintf(buffer, 100, "\t %10.5f %10.5f %10.5f\n",
	       orthoBasis[3], orthoBasis[4], orthoBasis[5]);
      str += buffer;
      snprintf(buffer, 100, "\t %10.5f %10.5f %10.5f]\n\n",
	       orthoBasis[6], orthoBasis[7], orthoBasis[8]);
      str += buffer;

      snprintf(buffer, 100, "Prefactor:\n\t%10.5f\nInverse:\n\t%10.5f\n",
	       prefactor, 1./prefactor);
      str += buffer;

      return str;

  } // end summary()


  double   get_rho()        {return rho;}
  double   get_T()          {return T;}
  double*  get_V()          {return V;}
  double   get_phi()        {return phi;}
  double   get_theta()      {return theta;}
  double*  get_normal()     {return normal;}
  double*  get_orthoBasis() {return orthoBasis;}
  double   get_a()          {return a;}
  double   get_prefactor()  {return prefactor;}
  
}; // end class MaxwellInflowGen

#endif