#include "Module.h"

#include "Peep.h"

#include "Path.h"
#include "Project.h"
#include "Model.h"
#include "Model_kit.h"
#include "Atom_kit.h"

#include "Abstracts.h"
#include "Create.h"
#include "Panel_profiler.h"

#include "System_of_Units.h"
#include "Temperature.h"
#include "Pressure.h"
#include "Charge_unit.h"
#include "Dipole_moment_unit.h"
#include "Selected_atoms.h"
#include "Defs.h"
#include "anint.h"

#include <math.h>
#include <fstream>

namespace MM
{

class Module_analyse : public Module
{
    static Module_analyse   singleton_;
                            Module_analyse() ;
    void                    init ();

    void                    calc_RDF        ();
    void                    calc_temperature();
    void                    calc_pressure   ();
    void                    calc_dipole     ();
    void                    profiler        ();
};
Module_analyse  Module_analyse::singleton_;


Module_analyse::Module_analyse () : Module ("Analyse") { }

void  Module_analyse::
init ()
{
    adopt_command (prototype <Show_information>().clone())
        //.set_icon ("Analyse", "info.png")
        .set_menu ("&Analyse", "&Information");

    adopt_command (create_command (*this, &Module_analyse::calc_RDF))
        .set_menu ("&Analyse", "RDF");

    adopt_command (create_command (*this, &Module_analyse::calc_temperature))
        .set_menu        ("&Analyse", "&Temperature");
    
    adopt_command (create_command (*this, &Module_analyse::calc_pressure))
        .set_menu        ("&Analyse", "Pressure");
    
    adopt_command (create_command (*this, &Module_analyse::calc_dipole))
        .set_menu        ("&Analyse", "Dipole moment");
    
    adopt_command (create_command (*this, &Module_analyse::profiler))
        .set_menu        ("&Analyse", "Profiler");
}

void Module_analyse::
calc_RDF ()
{
    Model & model = Project::singleton().current_model();

    if (model.kit().boundary_conditions ().type_name() != "Box")
    {
        to_user().fatal_error ("Periodic box should be defined.");
        return;
    }

    Selected_atoms selected (model);
    int fine_type_1, fine_type_2, i, j, k;

    if (selected.atom_count() == 1)
    {
        fine_type_1 = fine_type_2 = selected.atom (0).kit().fine_type_id();
    }
    else if (selected.atom_count() == 2)
    {
        fine_type_1 = selected.atom (0).kit().fine_type_id();
        fine_type_2 = selected.atom (1).kit().fine_type_id();
    }
    else 
    {
        to_user().fatal_error ("To calculate RDF select one or two appropriate atoms.");
        return;
    }

    int bin [100];
    for (i=0;  i<100;  ++i)
        bin[i] = 0;

    double L = model.kit().boundary_conditions ().box ().a_mod_rectangle();//fix
    double lIL = 1. / L;
    int count = model.atom_count(), N=0;

    for (i=0;  i<count;  ++i)
    {
        Atom const& atom_1 = model.atom (i);

        if (atom_1.kit().fine_type_id() != fine_type_1)
            continue;

        for (j=0;  j<count;  ++j)
        {
            Atom const& atom_2 = model.atom (j);

            if (atom_2.kit().fine_type_id() != fine_type_2 || i==j)
                continue;

            double x = atom_1.x() - atom_2.x();
            double y = atom_1.y() - atom_2.y();
            double z = atom_1.z() - atom_2.z();

            x -= L * anint (x * lIL);
            y -= L * anint (y * lIL);
            z -= L * anint (z * lIL);

            double r = sqrt (x*x + y*y + z*z);
            
            k = int(r*10);
            if (k < 100)
                ++bin [k];

            ++N;
        }
    }
    N = sqrt(double(N));

    std::ofstream out ((Path::report() + "RDF.xls").c_str());
    
    double V = model.kit().boundary_conditions ().box ().volume();
    double K = 4. * __Pi__ * (N / V ) / 3.;

    for (k=0;  k<100;  ++k)
    {
        double r_lower = k/10.;
        double r_upper = r_lower + 0.1;

        out << (r_lower+0.05) << "\t" 
            << (double(bin[k]) / N) / 
                (K * (r_upper*r_upper*r_upper - r_lower*r_lower*r_lower)) 
            << std::endl;
    }

    to_user().status ("RDF saved in Report/RDF.xls");
}

void Module_analyse::
calc_temperature ()
{
    Model       & model       = Project::singleton().current_model();
    Temperature temperature (model);
    Temperature_unit value;
    System_of_Units::singleton().set_temperature (value, temperature.value());
    to_user(). info (Text("temperature ") + value.K() + " K    " + value.deg() + " degrees");
}

void Module_analyse::
calc_pressure ()
{
    Model       & model       = Project::singleton().current_model();
    Interaction & interaction = model.kit().interaction();
    model.kit().strip (true);
    interaction.add_force ();
    Pressure pressure (model);
    Pressure_unit value;
    System_of_Units::singleton().set_pressure (value, pressure.value());
    to_user(). info (Text("pressure ") + value.atm() + " atm    " + value.Pa() + " Pa");
}

void Module_analyse::
calc_dipole ()
{
    double x=0, y=0, z=0;

    Model       & model       = Project::singleton().current_model();
    for (int i=0;  i<model.atom_count();  ++i)
    {
        Atom const& atom = model.atom (i);
        double q = atom.charge();
        x += q * atom.x();
        y += q * atom.y();
        z += q * atom.z();
    }
    Charge_unit charge_unit;
    charge_unit.set_electron(1);
    double C = charge_unit.C();
    Dipole_moment_unit dipole_moment_unit;
    dipole_moment_unit.set_Cm (1);
    double D = dipole_moment_unit.D();
    x *= C * D * 1e-10;
    y *= C * D * 1e-10;
    z *= C * D * 1e-10;
    double result = sqrt (x*x + y*y + z*z);

    to_user(). info (Text("Dipole_moment ") + result + " Debye" +
    "  X " + x + "  Y " + y + "  Z " + z);
}

void Module_analyse::
profiler ()
{
    Panel_profiler::singleton()->popup ();
}

}//MM
