source: Sophya/trunk/SophyaLib/SkyMap/spherehealpix.cc@ 2808

#include "sopnamsp.h"
#include "machdefs.h"
#include <math.h>
#include <complex>

#include "pexceptions.h"
#include "fiondblock.h"
#include "spherehealpix.h"
#include "strutil.h"

extern "C" 
{
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
}
      
using namespace SOPHYA;

/*!
  \class SOPHYA::SphereHEALPix
  \ingroup SkyMap
  \brief Spherical maps in HEALPix pixelisation scheme.

  Class implementing spherical maps, in the HEALPix pixelisation scheme, 
  with template data types (double, float, complex, ...)


\verbatim

    -----------------------------------------------------------------------
     version 0.8.2  Aug97 TAC  Eric Hivon, Kris Gorski
    -----------------------------------------------------------------------

    the sphere is split in 12 diamond-faces containing nside**2 pixels each
    the numbering of the pixels (in the nested scheme) is similar to
    quad-cube
    In each face the first pixel is in the lowest corner of the diamond
    the faces are                    (x,y) coordinate on each face

        .   .   .   .   <--- North Pole
       / \ / \ / \ / \                          ^        ^     
      . 0 . 1 . 2 . 3 . <--- z = 2/3             \      / 
       \ / \ / \ / \ /                        y   \    /  x  
      4 . 5 . 6 . 7 . 4 <--- equator               \  /      
       / \ / \ / \ / \                              \/      
      . 8 . 9 .10 .11 . <--- z = -2/3              (0,0) : lowest corner 
       \ / \ / \ / \ /      
        .   .   .   .   <--- South Pole
\endverbatim
    phi:0               2Pi                                 

   in the ring scheme pixels are numbered along the parallels
   the first parallel is the one closest to the north pole and so on  
   on each parallel, pixels are numbered starting from the one closest
    to phi = 0

    nside MUST be a power of 2 (<= 8192)

*/

/* --Methode-- */

template<class T>
SphereHEALPix<T>::SphereHEALPix() : pixels_(), sliceBeginIndex_(), 
                                                sliceLenght_()


{
  InitNul();
}

/*! \fn   SOPHYA::SphereHEALPix::SphereHEALPix(int_4 m)
 
  \param <m> : "nside" of the Healpix algorithm

  The total number of pixels will be Npix =  12*nside**2

  nside MUST be a power of 2 (<= 8192)
*/

template<class T>
SphereHEALPix<T>::SphereHEALPix(int_4 m)

{

  if(m <= 0 || m > 8192) 
    {
      cout << "SphereHEALPix : m hors bornes [0,8192], m= " << m << endl;
      throw RangeCheckError("SphereHEALPix<T>::SphereHEALPix() - Out of bound nside (< 8192)!");
    }
  // verifier que m est une puissance de deux 
  int x= m;
  while(x%2 == 0) x/=2;
  if(x != 1) 
    {  
      cout<<"SphereHEALPix: m doit etre une puissance de deux, m= "<<m<<endl;
      throw ParmError("SphereHEALPix<T>::SphereHEALPix() - nside != 2^n !");
    }
  InitNul();
  Pixelize(m);
  SetThetaSlices();
}
//++
template<class T>
SphereHEALPix<T>::SphereHEALPix(const SphereHEALPix<T>& s, bool share)
  : pixels_(s.pixels_, share), sliceBeginIndex_(s.sliceBeginIndex_, share), 
                                sliceLenght_(s.sliceLenght_, share)
//    copy constructor
//--
{
  nSide_= s.nSide_;
  nPix_ = s.nPix_;
  omeg_ = s.omeg_;
  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
}
//++
template<class T>
SphereHEALPix<T>::SphereHEALPix(const SphereHEALPix<T>& s)
  : pixels_(s.pixels_), sliceBeginIndex_(s.sliceBeginIndex_), 
                                sliceLenght_(s.sliceLenght_)
//    copy constructor
//--
{
  nSide_= s.nSide_;
  nPix_ = s.nPix_;
  omeg_ = s.omeg_;
  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
  //  CloneOrShare(s);
}

