#ifndef ARRAY_HPP
#define ARRAY_HPP
#include <iostream>
#include "ADIMPM.hpp"

template <class Type>
class Array2D
{
  public:
     Array2D(const int a_i, const int a_j);
     ~Array2D();
     int Index(const int a_i, const int a_j);
     Type& Data(const int a_i, const int a_j);
  private:
     Type* _array;
     int* _idx;
     int _ni, _nj;

};


template <class Type>
inline Array2D<Type>::Array2D(const int a_i, const int a_j):
_ni(a_i),
_nj(a_j)
{
  _array = new Type [_ni*_nj];

  _idx = new int [_ni];
  _idx[0] = 0;
  for(int i=1;i<_ni;i++)
    _idx[i] = _idx[i-1] + _nj;
};

template <class Type>
inline Array2D<Type>::~Array2D()
{
  delete [] _idx;
  delete [] _array;
};


template <class Type>
inline int Array2D<Type>::Index(const int a_i, const int a_j)
{
  assert(a_i>=0 && a_i<_ni && a_j>=0 && a_j<_nj);
  return _idx[a_i]+a_j;
};

template <class Type>
inline Type& Array2D<Type>::Data(const int a_i, const int a_j)
{
  return _array[Index(a_i,a_j)];
};


/* -------------------------------------------------------- */
template <class Type>
class Array3D
{
  public:
     Array3D(const int a_i, const int a_j, const int a_k);
     ~Array3D();
     Type& Data(const int a_i, const int a_j, const int a_k);
     int Index(const int a_i, const int a_j, const int a_k);
  private:
     Type* _array;
     int  *_idxi, *_idxj;
     const int _ni, _nj, _nk;
};


template <class Type>
inline Array3D<Type>::Array3D(const int a_i, const int a_j, const int a_k):
_ni(a_i),
_nj(a_j),
_nk(a_k)
{
  _array = new Type [_ni*_nj*_nk];

  _idxi = new int [_ni];
  _idxi[0] = 0;
  for(int i=1;i<_ni;i++)
  {
    _idxi[i] = _idxi[i-1] + _nj*_nk;
  }
  _idxj = new int [_nj];
  _idxj[0] = 0;
  for(int j=1;j<_nj;j++)
  {
    _idxj[j] = _idxj[j-1] + _nk;
  }
};

template <class Type>
inline Array3D<Type>::~Array3D()
{
  delete [] _idxi;
  delete [] _idxj;
  delete [] _array;
};

template <class Type>
inline int Array3D<Type>::Index(const int a_i, const int a_j, const int a_k)
{
   if(!(a_i>=0 && a_i<_ni && a_j>=0 && a_j<_nj && a_k>=0 && a_k<_nk))
 { 
      std::cout<<"ai-aj-ak "<<a_i<<" "<<a_j<<" "<<a_k<<std::endl;
      std::cout<<"ni-nj-nk "<<_ni<<" "<<_nj<<" "<<_nk<<std::endl;
      std::cout<<"arraymode  "<<arraymode<<std::endl;
  }
  
  //assert(a_i>=0 && a_i<_ni && a_j>=0 && a_j<_nj && a_k>=0 && a_k<_nk);
  
  return _idxi[a_i] +_idxj[a_j] + a_k;
};

template <class Type>
inline Type& Array3D<Type>::Data(const int a_i, const int a_j, const int a_k)
{
  return _array[Index(a_i,a_j,a_k)];
};


/* ---------------------------------------------------------*/
template <class Type>
class Array4D
{
  public:
     Array4D(const int a_i, const int a_j, const int a_k, const int a_v);
     ~Array4D();
     Type& Data(const int a_i, const int a_j, const int a_k, const int a_v);
     int Index(const int a_i, const int a_j, const int a_k, const int a_v);
  private:
     Type* _array;
     int  *_idxi, *_idxj, *_idxk;
     const int _ni, _nj, _nk, _nv;
};


template <class Type>
inline Array4D<Type>::Array4D(
const int a_i, const int a_j, const int a_k, const int a_v):
_ni(a_i),
_nj(a_j),
_nk(a_k),
_nv(a_v)
{
  _array = new Type [_ni*_nj*_nk*_nv];

  _idxi = new int [_ni];
  _idxi[0] = 0;
  for(int i=1;i<_ni;i++)
  {
    _idxi[i] = _idxi[i-1] + _nj*_nk*_nv;
  }
  _idxj = new int [_nj];
  _idxj[0] = 0;
  for(int j=1;j<_nj;j++)
  {
    _idxj[j] = _idxj[j-1] + _nk*_nv;
  }
  _idxk = new int [_nk];
  _idxk[0] = 0;
  for(int k=1;k<_nk;k++)
  {
    _idxk[k] = _idxk[k-1] + _nv;
  }  
  
};

template <class Type>
inline Array4D<Type>::~Array4D()
{
  delete [] _idxi;
  delete [] _idxj;
  delete [] _idxk;
  delete [] _array;
};

template <class Type>
inline int Array4D<Type>::Index(
const int a_i, const int a_j, const int a_k, const int a_v
)
{
   /*
  assert(a_i>=0 && a_i<_ni && 
         a_j>=0 && a_j<_nj && 
         a_k>=0 && a_k<_nk &&
         a_v>=0 && a_v<_nv );
  */
 // assert(a_i>=0 && a_i<_ni);
  if( !(a_i>=0 && a_i<_ni)) 
    std::cout<<"Error i "<<a_i<<std::endl;
    
  if( !(a_j>=0 && a_j<_nj)) 
    std::cout<<"Error J "<<a_j<<std::endl;
 // assert(a_j>=0 && a_j<_nj);

  
    if( !(a_k>=0 && a_k<_nk)) 
    std::cout<<"Error k "<<a_k<<std::endl;
    assert(a_k>=0 && a_k<_nk);
  assert(a_v>=0 && a_v<_nv);
  return _idxi[a_i] +_idxj[a_j] + _idxk[a_k] + a_v;
};

template <class Type>
inline Type& Array4D<Type>::Data(
const int a_i, const int a_j, const int a_k, const int a_v
)
{
  return _array[Index(a_i,a_j,a_k,a_v)];
};
#endif