#include "Measure_geometry_impl.h"

#include "Log.h"

#include "Flaw.h"
#include "Dihedral_angle.h"
#include "Angle.h"
#include "Mass_center.h"
#include "Project.h"
#include "Atom.h"
#include "Atom_kit.h"
#include "Molecule.h"
#include "Vector_3D_impl.h"
#include "Point_3D_impl.h"
#include "Branch.h"
#include "Default_builder.h"

#include "Mach_eps.h"
#include <math.h>

#include <iostream>

#ifdef WIN32
#pragma warning(disable : 4800)
#endif

namespace MM
{
Measure_geometry_impl::
Measure_geometry_impl ()
:
//    last_point (0,0,0)
    last_movement_(1,0,0)
{
    last_length_.        set_current_type ("A");
    last_angle_.         set_current_type ("deg");
    last_dihedral_angle_.set_current_type ("deg");

    four_points_.push_back (Point(0, "Atom", true));
    four_points_.push_back (Point(1, "Atom", true));
    four_points_.push_back (Point(2, "Atom", false));
    four_points_.push_back (Point(3, "Atom", false));
}

void Measure_geometry_impl::
prepend_point ()
{
    INVARIANT;

    int number = four_points_[0].number() - 1;
    number = number < 0 ? 0 : number;

    four_points_.pop_back ();
    four_points_.push_front (Point(number, "Atom", false));

    INVARIANT;
}

void Measure_geometry_impl::
append_point ()
{
    INVARIANT;

    //bool is_active = four_points_[0].is_active();
    int atom_count = Project::singleton().current_model().atom_count();
    int number = four_points_[3].number() + 1;
    number = number >= atom_count ? atom_count-1 : number;

    four_points_.pop_front ();
    four_points_.push_back (Point(number, "Atom", false));
    //four_points_[3].set_active (is_active);

    INVARIANT;
}

void Measure_geometry_impl::
set_atom_type (int point, int atom_number)
{
    REQUIRE ("Four points.", point >=0 && point <4);

    four_points_[point].set_number (atom_number);
    four_points_[point].set_space_position 
        (Project::singleton ().current_model().atom (atom_number).position());
    four_points_[point].set_type ("Atom");
    four_points_[point].set_active (true);
}

void Measure_geometry_impl::
set_space_position_type (int point, Point_3D const &space_position)
{
    REQUIRE ("Four points.", point >=0 && point <4);

    four_points_[point].set_number (-1);
    four_points_[point].set_space_position (space_position);
    four_points_[point].set_type ("Space position");
    four_points_[point].set_active (true);
}

void Measure_geometry_impl::
set_mass_center_type (int point, int mol_number) 
{
    REQUIRE ("Four points.", point >=0 && point <4);

    four_points_[point].set_number (mol_number);
    four_points_[point].set_type ("Mass center");
    four_points_[point].set_active (true);

    if (is_active (point))
        calc_space_position (four_points_[point]);
}

void Measure_geometry_impl::
get_space_position (int point, Point_3D & result) const
{
    result.affine (four_points_[point].space_position(), 1);
}

void Measure_geometry_impl::
set_space_position (int pnt, Point_3D const& position)
{
    Point & point = four_points_[pnt];

    //if (!point.is_active())
    if (!is_active(pnt))
        return;

    if (point.type() == "Atom")
    {
        Project::singleton().current_model().atom(point.number()).
            set_position (position);
    }

    else if (point.type() == "Mass center")
    {
        Mass_center mc (Project::singleton().current_model().
            molecule(point.number()));

        mc. set (position);
    }

    else if (point.type () == "Space position")
    {
        point.set_space_position (position);
    }

    else 
        FLAW ("Unknown type.");
}

int Measure_geometry_impl::
number (int point) const
{
    return four_points_[point].number();
}

void Measure_geometry_impl::
set_number (int point, int number)
{
    four_points_[point].set_number (number); 
}

const char * Measure_geometry_impl::
element (Point const* point) const
{
    if (point->type () == "Atom")
        return Project::singleton().current_model().
            atom (point->number()).element().c_str();
    
    else if (point->type () == "Mass center")
        return "mc";

    else if (point->type () == "Space position")
        return "sp";

    FLAW ("Unknown type.");
    return 0;
}

const char * Measure_geometry_impl::
mm_type_c_str (Point const* point) const
{
    if (point->type () == "Atom")
        return Project::singleton().current_model().atom 
            (point->number()).kit().mm_type_c_str();
    
    return "x";
}

Text const & Measure_geometry_impl::
type (int point) const 
{
    REQUIRE ("Four points.", point >=0 && point <4);

    return four_points_[point].type ();
}

void Measure_geometry_impl::
set_active (int point, bool on)
{
    REQUIRE ("Four points.",            point >=0 && point <4);

    four_points_[point].set_active (on);
}

bool Measure_geometry_impl::
is_active (int point_N) const
{
    REQUIRE ("Four points.", point_N >=0 && point_N <4);

    bool result = true;
    Point const& point = four_points_[point_N];

    result &= point.is_active ();

    if (point.type () == "Atom")
        result &= (point.number() 
            < Project::singleton().current_model().atom_count());

    else if (point.type () == "Mass center")
        result &= (point.number() 
            < Project::singleton().current_model().molecule_count());
    
    else if (point.type () == "Space position")
        ;

    else 
        FLAW ("Unknown point type.");


    //fix

    //for (int k=0;  k<4;  ++k)
    //{
    //    four_points_[k].is_active();
    //    four_points_[k].number();
    //}

    if (point.type () != "Space position")
        if (
            point_N != 0 && four_points_[0].is_active() &&  //fix why point_N != 0
            point.type () == four_points_[0].type() && 
            point.number() == four_points_[0].number()
            ||
            point_N != 1 && four_points_[1].is_active() && 
            point.type () == four_points_[1].type() && 
            point.number() == four_points_[1].number()
            ||
            point_N != 2 && four_points_[2].is_active() && 
            point.type () == four_points_[2].type() && 
            point.number() == four_points_[2].number()
            ||
            point_N != 3 && four_points_[3].is_active() && 
            point.type () == four_points_[3].type() && 
            point.number() == four_points_[3].number()
           )
            result = false;

    return result;
}

int Measure_geometry_impl::
active_count ()
{
    return active_points ().size();
}

void Measure_geometry_impl::
get_info (Text & result)
{
    INVARIANT;

    calc_space_positions ();

    // sp-mc-C-H
    // x-x-CT-H
    // length  angle   dihedral
    //    A     rad      rad   
    // 1.587        
    //   a0     deg      deg   
    // 3.00 

    Array <Point const*> points = active_points ();
    int number_of_active_points = points.size();

    result = "";

    Text elements, types, measures, units_1, values_1, units_2, values_2
        , space(" "), dash("-");
    
    switch (number_of_active_points)
    {
        case 0:
            return;

        case 4:
        {
            //Point_3D_impl p0 (points[0]->space_position ());
            //Point_3D_impl p1 (points[1]->space_position ());
            //Point_3D_impl p2 (points[2]->space_position ());
            //Point_3D_impl p3 (points[3]->space_position ());

            Dihedral_angle  dihedral_angle (points[0]->space_position (),
                                            points[1]->space_position (),
                                            points[2]->space_position (),
                                            points[3]->space_position ());

            last_dihedral_angle_.set_rad (dihedral_angle.value());

            elements = dash + element       (points[3]);
            types    = dash + mm_type_c_str (points[3]);
            measures = "dihedral";
            units_1  = "  rad   ";
            values_1 = last_dihedral_angle_.rad();
            units_2  = "  deg   ";
            values_2 = last_dihedral_angle_.deg();
        }

        case 3:
        {
            Angle  angle 
                (points[number_of_active_points-3]->space_position (),
                 points[number_of_active_points-2]->space_position (),
                 points[number_of_active_points-1]->space_position ());

            last_angle_.set_rad (angle.value());

            elements = dash + element       (points[2]) + elements;
            types    = dash + mm_type_c_str (points[2]) + types;
            measures = " angle  "                       + measures;
            units_1  = "  rad   "                       + units_1;
            values_1 = last_angle_.rad() + space        + values_1;
            units_2  = "  deg   "                       + units_2;
            values_2 = last_angle_.deg() + space        + values_2;
        }

        case 2:
        {
            Length  length
                (points[number_of_active_points-2]->space_position (),
                 points[number_of_active_points-1]->space_position ());

            last_length_.set_A (length.value());

            elements = dash + element       (points[1]) + elements;
            types    = dash + mm_type_c_str (points[1]) + types;
            measures = "length  "                       + measures;
            units_1  = "   A    "                       + units_1;
            values_1 = last_length_.A() + space         + values_1;
            units_2  = " bohr   "                       + units_2;
            values_2 = last_length_.bohr() + space      + values_2;
        }

        case 1:
        {
            elements = element       (points[0]) + elements;
            types    = mm_type_c_str (points[0]) + types;
        }

        {   
            result  = elements;
            result += "    (";
            result += types;
            result += ")\n";
            result += measures;
            result += "\n";
            result += units_1;
            result += "\n";
            result += values_1;
            result += "\n";
            result += "\n";
            result += units_2;
            result += "\n";
            result += values_2;
            result += "\n";
            return;
        }

        default: 
            FLAW ("wrong number_of_active_points");
    }
}

void Measure_geometry_impl::
get_points_info (Text & result)
{
    int i;
    Array <Point const*> point = active_points ();

    for (i=0;  i<point.size();  ++i)
    {
        if (i != 0)
            result += "-";
        
        result += element (point[i]);
    }

    result += "   (";

    for (i=0;  i<point.size();  ++i)
    {
        if (i != 0)
            result += "-";
        
        result += mm_type_c_str (point[i]);
    }
    result += ")";
}

Text const & Measure_geometry_impl::
parameter_type_name () const
{
    static const Text no        = "No parameters";
    //static const Text point     = "Point";
    static const Text point     = no;
    static const Text distance  = "Distance";
    static const Text angle     = "Angle";
    static const Text dihedral  = "Dihedral angle";

    switch (active_points ().size())
    {
    case 0: return no; 
    case 1: return point; 
    case 2: return distance; 
    case 3: return angle; 
    case 4: return dihedral; 

    default: FLAW ("Active points out of range.");
    }
}

Unit const & Measure_geometry_impl::
parameter () const
{
    int active = active_points ().size();

    switch (active)
    {
    //case 0: return 0; 
    //case 1: return 0; 
    case 2: return last_length         (); 
    case 3: return last_angle          (); 
    case 4: return last_dihedral_angle (); 

    default: FLAW ("Active points out of range.");
    }
}

void Measure_geometry_impl::
set_parameter (double new_value)
{
    if (active_points().size() < 2)
        return;

    Text info;
    get_info (info);

    Unit const & old_value = parameter ();
    double delta = new_value - old_value.value();

    add_parameter (delta);
}

void Measure_geometry_impl::
set_parameter (Unit const & new_value)
{
    if (active_points().size() < 2)
        return;

    Text info;
    get_info (info);

    Unit const & old_value = parameter ();
    
    REQUIRE ("Parameters of the same type.", 
        typeid (new_value) == typeid (old_value));
    
    //REQUIRE ("Values in the same units.", 
    //    new_value.current_type () == old_value.current_type ());
    old_value.set_current_type (new_value.current_type ());

    //Text info;
    //get_info (info);
    double delta = new_value.value() - old_value.value();
    Own <Unit> parameter_delta (new_value.clone());
    parameter_delta().set_value (delta);
    add_parameter (parameter_delta ());
    
    /*#ifndef NDEBUG
    get_info (info);
    double result = parameter().value();
    double target = new_value.value();
    log () << "result " << result << "   target " << target << "\n"; 
    #endif//*/
    //ENSURE ("", fabs (parameter().value() - new_value.value()) < f_mach_eps_2);
}

void Measure_geometry_impl::
add_parameter (double delta)
{
    Text info;

    Array <Measure_geometry_impl::Point const*> point = active_points ();
    int active_count = point.size ();
    //CHECK ("", active_count > 1);
    if (active_count <= 1)
        return;

    Point const* last_point = point.back ();

    // Define moving object =================================================
    Atom_group * moving_object;
    Model & model = Project::singleton ().current_model();
    
    if      (last_point->type() == "Atom")
        moving_object = &model.atom (last_point->number()).kit().in_molecule();

    else if (last_point->type() == "Mass center")
        moving_object = &model.molecule (last_point->number());

    else if (last_point->type() == "Space position")
        return;
    
    else 
        FLAW (last_point->type() + " - unknown type.");
    // Define moving object -------------------------------------------------

    switch (active_count)
    {
        case 2: // length ===================================================
        {
            //Own <Length_unit> delta_unit (last_length_.clone());
            //delta_unit().set_value (delta);
            last_length_.set_value (delta);

            Vector_3D_impl movement (0,0,0);

            point[1]->space_position(). difference 
                (point[0]->space_position(), &movement);

            if (movement.length_2() < d_mach_eps_4)
            {
                movement.set_x (last_movement_.x());
                movement.set_y (last_movement_.y());
                movement.set_z (last_movement_.z());
            }

            movement.norm() *= last_length_.A();

            last_movement_.set_x (movement.x());
            last_movement_.set_y (movement.y());
            last_movement_.set_z (movement.z());

            if (point[0]->type() == "Atom" && 
                point[1]->type() == "Atom")
            {
                Atom & atom_0 = model.atom (point[0]->number());
                Atom & atom_1 = model.atom (point[1]->number());

                if (atom_0.bound_with (atom_1))
                {
                    Branch branch (atom_1, atom_0);

                    for (int i=0;  i<branch.atom_count();  ++i)
                    {
                        branch.atom(i).move (movement);
                    }
                    break;
                }
            }

            int count = moving_object->atom_count ();
            for (int i=0;  i<count;  ++i)
                moving_object->atom(i).move (movement);

            break;
        }
        case 3: // angle ====================================================
        {
            //Angle_unit const& delta_unit = 
            //    dynamic_cast <Angle_unit const&> (delta);
            last_angle_.set_value (delta);

            Vector_3D_impl a (0,0,0);
            Vector_3D_impl b (0,0,0);
            Vector_3D_impl n (0,0,0);

            point[0]->space_position(). difference 
                (point[1]->space_position(), &a);

            point[2]->space_position(). difference 
                (point[1]->space_position(), &b);

            a.cross (b, &n);

            if (point[1]->type() == "Atom" && 
                point[2]->type() == "Atom")
            {
                Atom & atom_1 = model.atom (point[1]->number());
                Atom & atom_2 = model.atom (point[2]->number());

                if (atom_1.bound_with (atom_2))
                {
                    Branch branch (atom_2, atom_1);

                    for (int i=0;  i<branch.atom_count();  ++i)
                    {
                        branch.atom(i).
                            rotate (n, point[1]->space_position(), 
                                    last_angle_.rad());
                    }
                    break;
                }
            }

            int count = moving_object->atom_count ();
            for (int i=0;  i<count;  ++i)
                moving_object->atom(i).
                    rotate (n, point[1]->space_position(), last_angle_.rad());

            break;
        }
        case 4: // dihedral_angle ===========================================
        {
            //Angle_unit const& delta_unit = 
            //    dynamic_cast <Angle_unit const&> (delta);
            last_dihedral_angle_.set_value (delta);

            Vector_3D_impl a (0,0,0);

            point[2]->space_position(). difference 
                (point[1]->space_position(), &a);

            if (point[1]->type() == "Atom" && 
                point[2]->type() == "Atom")
            {
                Atom & atom_1 = model.atom (point[1]->number());
                Atom & atom_2 = model.atom (point[2]->number());

                if (atom_1.bound_with (atom_2))
                {
                    Branch branch (atom_2, atom_1);

                    for (int i=0;  i<branch.atom_count();  ++i)
                    {
                        branch.atom(i).
                            rotate (a, point[1]->space_position(), 
                                    last_dihedral_angle_.rad());
                    }
                    break;
                }
            }

            int count = moving_object->atom_count ();
            for (int i=0;  i<count;  ++i)
                moving_object->atom(i).
                    rotate (a, point[1]->space_position(), last_dihedral_angle_.rad());

            break;
        }

        default: FLAW (Text(active_count) + " active points are out of range.");
    }

    get_info (info);
}

void Measure_geometry_impl::
add_parameter (Unit const & delta)
{
    Text info;

    Array <Measure_geometry_impl::Point const*> point = active_points ();
    CHECK ("", point.size () > 1);
    Point const* last_point = point.back ();

    // Define moving object =================================================
    Atom_group * moving_object;
    Model & model = Project::singleton ().current_model();
    
    if      (last_point->type() == "Atom")
        moving_object = &model.atom (last_point->number()).kit().in_molecule();

    else if (last_point->type() == "Mass center")
        moving_object = &model.molecule (last_point->number());

    else if (last_point->type() == "Space position")
        return;
    
    else 
        FLAW (last_point->type() + " - unknown type.");
    // Define moving object -------------------------------------------------

    switch (point.size())
    {
        case 2: // length ===================================================
        {
            Length_unit const& delta_unit = 
                dynamic_cast <Length_unit const&> (delta);

            Vector_3D_impl movement (0,0,0);

            point[1]->space_position(). difference 
                (point[0]->space_position(), &movement);

            if (movement.length_2() < d_mach_eps_4)
            {
                movement.set_x (last_movement_.x());
                movement.set_y (last_movement_.y());
                movement.set_z (last_movement_.z());
            }

            movement.norm() *= delta_unit.A();

            last_movement_.set_x (movement.x());
            last_movement_.set_y (movement.y());
            last_movement_.set_z (movement.z());

            if (point[0]->type() == "Atom" && 
                point[1]->type() == "Atom")
            {
                Atom & atom_0 = model.atom (point[0]->number());
                Atom & atom_1 = model.atom (point[1]->number());

                if (atom_0.bound_with (atom_1))
                {
                    Branch branch (atom_1, atom_0);

                    for (int i=0;  i<branch.atom_count();  ++i)
                    {
                        branch.atom(i).move (movement);
                    }
                    break;
                }
            }

            int count = moving_object->atom_count ();
            for (int i=0;  i<count;  ++i)
                moving_object->atom(i).move (movement);

            break;
        }
        case 3: // angle ====================================================
        {
            Angle_unit const& delta_unit = 
                dynamic_cast <Angle_unit const&> (delta);

            Vector_3D_impl a (0,0,0);
            Vector_3D_impl b (0,0,0);
            Vector_3D_impl n (0,0,0);

            point[0]->space_position(). difference 
                (point[1]->space_position(), &a);

            point[2]->space_position(). difference 
                (point[1]->space_position(), &b);

            a.cross (b, &n);

            if (point[1]->type() == "Atom" && 
                point[2]->type() == "Atom")
            {
                Atom & atom_1 = model.atom (point[1]->number());
                Atom & atom_2 = model.atom (point[2]->number());

                if (atom_1.bound_with (atom_2))
                {
                    Branch branch (atom_2, atom_1);

                    for (int i=0;  i<branch.atom_count();  ++i)
                    {
                        branch.atom(i).
                            rotate (n, point[1]->space_position(), 
                                    delta_unit.rad());
                    }
                    break;
                }
            }

            int count = moving_object->atom_count ();
            for (int i=0;  i<count;  ++i)
                moving_object->atom(i).
                    rotate (n, point[1]->space_position(), delta_unit.rad());

            break;
        }
        case 4: // dihedral_angle ===========================================
        {
            Angle_unit const& delta_unit = 
                dynamic_cast <Angle_unit const&> (delta);

            Vector_3D_impl a (0,0,0);

            point[2]->space_position(). difference 
                (point[1]->space_position(), &a);

            if (point[1]->type() == "Atom" && 
                point[2]->type() == "Atom")
            {
                Atom & atom_1 = model.atom (point[1]->number());
                Atom & atom_2 = model.atom (point[2]->number());

                if (atom_1.bound_with (atom_2))
                {
                    Branch branch (atom_2, atom_1);

                    for (int i=0;  i<branch.atom_count();  ++i)
                    {
                        branch.atom(i).
                            rotate (a, point[1]->space_position(), 
                                    delta_unit.rad());
                    }
                    break;
                }
            }

            int count = moving_object->atom_count ();
            for (int i=0;  i<count;  ++i)
                moving_object->atom(i).
                    rotate (a, point[1]->space_position(), delta_unit.rad());

            break;
        }

        default: FLAW ("Active points are out of range.");
    }

    get_info (info);
}

Array <Measure_geometry_impl::Point const*> Measure_geometry_impl::
active_points () const
{
    INVARIANT;

    Array <Point const*> result;
    for (int i=0;  i<4;  ++i)
        if (is_active(i))
            result.push_back (&four_points_[i]);

    return result;
}

int Measure_geometry::
type_counter ()
{
    return 3;
}

Text const & Measure_geometry::
type_name (int index)
{
    static const Text name_0("Atom");
    static const Text name_1("Mass center");
    static const Text name_2("Space position");

    switch (index)
    {
    case 0:  return name_0;
    case 1:  return name_1;
    case 2:  return name_2;
    default: FLAW ("Out of range.");
    }
}

int Measure_geometry::
type_index (Text const & name)
{
    if       (name == "Atom")
        return 0;
    else  if (name == "Mass center")
        return 1;
    else  if (name == "Space position")
        return 2;
    else
        FLAW (Text("Wrong name - \"") + name + "\".");
}

void Measure_geometry_impl::
calc_space_positions ()
{
    if (is_active(0))
        calc_space_position (four_points_[0]);

    if (is_active(1))
        calc_space_position (four_points_[1]);

    if (is_active(2))
        calc_space_position (four_points_[2]);

    if (is_active(3))
        calc_space_position (four_points_[3]);
}

void Measure_geometry_impl::
calc_space_position (Point & point)
{
    if (!point.is_active())
        return;

    if (point.type() == "Atom")
    {
        point.set_space_position 
            (Project::singleton ().current_model().
                atom (point.number()).position());
    }

    else if (point.type() == "Mass center")
    {
        Point_3D_impl center(0,0,0);
        Mass_center mc(Project::singleton().current_model().
            molecule(point.number()));
        mc.evaluate (center);
        point.set_space_position (center);
    }

    else if (point.type () == "Space position")
        ;

    else 
        FLAW ("Unknown type.");
}

bool Measure_geometry_impl::
support_restraint () const
{
    Array <Point const*> point = active_points ();

    if (point.size() == 2 &&
        point[0]->type() == "Atom" &&
        point[1]->type() == "Atom")

        return true;

    return false;
}

void Measure_geometry_impl::
new_distance_restraint (double R_eqv, double K_r) const
{
    Array <Point const*> points = active_points ();
    Model & model = Project::singleton ().current_model();
    Selector & selector = model.kit().selector();
    int n0 = points[0]->number();
    int n1 = points[1]->number();
    Text type; 
    type += element (points[0]);
    type += " ";
    type += n0;
    type += " - ";
    type += element (points[1]);
    type += " ";
    type += n1;
    Atom & atom0 = model.atom (n0);
    Atom & atom1 = model.atom (n1);
    selector.select (atom0, true);
    selector.select (atom1, true);

    Default_builder().new_distance_restraint 
        (model, atom0, atom1, type, R_eqv, K_r);
}

}//MM
