#include "sphereecp.h"
#include <math.h>

namespace SOPHYA {

template <class T>
SphereECP<T>::SphereECP()
  : _pixels()
{
  _partial = false;
  _theta1 = 0.;
  _theta2 = M_PI;
  _phi1 = 0.;
  _phi2 = 2*M_PI;
  _outofmapidx = -99;
  _outofmappix = 0;
  _outofmapval = 0;
  _dtheta = 0.;
  _dphi = 0.;
}

template <class T>
SphereECP<T>::SphereECP(int m)
  : _pixels(2*m,m,0,0,0,true)
{
  _partial = false;
  _theta1 = 0.;
  _theta2 = M_PI;
  _phi1 = 0.;
  _phi2 = 2*M_PI;
  _outofmapidx = -99;
  _outofmappix = 0;
  _outofmapval = 0;
  _dtheta = _dphi = M_PI/m;
}

template <class T>
SphereECP<T>::SphereECP(int ntet, int nphi)
  : _pixels(nphi,ntet,0,0,0,true)
{
  _partial = false;
  _theta1 = 0.;
  _theta2 = M_PI;
  _phi1 = 0.;
  _phi2 = 2*M_PI; 
  _outofmapidx = -99;
  _outofmappix = 0;
  _outofmapval = 0;
  _dtheta = M_PI/ntet;
  _dphi = 2*M_PI/nphi;
}

template <class T>
SphereECP<T>::SphereECP(r_8 tet1, r_8 tet2, int ntet, r_8 phi1, r_8 phi2, int nphi)
  : _pixels(nphi,ntet,0,0,0,true)
{
  if( (tet1 > M_PI) || (tet1 < 0.) || (tet2 > M_PI) || (tet2 < 0.) || 
      (phi1 < 0.) || (phi1 > 2*M_PI) || (phi1 < 0.) || (phi1 > 2*M_PI) ||
      (tet2 <= tet1) || (phi2 <= phi1) )
    throw  ParmError("SphereECP::SphereECP() Bad range for theta/phi limits");

   _partial = true;
  _theta1 = tet1;
  _theta2 = tet2;
  _phi1 = phi1;
  _phi2 = phi2; 
  _outofmapidx = _pixels.Size()+659;
  _outofmappix = 0;
  _outofmapval = 0;
  _dtheta = (_theta2-_theta1)/ntet;
  _dphi = (_phi2-_phi1)/nphi;
}

template <class T>
SphereECP<T>::SphereECP(const SphereECP<T>& a, bool share)
  : _pixels(a._pixels, share)
{
  _partial = a._partial;
  _theta1 = a._theta1;
  _theta2 = a._theta2;
  _phi1 = a._phi1;
  _phi2 = a._phi2;   
  _outofmapidx = a._outofmapidx;
  _outofmappix = a._outofmappix;
  _outofmapval = a._outofmapval;
  _dtheta = a._dtheta;
  _dphi = a._dphi;
}

template <class T>
SphereECP<T>::SphereECP(const SphereECP<T>& a)
  : _pixels(a._pixels)
{
  _partial = a._partial;
  _theta1 = a._theta1;
  _theta2 = a._theta2;
  _phi1 = a._phi1;
  _phi2 = a._phi2;   
  _outofmapidx = a._outofmapidx;
  _outofmappix = a._outofmappix;
  _outofmapval = a._outofmapval;
  _dtheta = a._dtheta;
  _dphi = a._dphi;
}

template <class T>
SphereECP<T>::~SphereECP()
{
}

template <class T>
string SphereECP<T>::TypeOfMap() const
{
  return string("TETAFI");
}

template <class T>
void SphereECP<T>::SetTemp(bool temp) const
{
  _pixels.SetTemp(temp);
}

template <class T>
int_4 SphereECP<T>::NbPixels() const
{
  return _pixels.Size();
}

template <class T>
T& SphereECP<T>::PixVal(int_4 k)
{
  if ((_partial) && (k==_outofmapidx))   return _outofmappix;
  if((k < 0) || (k >= _pixels.Size()) )  
     throw RangeCheckError("SphereECP::PIxVal() Pixel index out of range");
  return _pixels.DataBlock()(k);
}

template <class T>
T const& SphereECP<T>::PixVal(int_4 k)  const
{
  if ((_partial) && (k==_outofmapidx))   return _outofmapval;
  if((k < 0) || (k >= _pixels.Size()) )  
     throw RangeCheckError("SphereECP::PIxVal() Pixel index out of range");
  //  return _pixels.DataBlock()(k);
  return _pixels.Data()[k];
}

template <class T>
bool  SphereECP<T>::ContainsSph(double theta, double phi) const
{
  if (_partial) {
    if ( (theta >= _theta1) && (theta <= _theta2) &&
	 (phi >= _phi1) && (phi <= _phi2) )  return true;
    else return false;
  }
  else return false;
}

template <class T>
int_4  SphereECP<T>::PixIndexSph(double theta, double phi) const
{
  if((theta > M_PI) || (theta < 0.)) return(-1);
  if((phi < 0.) || (phi > 2*M_PI)) return(-1);
  if (_partial) {
    if ( (theta < _theta1) || (theta > _theta2) ||
	 (phi < _phi1) || (phi > _phi2) )  return _outofmapidx;
  }
  int_4 it = (theta-_theta1)/_dtheta;
  int_4 jp = (phi-_phi1)/_dphi;
  return (it*_pixels.SizeX()+jp);
}

template <class T>
void  SphereECP<T>::PixThetaPhi(int_4 k, double& theta, double& phi) const
{
  if ((_partial) && (k==_outofmapidx)) {
    theta = M_PI;
    phi = 2*M_PI;
    return;
  }
  if((k < 0) || (k >= _pixels.Size()) )  
     throw RangeCheckError("SphereECP::PixThetaPhi() Pixel index out of range");

  int_4 it = k / _pixels.SizeX();
  int_4 jp = k % _pixels.SizeX();
  theta = _theta1+it*_dtheta;
  phi = _phi1*jp*_dphi;
  return;
}

template <class T>
T  SphereECP<T>::SetPixels(T v)
{
  _pixels = v;
  return v;
}

template <class T>
double  SphereECP<T>::PixSolAngle(int_4 k) const
{
  if ((k < 0) || (k >= _pixels.Size()))  return (_dtheta*_dphi);
  int_4 it = k / _pixels.SizeX();
  double theta = _theta1+it*_dtheta;
  return (_dtheta*_dphi*sin(theta));
}

template <class T>
int_4  SphereECP<T>::SizeIndex() const
{
  return _pixels.SizeY();
}

template <class T>
void  SphereECP<T>::Resize(int_4 m)
{
  if ( (m <= 0) || ( m == _pixels.SizeY()) ) { 
    cout << " SphereECP<T>::Resize(int_4 " << m << ") m<0 ou m=NTheta - Ne fait rien " << endl;
    return;
  }
  int mphi = m;
  if (_pixels.Size() > 0)  mphi = m*_pixels.SizeX()/_pixels.SizeY();
  cout << " SphereECP<T>::Resize(" << m 
       << ") -> _pixels.Resize(" << mphi << "," << m << ")" << endl; 
  sa_size_t sz[5] = {0,0,0,0,0};
  sz[0] = mphi;  sz[1] = m;
  _pixels.ReSize(2,sz);
  _outofmapidx = _pixels.Size()+659;
  _dtheta = (_theta2-_theta1)/m;
  _dphi = (_phi2-_phi1)/mphi;
  return;
}

template <class T>
uint_4  SphereECP<T>::NbThetaSlices() const
{
  return _pixels.SizeY();
}

template <class T>
void  SphereECP<T>::GetThetaSlice(int_4 index,r_8& theta, 
				  TVector<r_8>& phi, TVector<T>& value) const 
{
  if( (index < 0) || (index >= _pixels.SizeY()) )  
     throw RangeCheckError("SphereECP::GetThetaSlice() theta index out of range");

  theta = _theta1 + index*_dtheta;
  int_4 nphit = 2.*M_PI/_dphi;
  phi.ReSize(nphit); 
  value.ReSize(nphit); 
  int_4 ioff = _phi1/_dphi;
  int_4 i;
  for(i=0; i<ioff; i++) {
    phi(i) = _phi1 + (i-ioff)*_dphi;
    value(i) = _outofmapval;
  }
  for(i=ioff; i<ioff+_pixels.SizeX(); i++) {
    phi(i) = _phi1 + (i-ioff)*_dphi;
    value(i) = _pixels(i-ioff,index,0);
  }
  for(i=ioff+_pixels.SizeX(); i<nphit; i++) {
    phi(i) = _phi1 + (i-ioff)*_dphi;
    value(i) = _outofmapval;
  }
  return;
}

template <class T>
void  SphereECP<T>::GetThetaSlice(int_4 index, r_8& theta, r_8& phi0, 
				  TVector<int_4>& pixidx, TVector<T>& value) const
{
  if( (index < 0) || (index >= _pixels.SizeY()) )  
    throw RangeCheckError("SphereECP::GetThetaSlice() theta index out of range");
  
  theta = _theta1 + index*_dtheta;
  int_4 ioff = _phi1/_dphi;
  phi0 = _phi1 - ioff*_dphi;

  int_4 nphit = 2.*M_PI/_dphi;
  pixidx.ReSize(nphit); 
  value.ReSize(nphit);
 
  int_4 i;
  for(i=0; i<ioff; i++) {
    pixidx(i) =  _outofmapidx;
    value(i) = _outofmapval;
  }
  for(i=ioff; i<ioff+_pixels.SizeX(); i++) {
    pixidx(i) = index*_pixels.SizeX()+(i-ioff);
    value(i) = _pixels(i-ioff,index,0);
  }
  for(i=ioff+_pixels.SizeX(); i<nphit; i++) {
    pixidx(i) =  _outofmapidx;
    value(i) = _outofmapval;
  }
  return;

}

template <class T>
void  SphereECP<T>::Print(ostream& os) const
{
  if (_partial)
    os << "SphereECP<T>::Print() - Partial ECP Map NPhi=" << _pixels.SizeX() 
       << " x NTheta= " <<  _pixels.SizeX() << " SolidAngle=" << PixSolAngle() << "\n"
       << "ThetaLimits= " << _theta1 << " .. " << _theta2 
       << "  PhiLimits= " << _phi1 << " .. " << _phi2 << endl;
  else 
    os << "SphereECP<T>::Print() - Full ECP Map NPhi=" << _pixels.SizeX() 
       << " x NTheta= " <<  _pixels.SizeY() << " SolidAngle=" << PixSolAngle() << endl;
  os << _pixels;
}

template <class T>
SphereECP<T>&  SphereECP<T>::Set(const SphereECP<T>& a)
{
  _pixels.Set(a._pixels);
  _partial = a._partial;
  _theta1 = a._theta1;
  _theta2 = a._theta2;
  _phi1 = a._phi1;
  _phi2 = a._phi2;   
  _outofmapidx = a._outofmapidx;
  _outofmappix = a._outofmappix;
  _outofmapval = a._outofmapval;
  _dtheta = a._dtheta;
  _dphi = a._dphi; 
  return *this ;
}

template <class T>
SphereECP<T>&  SphereECP<T>::SetCst(T x)
{
  _pixels.SetCst(x);
  return *this ;
}

template <class T>
SphereECP<T>&  SphereECP<T>::AddCst(T x)
{
  _pixels += x;
  return *this ;
}

template <class T>
SphereECP<T>&  SphereECP<T>::MulCst(T x)
{
  _pixels *= x;
  return *this ;
}





#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template SphereECP<int_4>
#pragma define_template SphereECP<r_4>
#pragma define_template SphereECP<r_8>
#pragma define_template SphereECP< complex<r_4> >
#pragma define_template SphereECP< complex<r_8> >
#endif
#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class SphereECP<int_4>;
template class SphereECP<r_4>;
template class SphereECP<r_8>;
template class SphereECP< complex<r_4> >;
template class SphereECP< complex<r_8> >;
#endif

}// Fin du namespace