#ifndef ANDERSEN_THERMOSTAT_H
#define ANDERSEN_THERMOSTAT_H

#ifndef THERMOSTAT_H
#include "Thermostat.h"
#endif

#ifndef VECTOR_3D_IMPL_H
#include "Vector_3D_impl.h"
#endif

#ifndef HYDROGEN_ATOMS_H
#include "Hydrogen_atoms.h"
#endif

//#include "__rand.h"

#ifdef WIN32
#pragma warning(disable : 4800)
#pragma warning(disable : 4996)
#endif
#include "MersenneTwister.h"

#include "Flaw.h"

namespace MM
{

class Andersen_thermostat : public Thermostat
{
    double          nu_;
    //T_Base_random   rand_;
    MTRand          distr_;

public:
    Andersen_thermostat (Model & model, Atom_group & atom_group, MD & simulation)
    :   
        Thermostat (model, atom_group, simulation), 
            nu_(0.1), 
            //nu_(0.25), 

            //rand_(0., 1.),
            distr_(1966)
    {
        //rand_.set_method (T_Base_random::Mix);
    }

    char const* title () const      {return "Andersen";}

    void            set_nu (double nu) {nu_ = nu;}

    virtual void    handle_execute ()
    {
        MD_simulation_.algorithm().description ("Andersen thermostat");
        //Hydrogen_atoms group (atom_group_);
        Atom_group & group = atom_group_;

        double dt    = MD_simulation_.step_length ();                           // time step of integration algorithm
        double T0    = T0_ * System_of_Units::singleton().gas_constant();
        int    count = group.atom_count();

        for (int i=0;  i<count;  ++i)
        {
            //if (rand_.next() < nu_*dt)
            if (distr_.rand() < nu_*dt)
            {
                Atom & atom  = group.atom (i);
                double mass  = atom.mass();
                double sigma = sqrt (T0 / mass);
                double x     = sigma * gauss ();
                double y     = sigma * gauss ();
                double z     = sigma * gauss ();
                atom.kit().set_velocity (Vector_3D_impl (x, y, z));
            }
        }
    }

    double gauss ()
    {
        return distr_.randNorm ();

        //double R_2 = 2., v1, v2;

        //while (R_2 >= 1.)
        //{
        //    v1 = 2. * rand_.next() - 1.;
        //    v2 = 2. * rand_.next() - 1.;
        //    R_2 = v1*v1 + v2*v2;
        //}

        //double result = v1 * sqrt (-2. * log(R_2) / R_2);
        //return result;
    }

    double              response_constant () const       {return nu_;}
    void                set_response_constant (double v) {nu_ = v;}
};

}//MM

#endif //ANDERSEN_THERMOSTAT_H
