#include "Box_impl.h"

#include "Mach_eps.h"
#include "Model.h"
#include "Atom.h"
#include "Defs.h"
#include "Log.h"

#include <math.h>

namespace
{
    double nil () {return ::d_mach_eps_2;}
}

namespace MM
{
const double Box_impl::default_L_ = 40.;

Box_impl::
Box_impl ()
:
    o_(0.,0.,0.),
    a_(default_L_,0.,0.),
    b_(0.,default_L_,0.),
    c_(0.,0.,default_L_)
{
}

bool Box_impl::is_default () const
{
    return is_rectangular () &&
           a_mod() == default_L_ &&
           b_mod() == default_L_ &&
           c_mod() == default_L_;}

bool Box_impl::
is_rectangular () const
{
    return  fabs (alpha () - 90) < nil() &&
            fabs (beta  () - 90) < nil() &&
            fabs (gamma () - 90) < nil();
}

double Box_impl::
volume () const
{
    Vector_3D_impl v (a());
    v.cross (b());
    double result = v.dot (c());
    return fabs (result);
}

double Box_impl::
b_mod_rectangle () const
{
    return b_.length() * sin (__Pi__ * gamma() / 180);
}

double Box_impl::
c_mod_rectangle () const
{
    double cos_alpha = cos (__Pi__ * alpha () / 180);
    double cos_beta  = cos (__Pi__ * beta  () / 180);
    double cos_gamma = cos (__Pi__ * gamma () / 180);

    double V =
        a_mod()*b_mod()*c_mod()* sqrt(
            (1.- cos_alpha*cos_alpha - cos_beta*cos_beta - cos_gamma*cos_gamma)
            + 2. * (cos_alpha*cos_beta*cos_gamma)
        );

    return V / (a_mod_rectangle() * b_mod_rectangle());
}

double Box_impl::
alpha () const
{
    INVARIANT;
    return 180. * acos (b_.dot (c_) / (b_.length() * c_.length())) / __Pi__;
}

double Box_impl::
beta () const
{
    INVARIANT;
    return 180. * acos (a_.dot (c_) / (a_.length() * c_.length())) / __Pi__;
}

double Box_impl::
gamma () const
{
    INVARIANT;
    return 180. * acos (b_.dot (a_) / (b_.length() * a_.length())) / __Pi__;
}

void Box_impl::
set_o (Point_3D const & new_o)
{
    o_.set_zero();
    o_.affine (new_o, 1.);
}

void Box_impl::
set (Vector_3D const & new_a, Vector_3D const & new_b, Vector_3D const & new_c)
{
    REQUIRE ("", new_a.length() > nil() && 
                 new_b.length() > nil() && 
                 new_c.length() > nil());

    a_.clear();     b_.clear();     c_.clear();
    a_ += new_a;    b_ += new_b;    c_ += new_c;
    was_changed_.publish();
}

void Box_impl::
set (double new_a, double new_b, double new_c,
     double new_alpha, double new_beta, double new_gamma)
{
    REQUIRE ("", new_a     > nil() && new_b    > nil() && new_c     > nil());
    REQUIRE ("", new_alpha > nil() && new_beta > nil() && new_gamma > nil());

    new_alpha *= __Pi__ / 180.;
    new_beta  *= __Pi__ / 180.;
    new_gamma *= __Pi__ / 180.;

    a_.clear(); 
    a_ += Vector_3D_impl (new_a, 0., 0.);
    
    b_.clear(); 
    b_ += Vector_3D_impl (new_b * cos(new_gamma), new_b * sin(new_gamma), 0.);

    double cx = new_c * cos (new_beta);
    double cy = new_c * (cos(new_alpha) - cos(new_beta) * cos(new_gamma)) 
                / sin(new_gamma);
    double cz = sqrt (new_c*new_c - cx*cx - cy*cy);
    
    c_.clear(); 
    c_ += Vector_3D_impl (cx, cy, cz);

    was_changed_.publish();

    //log () << "a_mod()=" << a_mod();
    ENSURE ("", fabs (a_mod() - new_a)                     < nil());
    ENSURE ("", fabs (b_mod() - new_b)                     < nil());
    ENSURE ("", fabs (c_mod() - new_c)                     < nil());
    ENSURE ("", fabs (alpha() - new_alpha * 180. / __Pi__) < nil());
    ENSURE ("", fabs (beta () - new_beta  * 180. / __Pi__) < nil());
    ENSURE ("", fabs (gamma() - new_gamma * 180. / __Pi__) < nil());
}

void Box_impl::
set (double new_a, double new_b, double new_c)
{
    set (new_a, new_b, new_c, alpha(), beta(), gamma());
}

bool Box_impl::
invariant () const
{
    bool result =
        a_.length () > nil() &&
        b_.length () > nil() &&
        c_.length () > nil();

    if (!result) log ()
        << "\n"
        << "Error:"                                  << "\n"
        << "a_.length () > nil()   " << a_.length () << "\n"
        << "b_.length () > nil()   " << b_.length () << "\n"
        << "c_.length () > nil()   " << c_.length () << "\n"
        << "\n";

    return result;//*/
}

void Box_impl::
fit (Atom_group const & model, Shape shape, double margin)
{
    const double min_margin = 2.;
    int count = model.atom_count();
    a_.clear();     b_.clear();     c_.clear();

    if (count == 0)
    {
        margin = margin < min_margin ? min_margin : margin;
        o_.set_x (-margin); o_.set_y (-margin); o_.set_z (-margin);
        a_.set_x (margin);  b_.set_y (margin);  c_.set_z (margin);
    }
    else
    {
        double x_max = model.atom (0).position().x();
        double y_max = model.atom (0).position().y();
        double z_max = model.atom (0).position().z();
        double x_min = x_max;
        double y_min = y_max;
        double z_min = z_max;

        for (int i=1;  i<count;  ++i)
        {
            Point_3D const& atom = model.atom (i).position();
            double x = atom.x(), y = atom.y(), z = atom.z();
            x_min = x_min < x ? x_min : x;
            y_min = y_min < y ? y_min : y;
            z_min = z_min < z ? z_min : z;
            x_max = x_max > x ? x_max : x;
            y_max = y_max > y ? y_max : y;
            z_max = z_max > z ? z_max : z;
        }

        double x = x_max - x_min + 2*margin;
        double y = y_max - y_min + 2*margin;
        double z = z_max - z_min + 2*margin;
        
        if (x > 2*min_margin)
        {
            o_.set_x (x_min - margin);   
            a_.set_x (x);
        }
        else
        {
            o_.set_x (x_min - min_margin);
            a_.set_x (2*min_margin);
        }
        
        if (y > 2*min_margin)
        {
            o_.set_y (y_min - margin);   
            b_.set_y (y);
        }
        else
        {
            o_.set_y (y_min - min_margin);
            b_.set_y (2*min_margin);
        }
        
        if (z > 2*min_margin)
        {
            o_.set_z (z_min - margin);   
            c_.set_z (z);
        }
        else
        {
            o_.set_z (z_min - min_margin);
            c_.set_z (2*min_margin);
        }
    }

    if (shape == Cube)
    {
        double max = a_mod() > b_mod() ? a_mod() : b_mod();
               max = c_mod() > max     ? c_mod() : max;
        set (max, max, max, 90, 90, 90);
    }

    was_changed_.publish();
    INVARIANT;
}

}//MM