//! Clone if \b a is not temporary, share if temporary
/*! \sa NDataBlock::CloneOrShare(const NDataBlock<T>&) */
template<class T>
void SphereHEALPix<T>::CloneOrShare(const SphereHEALPix<T>& a)
{
  nSide_= a.nSide_;
  nPix_ = a.nPix_;
  omeg_ = a.omeg_;
  pixels_.CloneOrShare(a.pixels_);
  sliceBeginIndex_.CloneOrShare(a.sliceBeginIndex_);
  sliceLenght_.CloneOrShare(a.sliceLenght_);
  if (mInfo_) {delete mInfo_; mInfo_ = NULL;}
  if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
}

//! Share data with a
template<class T>
void SphereHEALPix<T>::Share(const SphereHEALPix<T>& a)
{
  nSide_= a.nSide_;
  nPix_ = a.nPix_;
  omeg_ = a.omeg_;
  pixels_.Share(a.pixels_);
  sliceBeginIndex_.Share(a.sliceBeginIndex_);
  sliceLenght_.Share(a.sliceLenght_);
  if (mInfo_) {delete mInfo_; mInfo_ = NULL;}
  if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
}

////////////////////////// methodes de copie/share
template<class T>
SphereHEALPix<T>& SphereHEALPix<T>::Set(const SphereHEALPix<T>& a)
{
  if (this != &a) 
    { 
      
      if (a.NbPixels() < 1) 
        throw RangeCheckError("SphereHEALPix<T>::Set(a ) - SphereHEALPix a not allocated ! ");
      if (NbPixels() < 1) CloneOrShare(a);
      else CopyElt(a);
      
      
      if (mInfo_) delete mInfo_;
      mInfo_ = NULL;
      if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
    }
  return(*this);
}

template<class T>
SphereHEALPix<T>& SphereHEALPix<T>::CopyElt(const SphereHEALPix<T>& a)
{
  if (NbPixels() < 1) 
    throw RangeCheckError("SphereHEALPix<T>::CopyElt(const SphereHEALPix<T>& )  - Not Allocated SphereHEALPix ! ");
  if (NbPixels() != a.NbPixels()) 
    throw(SzMismatchError("SphereHEALPix<T>::CopyElt(const SphereHEALPix<T>&) SizeMismatch")) ;
  nSide_= a.nSide_;
  nPix_ = a.nPix_;
  omeg_ = a.omeg_;
  int k;
  for (k=0; k< nPix_; k++) pixels_(k) = a.pixels_(k);
  for (k=0; k< a.sliceBeginIndex_.Size(); k++) sliceBeginIndex_(k) = a.sliceBeginIndex_(k);
  for (k=0; k< a.sliceLenght_.Size(); k++) sliceLenght_(k) = a.sliceLenght_(k);
  return(*this);
}



//++
// Titre        Destructor
//--
//++
template<class T>
SphereHEALPix<T>::~SphereHEALPix()

//--
{
}

//++
// Titre        Public Methods
//--

/*!  \fn   void SOPHYA::SphereHEALPix::Resize(int_4 m)
  \param <m>   "nside" of the Gorski algorithm

  The total number of pixels will be Npix =  12*nside**2

  nside MUST be a power of 2 (<= 8192)
*/
template<class T>
void SphereHEALPix<T>::Resize(int_4 m)


{
  if (m<=0 || m> 8192) {
    cout << "SphereHEALPix : m hors bornes [0,8192], m= " << m << endl;
    exit(1);
  }
  // verifier que m est une puissance de deux 
  int x= m;
  while (x%2==0) x/=2;
  if(x != 1) 
    {  
      cout<<"SphereHEALPix: m doit etre une puissance de deux, m= "<<m<<endl;
      exit(1);
    }
  InitNul();
  Pixelize(m);
  SetThetaSlices();
}

template<class T>
void  SphereHEALPix<T>::Pixelize( int_4 m) 

//    prépare la pixelisation Gorski (m a la même signification 
//    que pour le constructeur)
//
//    
{
  // On memorise les arguments d'appel
  nSide_= m;  

  // Nombre total de pixels sur la sphere entiere
  nPix_= 12*nSide_*nSide_;

  // pour le moment les tableaux qui suivent seront ranges dans l'ordre 
  // de l'indexation GORSKY "RING"
  // on pourra ulterieurement changer de strategie et tirer profit
  // de la dualite d'indexation GORSKY (RING et NEST) : tout dependra 
  // de pourquoi c'est faire

  // Creation et initialisation du vecteur des contenus des pixels 
  pixels_.ReSize(nPix_);
  pixels_.Reset();

  // solid angle per pixel   
  omeg_= 4.0*Pi/nPix_;
}

