#include "SHRiMPS/Eikonals/Grid.H"
#include "ATOOLS/Org/Message.H"

using namespace SHRIMPS;
using namespace ATOOLS;

Grid::Grid(const double & Ymax) : m_Ymax(Ymax) {}

Grid::~Grid() {
  for (int i=0;i<m_ff1steps+1;i++) {
    for (int j=0;j<m_ff2steps+1;j++) {
      m_grid[i][j].clear();
    }
    m_grid[i].clear();
  }
  m_grid.clear();
}

void Grid::Initialise(const size_t & ff1steps,const size_t & ff2steps,
		      const double & ff1max,const double & ff2max) 
{
  m_ff1steps = ff1steps-1;
  m_ff2steps = ff2steps-1;
  m_grid.resize(m_ff1steps+1);
  for (int i=0;i<m_ff1steps+1;i++) m_grid[i].resize(m_ff2steps+1);
  m_ff1max   = ff1max;
  m_ff2max   = ff2max;
  m_deltaff1 = m_ff1max/double(m_ff1steps);  
  m_deltaff2 = m_ff2max/double(m_ff2steps);  
}

void Grid::InsertValues(const size_t & i,const size_t & j,
			const std::vector<double> & values) {
  //msg_Out()<<METHOD<<"("<<i<<", "<<j<<" -> "<<values.size()<<") : "
  //	   <<values[0]<<" "<<values[values.size()/2]<<" "
  //	   <<values[values.size()-1]<<".\n";
  m_grid[i][j].resize(values.size());
  copy(values.begin(),values.end(),m_grid[i][j].begin());
  m_deltay = 2.*m_Ymax/double(m_grid[i][j].size()-1);
}

double Grid::
operator()(const double & ff1,const double & ff2,const double & y) {
  if (ff1==0. || ff2==0.) return 0.;
  if (!FixBins(ff1,ff2))  return 0.;
  if (y<=-m_Ymax) return ValueAtLowerYEdge();
  if (y>=m_Ymax)  return ValueAtUpperYEdge();
  return Value(y);
}

bool Grid::FixBins(const double & ff1,const double & ff2) {
  m_ff1bin = (m_ff1max-ff1)/m_deltaff1;
  m_ff2bin = (m_ff2max-ff2)/m_deltaff2;
  if (m_ff1bin>=m_grid.size()-1 || m_ff2bin>=m_grid[0].size()-1) {
    msg_Error()<<"Error in "<<METHOD<<"("<<ff1<<", "<<ff2<<"):\n"
	       <<"   "<<m_ff1bin<<"/"<<m_ff2bin<<" from "
	       <<m_ff1max<<"/"<<m_ff2max<<" and "
	       <<m_deltaff1<<"/"<<m_deltaff2<<" vs. sizes "
	       <<m_grid.size()<<"/"<<m_grid[0].size()<<".\n";
    return false;
  }
  m_ff1low = m_ff1max-m_ff1bin*m_deltaff1; 
  m_ff1up  = m_ff1max-(m_ff1bin+1)*m_deltaff1;
  m_ff2low = m_ff2max-m_ff2bin*m_deltaff2;
  m_ff2up  = m_ff2max-(m_ff2bin+1)*m_deltaff2;
  m_d1up   = m_ff1up-ff1;
  m_d1low  = ff1-m_ff1low;
  m_d2up   = m_ff2up-ff2;
  m_d2low  = ff2-m_ff2low;
  return true;
}

double Grid::ValueAtLowerYEdge() {
  return
    (m_d1low * m_d2low * m_grid[m_ff1bin+1][m_ff2bin+1][0]+
     m_d1low * m_d2up  * m_grid[m_ff1bin+1][m_ff2bin+0][0]+
     m_d1up  * m_d2low * m_grid[m_ff1bin+0][m_ff2bin+1][0]+
     m_d1up  * m_d2up  * m_grid[m_ff1bin+0][m_ff2bin+0][0])/
    (m_deltaff1*m_deltaff2);
}

double Grid::ValueAtUpperYEdge() {
  size_t ylast(m_grid[0][0].size()-1);
  return
    (m_d1low * m_d2low * m_grid[m_ff1bin+1][m_ff2bin+1][ylast]+
     m_d1low * m_d2up  * m_grid[m_ff1bin+1][m_ff2bin+0][ylast]+
     m_d1up  * m_d2low * m_grid[m_ff1bin+0][m_ff2bin+1][ylast]+
     m_d1up  * m_d2up  * m_grid[m_ff1bin+0][m_ff2bin+0][ylast])/
    (m_deltaff1*m_deltaff2);
}

double Grid::Value(const double & y) {
  size_t ybin((y+m_Ymax)/m_deltay);
  double ylow(-m_Ymax+ybin*m_deltay), yup(-m_Ymax+(ybin+1)*m_deltay);
  double dyup(yup-y), dylow(y-ylow);
  if (m_ff1bin>=m_grid.size()-1 || m_ff2bin>=m_grid[0].size()-1 || 
      ybin>=m_grid[0][0].size()-1) {
    msg_Error()<<"Error in "<<METHOD<<".\n";
    return 0.;
  }
  return 
    (m_d1low * m_d2low * dylow * m_grid[m_ff1bin+1][m_ff2bin+1][ybin+1]+
     m_d1low * m_d2low * dyup  * m_grid[m_ff1bin+1][m_ff2bin+1][ybin+0]+
     m_d1low * m_d2up  * dylow * m_grid[m_ff1bin+1][m_ff2bin+0][ybin+1]+
     m_d1low * m_d2up  * dyup  * m_grid[m_ff1bin+1][m_ff2bin+0][ybin+0]+
     m_d1up  * m_d2low * dylow * m_grid[m_ff1bin+0][m_ff2bin+1][ybin+1]+
     m_d1up  * m_d2low * dyup  * m_grid[m_ff1bin+0][m_ff2bin+1][ybin+0]+
     m_d1up  * m_d2up  * dylow * m_grid[m_ff1bin+0][m_ff2bin+0][ybin+1]+
     m_d1up  * m_d2up  * dyup  * m_grid[m_ff1bin+0][m_ff2bin+0][ybin+0])/
    (m_deltaff1*m_deltaff2*m_deltay);
}