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

#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        SphereHEALPix
//
// include      SphereHealpix.h strutil.h
//
//    Pixelisation Gorski  
//
//
//|    -----------------------------------------------------------------------
//|     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
//|
//    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)
//--
//++
//
// Links        Parents
//
//    SphericalMap 
//-- 

/* --Methode-- */
//++
// Titre        Constructors
//--
//++

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

//--
{
  InitNul();
  //  SetTemp(false);
}

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

//    m is the "nside" of the Gorski algorithm
//
//    The total number of pixels will be Npix =  12*nside**2
//
//    nside MUST be a power of 2 (<= 8192)
//--
{

  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();
  //  SetTemp(false);
  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);
}

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_);

  // pas forcement a conserver, pas forcement a cet endroit (GLM)
  //  if (a.IsTemp() ) SetTemp(true);
}

////////////////////////// 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 ) - Array 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 Array ! ");
  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
//--

//++
template<class T>
void SphereHEALPix<T>::Resize(int_4 m)

//    m is the "nside" of the Gorski algorithm
//
//    The total number of pixels will be Npix =  12*nside**2
//
//    nside MUST be a power of 2 (<= 8192)
//--
{
  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-- */
//++
template<class T>
int_4 SphereHEALPix<T>::NbPixels() const

//    Retourne le nombre de pixels du découpage
//--
{  
  return(nPix_);
}

//++
template<class T>
uint_4 SphereHEALPix<T>::NbThetaSlices() const 

//    Return number of slices in theta direction on the  sphere
//--
{
  uint_4 nbslices = uint_4(4*nSide_-1);
  if (nSide_<=0) 
    {
      nbslices = 0;
      throw PException(" sphere not pixelized, NbSlice=0 ");
    }
  return nbslices;
}

//++
template<class T>
void SphereHEALPix<T>::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 
//
//--
{

  if (index<0 || index >= NbThetaSlices()) 
    {
      // THROW(out_of_range("SphereHEALPix::PIxVal Pixel index out of range")); 
      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;
}
//++
//++

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

//    For a theta-slice with index 'sliceIndex', return : 
//
//    the corresponding "theta" 
//    the corresponding "phi" for first pixel of the slice 
//
//    a vector containing the indices of the pixels of the slice
//    (equally distributed in phi)
//
//    a vector containing the corresponding values of pixels 
//
//--
{

  if (sliceIndex<0 || sliceIndex >= NbThetaSlices()) 
    {
      // THROW(out_of_range("SphereHEALPix::PIxVal Pixel index out of range")); 
      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;
    }
}

/* --Methode-- */
//++
template<class T>
T& SphereHEALPix<T>::PixVal(int_4 k)

//    Return value of  pixel with "RING" index k 
//--
{
  if((k < 0) || (k >= nPix_)) 
    {
      throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range");
    }
  return pixels_(k);
}

/* --Methode-- */
//++
template<class T>
T const& SphereHEALPix<T>::PixVal(int_4 k) const

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

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

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

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

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

/* --Methode-- */
//++
template<class T>
bool SphereHEALPix<T>::ContainsSph(double /*theta*/, double /*phi*/) const
//--
{
return(true);
}

/* --Methode-- */
//++
template<class T>
int_4 SphereHEALPix<T>::PixIndexSph(double theta,double phi) const

//    Return "RING" index of the pixel corresponding to 
//    direction (theta, phi).
//--
{
  return ang2pix_ring(nSide_,theta,phi);
}

//++
template<class T>
int_4 SphereHEALPix<T>::PixIndexSphNest(double theta,double  phi) const

//    Return "NESTED" index of the pixel corresponding to 
//    direction (theta, phi).
//--
{
  return ang2pix_nest(nSide_,theta,phi);
}


/* --Methode-- */
//++
template<class T>
void SphereHEALPix<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const

//   Return (theta,phi) coordinates of middle of  pixel with "RING" index k
//--
{
  pix2ang_ring(nSide_,k,theta,phi);
}

template <class T>
T SphereHEALPix<T>::SetPixels(T v)
{
pixels_.Reset(v);
return(v);
}

//++
template<class T>
double SphereHEALPix<T>::PixSolAngle(int_4 /*dummy*/) const
//    Pixel Solid angle  (steradians)
//    All the pixels have the same solid angle. The dummy argument is
//    for compatibility with eventual pixelizations which would not 
//    fulfil this requirement.
//--
{
  return omeg_;
}

//++
template<class T>
void SphereHEALPix<T>::PixThetaPhiNest(int_4 k,double& theta,double& phi) const

//   Return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
//--
{   
  pix2ang_nest(nSide_,k,theta,phi);
}

//++
template<class T>
int_4 SphereHEALPix<T>::NestToRing(int_4 k) const

//    translation from NESTED index  into RING index 
//
//--
{
  return  nest2ring(nSide_,k);
}

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


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;

  os << endl;
  //const PIXELS_XY& PXY= PIXELS_XY::instance();

  //os << endl;  os << " contenu des tableaux conversions "<<endl;
  //for(int i=0; i < 5; i++)
  //  {
  //   os<<PXY.pix2x_(i)<<", "<<PXY.pix2y_(i)<<", "<<PXY.x2pix_(i)<<", "<<PXY.y2pix_(i)<<endl;
  // }
  os << endl;

}



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

#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template SphereHEALPix<uint_2>
#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<r_8>;
template class SphereHEALPix<r_4>;
template class SphereHEALPix< complex<r_4> >;
template class SphereHEALPix< complex<r_8> >;
#endif