template<class T>
void SphereHEALPix<T>::InitNul()
//
//    initialise à zéro les variables de classe 
{
  nSide_= 0;
  nPix_ = 0;
  omeg_ = 0.;
//  pixels_.Reset();  -  Il ne faut pas mettre les pixels a zero si share !
}

/* --Methode-- */
/* Nombre de pixels du decoupage */
/*! \fn int_4 SOPHYA::SphereHEALPix::NbPixels() const

   Return number of  pixels of the  splitting 
*/
template<class T>
int_4 SphereHEALPix<T>::NbPixels() const
{  
  return(nPix_);
}


/*! \fn uint_4 SOPHYA::SphereHEALPix::NbThetaSlices() const

   \return number of slices in theta direction on the  sphere 
*/
template<class T>
uint_4 SphereHEALPix<T>::NbThetaSlices() const 
{
  uint_4 nbslices = uint_4(4*nSide_-1);
  if (nSide_<=0) 
    {
      nbslices = 0;
      throw PException(" sphere not pixelized, NbSlice=0 ");
    }
  return nbslices;
}

/*!  \fn void SOPHYA::SphereHEALPix::GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const 

 For a theta-slice with index 'index', return :
 
   the corresponding "theta" 

   a vector containing the phi's of the pixels of the slice

   a vector containing the corresponding values of pixels 
*/
template<class T>
void SphereHEALPix<T>::GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const 

{

  if (index<0 || index >= NbThetaSlices()) 
    {
      cout << " SphereHEALPix::GetThetaSlice : Pixel index out of range" <<endl;
      throw RangeCheckError(" SphereHEALPix::GetThetaSlice : Pixel index out of range");
    }


  int_4 iring= sliceBeginIndex_(index);
  int_4 lring  =  sliceLenght_(index);

   phi.ReSize(lring);
   value.ReSize(lring);

  double TH= 0.;
  double FI= 0.;
  for(int_4 kk = 0; kk < lring;kk++) 
    {
      PixThetaPhi(kk+iring,TH,FI);
      phi(kk)= FI;
      value(kk)= PixVal(kk+iring);
    }
  theta= TH;
}
/*! \fn void SOPHYA::SphereHEALPix::GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const

   For a theta-slice with index 'index', return : 

   the corresponding "theta" 

   the corresponding "phi" for first pixel of the slice 

    a vector containing indices of the pixels of the slice

   (equally distributed in phi)

    a vector containing the corresponding values of pixels 
*/

template<class T>
void SphereHEALPix<T>::GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const 

{

  if (sliceIndex<0 || sliceIndex >= NbThetaSlices()) 
    {
      cout << " SphereHEALPix::GetThetaSlice : Pixel index out of range" <<endl;
      throw RangeCheckError(" SphereHEALPix::GetThetaSlice : Pixel index out of range");
    }
  int_4 iring= sliceBeginIndex_(sliceIndex);
  int_4 lring  =  sliceLenght_(sliceIndex);
  pixelIndices.ReSize(lring);
  value.ReSize(lring);

  for(int_4 kk = 0; kk < lring;kk++) 
    {
      pixelIndices(kk)= kk+iring;
      value(kk)= PixVal(kk+iring);
    }
  PixThetaPhi(iring, theta, phi0);
}

template<class T>
void SphereHEALPix<T>::SetThetaSlices()  
{
  sliceBeginIndex_.ReSize(4*nSide_-1);
  sliceLenght_.ReSize(4*nSide_-1);
  int sliceIndex;
  for (sliceIndex=0; sliceIndex<  nSide_-1; sliceIndex++)
    {
      sliceBeginIndex_(sliceIndex)  = 2*sliceIndex*(sliceIndex+1);
      sliceLenght_(sliceIndex) = 4*(sliceIndex+1);
    }
  for (sliceIndex= nSide_-1; sliceIndex<  3*nSide_; sliceIndex++)
    {
      sliceBeginIndex_(sliceIndex)  = 2*nSide_*(2*sliceIndex-nSide_+1);
      sliceLenght_(sliceIndex) = 4*nSide_;
    }
  for (sliceIndex= 3*nSide_; sliceIndex< 4*nSide_-1; sliceIndex++)
    {
      int_4 nc= 4*nSide_-1-sliceIndex;
      sliceBeginIndex_(sliceIndex)  = nPix_-2*nc*(nc+1); 
      sliceLenght_(sliceIndex) = 4*nc;
    }
}


