#ifndef _LATTICE_TYPE_H
#define _LATTICE_TYPE_H

struct lattice_type
{
/* must have serialize: AND in this order/syntax: http://www.boost.org/doc/libs/1_47_0/libs/serialization/doc/ */
#ifdef BOOST_MPI_HPP
private:

    friend class boost::serialization::access;
#endif
/* BOOST_MPI_HPP */

public:
    lattice_type* neighbors[ en - 1 ]; 

/* must have serialize; helps too w/IO? https://groups.google.com/group/boost-list/browse_thread/thread/95a8eb0570fd6ffe?hl=en */
#ifdef BOOST_MPI_HPP
protected:
    
    template<class Archive> //serializes (boost::mpi::packed_iarchive& ar) & deserializes (boost::mpi::packed_oarchive& ar)
    inline void serialize( Archive & ar , const unsigned int ) //http://www.boost.org/doc/libs/1_46_1/doc/html/mpi/tutorial.html
    {
        ar & neighbors ; // for 'packing' (& unpacking) for message-passing: put together in 'series' to exchange
    }

public:
#endif
/* BOOST_MPI_HPP */
  //KGK By defining elements/fields as static, it becomes unique (only 1 set for all elements in the vector)
  static int center_x;
  static int center_y;
  static int center_z;
  static int height;
  static int width;
  static int inlettype;
  static double inletradius;
  static double diameter;
  static xyz_type<int> dims;
  static std::vector<lattice_type>::iterator begin_iterator;
  static xyz_type<double> hh;       //KGK lattice length, width and height

  static int offset(const int x, const int y, const int z) 
  {
    return x * (dims.y * dims.z) + y * (dims.z) + z; //KGK x={0..dims.x-1}, y.., z..
  }
  
  xyz_type<int> idx;

  //fluid species[num_species]; 	// Fluid-specific kinematic, thermodynamic & electromagnetic properties
  double rho, T, p ;  //, s ; 
  PDF_type fi;
  PDF_type old_fi;
  xyz_dt u;  
  xyz_dt old_u;  
  xyz_dt gradrho;  //drho    (gradient of rho) used for FFT  - to determine electric field using spectral method
  xyz_dt J;        //J       (current density) used for FFT  - to determine magnetic field using spectral method
  xyz_dt curlJ;    //curlJ                     used for FFT  - to determine magnetic field using spectral method
	
  double Phi ;
  vec_PDF_type g;
  vec_PDF_type old_g;
  xyz_dt B;

  //KGK Variables used for the spectral method - START
  int mm,nn,rr;    //ID for the lattice in x,y,z directions
  int neighborID[ en - 1 ][3];
  int onBoundaryTag;  //If the cell is on the boundary then set onBoundary to 1
  xyz_dt A;        //Magnetic potential
  xyz_dt E;        //Electric field
  xyz_dt realFmnr; //Real part of Fmnr         used for IFFT - to determine magnetic/electric field using spectral method
  xyz_dt ImgFmnr;  //Imiginary part of Fmnr    used for IFFT - to determine magnetic/electric field using spectral method
  xyz_dt old_B;    // Testing FD solver for B bfunctions.h
  xyz_dt old_E;    // Testing FD solver for E bfunctions.h
  //KGK Variables used for the spectral method - END

  lattice_type() : u(0) { }

  bool is_on_face(int i) 
  {
    switch(i) 
    {
      case Xn: return (idx.x == 0);
      case Xp: return (idx.x == dims.x-1);
      case Yn: return (idx.y == 0);
      case Yp: return (idx.y == dims.y-1);
      case Zn: return (idx.z == 0);
      case Zp: return (idx.z == dims.z-1);
	  default: return false ; 
    } 
    return false ;
  }

  bool is_at_corner(void) const //KGK
  {
      return (((idx.x == 0 || idx.x == dims.x-1) &&
          (idx.y == 0 || idx.y == dims.y-1) &&
		  (idx.z == 0 || idx.z == dims.z-1)));	  
  }
	
  bool is_center_x(void) const 
  {
	return (idx.x == center_x);
  }
	
  bool is_center_y(void) const 
  {
    return (idx.y == center_y);
  }

  bool is_center_z(void) const
  {
    return (idx.z == center_z);
  }

  bool is_on_centerline(void) const
  {
    return (is_center_y() && is_center_z());
  }

  bool is_velinlet(void) const
  {
    const double distance_y = fabs(double(idx.y-center_y));
    const double distance_z = fabs(double(idx.z-center_z));

    if(inlettype == 0)
    {
    if (idx.x == 0 && distance_y<double(width/2.0) && distance_z<double(height/2.0))
    {
      return true;
    }
    return false;
    }

    else
    {
    if (idx.x == 0 && inletradius*inletradius > (distance_y*distance_y + distance_z*distance_z) )
    {
      return true;
    }
    return false;
    }
  }
  


};

int lattice_type::center_x       = 0;
int lattice_type::center_y       = 0;
int lattice_type::center_z       = 0;
int lattice_type::height         = 0;
int lattice_type::width          = 0;
double lattice_type::inletradius = 0.0;
int lattice_type::inlettype      = 0;
double lattice_type::diameter    = 0.0;
xyz_type<int> lattice_type::dims = 0;
std::vector<lattice_type>::iterator lattice_type::begin_iterator = 
  std::vector<lattice_type>::iterator();
xyz_type<double> lattice_type::hh = 0.0;

inline std::ostream& operator<< (std::ostream& s, const lattice_type &lattice)
{
  return s << lattice.u << "  " << std::setprecision(9) << std::setw(15) 
	   << lattice.rho << "  " << lattice.T << "  " << lattice.p << "  " << lattice.B << "  " << lattice.E; //<< lattice.s << "  " 
};

#endif