#include "Collect_molecular_mass_center_temperature.h"

#include "Model.h"
#include "Molecule.h"
#include "Mass_center.h"
#include "Path.h"
#include "Point_3D_impl.h"
#include "System_of_Units.h"

namespace MM
{

Collect_molecular_mass_center_temperature::
Collect_molecular_mass_center_temperature (Atom_group & atom_group, Text in)
:
    Collect_scalar (atom_group, in)
{
//    new_graph_presentation ("temperature");
    new_file_presentation (Path::report() +"dyn_Temperature_CM.xls");
}


Monitor * Collect_molecular_mass_center_temperature::
clone (Atom_group & atom_group, Text const & args, Text in) const
{
    return new Collect_molecular_mass_center_temperature (atom_group, in);
}

double Collect_molecular_mass_center_temperature::
extract_current_value ()
{
    Model & model = dynamic_cast<Model &> (atom_group());
    int     count = model.molecule_count();
    double  sum = 0.;

    const Mass_center mc (const_cast<Atom_group &>(atom_group()));
    double mass;
    Point_3D_impl   position(0,0,0);
    Vector_3D_impl  velocity(0,0,0);
    mc.evaluate (position, &mass, &velocity);

    double velocity_2 = velocity.length_2();
    double K_mc = mass * velocity_2;

    for (int i=0;  i<count;  ++i)
    {
        Molecule & molecule = model.molecule (i);
        Mass_center mc (molecule);
        double mass;
        Point_3D_impl   position(0,0,0);
        Vector_3D_impl  velocity(0,0,0);
        mc.evaluate (position, &mass, &velocity);

        double velocity_2 = velocity.length_2();
        sum += mass * velocity_2;
    }

    sum -= K_mc;
    double kB = System_of_Units::singleton().gas_constant();
    return  sum / ((3. * count - 3.) * kB);
}

}//MM