/*! \fn T& SOPHYA::SphereHEALPix::PixVal(int_4 k)

   \return value of  pixel with "RING" index k 
*/
template<class T>
T& SphereHEALPix<T>::PixVal(int_4 k)

{
  if((k < 0) || (k >= nPix_)) 
    {
      throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range");
    }
  return pixels_(k);
}

/*! \fn T const& SOPHYA::SphereHEALPix::PixVal(int_4 k) const

   \return value of  pixel with "RING" index k 
*/
template<class T>
T const& SphereHEALPix<T>::PixVal(int_4 k) const

{
  if((k < 0) || (k >= nPix_)) 
    {
      throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range");
    }
  return *(pixels_.Data()+k);
}

/*! \fn T& SOPHYA::SphereHEALPix::PixValNest(int_4 k)

   \return value of  pixel with "NESTED" index k 
*/
template<class T>
T& SphereHEALPix<T>::PixValNest(int_4 k)

//--
{
  if((k < 0) || (k >= nPix_)) 
    { 
      throw RangeCheckError("SphereHEALPix::PIxValNest Pixel index out of range");
    }
  return pixels_(nest2ring(nSide_,k));
}

/*!  \fn T const& SOPHYA::SphereHEALPix::PixValNest(int_4 k) const

  \return value of  pixel with "NESTED" index k 
*/

template<class T>
T const& SphereHEALPix<T>::PixValNest(int_4 k) const

{
  if((k < 0) || (k >= nPix_)) 
    { 
      throw RangeCheckError("SphereHEALPix::PIxValNest Pixel index out of range");
  }
  int_4 pix= nest2ring(nSide_,k);
  return *(pixels_.Data()+pix);
}

/*! \fn bool SOPHYA::SphereHEALPix::ContainsSph(double theta, double phi) const

\return true if teta,phi in map 
*/
template<class T>
bool SphereHEALPix<T>::ContainsSph(double theta, double phi) const
{
return(true);
}

/*! \fn int_4 SOPHYA::SphereHEALPix::PixIndexSph(double theta,double phi) const

 \return "RING" index of the pixel corresponding to direction (theta, phi).
 */
template<class T>
int_4 SphereHEALPix<T>::PixIndexSph(double theta,double phi) const

{
  return ang2pix_ring(nSide_,theta,phi);
}

/*! \fn int_4 SOPHYA::SphereHEALPix::PixIndexSphNest(double theta,double phi) const

 \return "NESTED" index of the pixel corresponding to direction (theta, phi).
 */
template<class T>
int_4 SphereHEALPix<T>::PixIndexSphNest(double theta,double  phi) const

{
  return ang2pix_nest(nSide_,theta,phi);
}


/* --Methode-- */
/*! \fn void SOPHYA::SphereHEALPix::PixThetaPhi(int_4 k,double& theta,double& phi) const
 \return (theta,phi) coordinates of middle of  pixel with "RING" index k
 */
template<class T>
void SphereHEALPix<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
{
  pix2ang_ring(nSide_,k,theta,phi);
}

/*! \fn T SOPHYA::SphereHEALPix::SetPixels(T v)

Set all pixels to value v 
*/
template <class T>
T SphereHEALPix<T>::SetPixels(T v)
{
pixels_.Reset(v);
return(v);
}



/*! \fn void SOPHYA::SphereHEALPix::PixThetaPhiNest(int_4 k,double& theta,double& phi) const

  \return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
 */
template<class T>
void SphereHEALPix<T>::PixThetaPhiNest(int_4 k,double& theta,double& phi) const

{   
  pix2ang_nest(nSide_,k,theta,phi);
}


/*! \fn int_4 SOPHYA::SphereHEALPix::NestToRing(int_4 k) const

   translation from NESTED index  into RING index 
*/
template<class T>
int_4 SphereHEALPix<T>::NestToRing(int_4 k) const

