#ifndef ATOM_GRID_H
#define ATOM_GRID_H

#include <list>
#include <math.h>

#ifndef GRID_3D_H
#include "Grid_3D.h"
#endif

#ifndef BOX_IMPL_H
#include "Box_impl.h"
#endif

#ifndef ATOM_H
#include "Atom.h"
#endif

#ifndef ATOM_GROUP_H
#include "Atom_group.h"
#endif

#ifndef ARBITRARY_ATOM_GROUP_H
#include "Arbitrary_atom_group.h"
#endif

#ifndef DBC_H
#include "DbC.h"
#endif

#include "Length.h"
#include "Mach_eps.h"

namespace MM
{

class Atom_grid : public Grid_3D <std::list <Atom*> >, protected DbC
{
    double              lIcell_size_;
    Box_impl            box_;

public:
    class Near;
    friend class Near;

    explicit Atom_grid (Atom_group & atoms, double cell_size = 2.5)
        : lIcell_size_(1. / cell_size)
    {
        box_.fit (atoms, Box::Cuboid, cell_size);

        int x_size = (int)ceil (box_.a_mod () * lIcell_size_); 
        int y_size = (int)ceil (box_.b_mod () * lIcell_size_); 
        int z_size = (int)ceil (box_.c_mod () * lIcell_size_); 
        CHECK("", x_size > box_.a_mod() * lIcell_size_ - d_mach_eps);
        CHECK("", y_size > box_.b_mod() * lIcell_size_ - d_mach_eps);
        CHECK("", z_size > box_.c_mod() * lIcell_size_ - d_mach_eps);

        resize(x_size, y_size, z_size);

        int count = atoms.atom_count();
        for (int i=0;  i<count;  ++i)
        {
            Point_3D const & position = atoms.atom(i).position();
            int x = (int)floor ((position.x () - box_.o().x()) * lIcell_size_);
            int y = (int)floor ((position.y () - box_.o().y()) * lIcell_size_);
            int z = (int)floor ((position.z () - box_.o().z()) * lIcell_size_);
            CHECK("", x >= 0 && x < x_size);
            CHECK("", y >= 0 && y < y_size);
            CHECK("", z >= 0 && z < z_size);

            operator () (x, y, z). push_back (&atoms.atom(i));
        }
    }

    double              cell_size   () const    {return 1. / lIcell_size_;}
    double              lIcell_size () const    {return lIcell_size_;}


    class Near : public Arbitrary_atom_group
    {
    public:
        Near (Atom_grid const & grid, Point_3D const & point, double distance)
        {
            int x = (int)floor ((point.x () - grid.box_.o().x()) * grid.lIcell_size_);
            int y = (int)floor ((point.y () - grid.box_.o().y()) * grid.lIcell_size_);
            int z = (int)floor ((point.z () - grid.box_.o().z()) * grid.lIcell_size_);

            int shift = 1 + (int)ceil (distance * grid.lIcell_size_);

            int ib = x-shift < 0 ? 0 : x-shift;
            int jb = y-shift < 0 ? 0 : y-shift;
            int kb = z-shift < 0 ? 0 : z-shift;

            int ie = x+shift < grid.x_size() ? x+shift : grid.x_size()-1;
            int je = y+shift < grid.y_size() ? y+shift : grid.y_size()-1;
            int ke = z+shift < grid.z_size() ? z+shift : grid.z_size()-1;

            for (int i=ib;  i<=ie;  ++i)
            for (int j=jb;  j<=je;  ++j)
            for (int k=kb;  k<=ke;  ++k)
            {   
                std::list <Atom*> list = grid(i, j, k);
                std::list <Atom*>::const_iterator it  = list.begin();
                std::list <Atom*>::const_iterator end = list.end  ();

//                int size = list.size();

                for (;  it!= end;  ++it)
                    add (**it);
            }
        }
    };

};

}//MM

#endif //ATOM_GRID_H