{
  return  nest2ring(nSide_,k);
}


/*!  \fn  int_4 SOPHYA::SphereHEALPix::RingToNest(int_4 k) const

 translation from  RING index  into NESTED index 
*/
template<class T>
int_4 SphereHEALPix<T>::RingToNest(int_4 k) const
{
  return  ring2nest(nSide_,k);
}

//   ...... Operations de calcul  ......

//! Fill a SphereHEALPix with a constant value \b a
template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::SetT(T a) 
{
  if (NbPixels() < 1) 
    throw RangeCheckError("SphereHEALPix<T>::SetT(T )  - SphereHEALPix not dimensionned ! ");
  pixels_ = a;
  return (*this);
}

/*! Add a constant value \b x to a SphereHEALPix */
template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::Add(T a)
 {
   cout << " c'est mon Add " << endl;
  if (NbPixels() < 1) 
    throw RangeCheckError("SphereHEALPix<T>::Add(T )  - SphereHEALPix not dimensionned ! ");
  //   pixels_ += a; 
   pixels_.Add(a); 
  return (*this);
}

/*! Substract a constant value \b a to a SphereHEALPix */
template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::Sub(T a,bool fginv) 
{
  if (NbPixels() < 1) 
    throw RangeCheckError("SphereHEALPix<T>::Sub(T )  - SphereHEALPix not dimensionned ! ");
  pixels_.Sub(a,fginv); 
  return (*this);
} 

/*! multiply a SphereHEALPix by a constant value \b a */
template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::Mul(T a) 
{
  if (NbPixels() < 1) 
    throw RangeCheckError("SphereHEALPix<T>::Mul(T )  - SphereHEALPix not dimensionned ! ");
  pixels_ *= a; 
  return (*this);
}

/*! divide a SphereHEALPix by a constant value \b a */
template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::Div(T a) 
{
  if (NbPixels() < 1) 
    throw RangeCheckError("SphereHEALPix<T>::Div(T )  - SphereHEALPix not dimensionned ! ");
  pixels_ /= a; 
  return (*this);
} 

//  >>>> Operations avec 2nd membre de type SphereHEALPix
//! Add two SphereHEALPix

template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::AddElt(const SphereHEALPix<T>& a)
{
  if (NbPixels() != a.NbPixels() )
    {
    throw(SzMismatchError("SphereHEALPix<T>::AddElt(const SphereHEALPix<T>&) SizeMismatch")) ;
    }
  pixels_ += a.pixels_;
  return (*this);
}

//! Substract two SphereHEALPix
template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::SubElt(const SphereHEALPix<T>& a)
{
  if (NbPixels() != a.NbPixels() )
    {
    throw(SzMismatchError("SphereHEALPix<T>::SubElt(const SphereHEALPix<T>&) SizeMismatch")) ;
    }
  pixels_ -= a.pixels_;
  return (*this);
}

//! Multiply two SphereHEALPix (elements by elements)
template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::MulElt(const SphereHEALPix<T>& a)
{
  if (NbPixels() != a.NbPixels() )
    {
    throw(SzMismatchError("SphereHEALPix<T>::MulElt(const SphereHEALPix<T>&) SizeMismatch")) ;
    }
  pixels_ *= a.pixels_;
  return (*this);
}

//! Divide two SphereHEALPix (elements by elements) - No protection for divide by 0
template <class T>
SphereHEALPix<T>& SphereHEALPix<T>::DivElt(const SphereHEALPix<T>& a)
{
  if (NbPixels() != a.NbPixels() )
    {
    throw(SzMismatchError("SphereHEALPix<T>::DivElt(const SphereHEALPix<T>&) SizeMismatch")) ;
    }
  pixels_ /= a.pixels_;
  return (*this);
}




template <class T>
void SphereHEALPix<T>::print(ostream& os) const
{
  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
  //
  os << " nSide_ = " << nSide_  << endl;
  os << " nPix_   = " << nPix_   << endl;
  os << " omeg_ = " << omeg_   << endl;

  os << " content of pixels : ";
  for(int i=0; i < nPix_; i++)
    {
      if(i%5 == 0) os << endl;
      os <<  pixels_(i) <<", ";
    }
  os << endl;


}



//*******************************************************************

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