source: Sophya/trunk/SophyaLib/SkyMap/spherehealpix.h@ 975

#ifndef SPHEREHEALPIX_SEEN
#define SPHEREHEALPIX_SEEN

#include "sphericalmap.h"
#include "tvector.h"
#include "ndatablock.h"

#include "anydataobj.h"
#include "ppersist.h"


namespace SOPHYA {


// ***************** CLASSE SphereHEALPix *****************************

//!  class SphereHEALPix
/*!
   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
\verbatim
        .   .   .   .   <--- 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)

*/


template<class T>
class FIO_SphereHEALPix;

template<class T>
class FITS_SphereHEALPix;

template<class T>
class SphereHEALPix : public SphericalMap<T>
{


  friend class FIO_SphereHEALPix<T>;
  friend class FITS_SphereHEALPix<T>;

public :

SphereHEALPix();
/*!    
  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)
*/
SphereHEALPix(int_4 m);
SphereHEALPix(const SphereHEALPix<T>& s, bool share);
SphereHEALPix(const SphereHEALPix<T>& s);
//!     Destructor
virtual ~SphereHEALPix();

  // Temporaire?
inline virtual bool   IsTemp(void) const {

    if (sliceBeginIndex_.IsTemp() != pixels_.IsTemp() || sliceLenght_.IsTemp() != pixels_.IsTemp() )
      throw PException(" l'etat 'temporaire' de la spherehealpix est incoherent"); 
    return pixels_.IsTemp();
}
/*! Setting blockdata to temporary (see ndatablock documentation) */
inline virtual void SetTemp(bool temp=false) const 
  {
    pixels_.SetTemp(temp);
    sliceBeginIndex_.SetTemp(temp);
    sliceLenght_.SetTemp(temp);
  };
// ------------------ Definition of PixelMap abstract methods

/* Nombre de pixels du decoupage */
/*!    Return number of  pixels of the  splitting */
virtual int_4 NbPixels() const;

/* Valeur du contenu du pixel d'indice "RING" k  */
/*!    Return value of  pixel with "RING" index k */
virtual T& PixVal(int_4 k);
virtual T const& PixVal(int_4 k) const;

/* Nombre de tranches en theta */ 
/*!    Return number of slices in theta direction on the  sphere */
uint_4 NbThetaSlices() 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 
*/
virtual void 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" 

   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 
*/
virtual void GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const ;

/* Return true if teta,phi in map  */
virtual bool ContainsSph(double theta, double phi) const;
/* Indice "RING" du pixel vers lequel pointe une direction definie par 
ses  coordonnees spheriques */                                    
/*!  Return "RING" index of the pixel corresponding to direction (theta, phi).
 */
virtual int_4 PixIndexSph(double theta,double phi) const;

/* Coordonnees spheriques du milieu du pixel d'indice "RING" k   */
virtual void PixThetaPhi(int_4 k,double& theta,double& phi) const;

/*! Set all pixels to value v */
virtual T SetPixels(T v);

/* Pixel Solid angle  (steradians) */
/*!    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.
*/
virtual double PixSolAngle(int_4 dummy=0) const;

/*  Acces to the DataBlock  */
inline       NDataBlock<T>& DataBlock()       {return pixels_;}
inline const NDataBlock<T>& DataBlock() const {return pixels_;}

// --------------- Specific methods 

/*!    
  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)
*/
virtual void Resize(int_4 m);

// pour l'instant le tableau est ordonne selon RING, uniquement
inline virtual char* TypeOfMap() const {return "RING";};


/* Valeur du contenu du pixel d'indice "NEST" k */
/*!    Return value of  pixel with "NESTED" index k */
virtual T& PixValNest(int_4 k);
/*!    Return value of  pixel with "NESTED" index k */
virtual T const& PixValNest(int_4 k) const;

/* Indice "NEST" du pixel vers lequel pointe une direction definie par 
ses  coordonnees spheriques */                                    
/*! Return "NESTED" index of the pixel corresponding to direction (theta, phi).
 */
virtual int_4 PixIndexSphNest(double theta,double phi) const;

/* Coordonnees spheriques du milieu du pixel d'indice "NEST" k       */
/*!   Return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
 */
virtual void PixThetaPhiNest(int_4 k,double& theta,double& phi) const;

/* algorithme de pixelisation */
void Pixelize(int_4); 

/* convertit index nested en ring  */
/*!    translation from NESTED index  into RING index */
int_4 NestToRing(int_4) const;

/* convertit index ring en nested" */
/*!    translation from  RING index  into NESTED index */
int_4 RingToNest(int_4) const;


/* retourne la valeur du parametre Gorski */
inline virtual int_4 SizeIndex() const {return(nSide_);}

/* impression */ 
void print(ostream& os) const;



inline  SphereHEALPix<T>& operator = (const SphereHEALPix<T>& a) 
                                                       {return Set(a);}
        
private :

// ------------- méthodes internes ----------------------
void InitNul();
void SetThetaSlices(); 

int_4  nest2ring(int_4 nside,int_4 ipnest) const;
int_4  ring2nest(int_4 nside,int_4 ipring) const;

int_4  ang2pix_ring(int_4 nside,double theta,double phi) const;
int_4  ang2pix_nest(int_4 nside,double theta,double phi) const;
void pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) const;
void pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) const;
inline void setParameters(int_4 nside, int_4 nbpixels, double solangle) 
  {
    nSide_= nside;
    nPix_= nbpixels;
    omeg_= solangle;
  }

 void CloneOrShare(const SphereHEALPix<T>& a);
 SphereHEALPix<T>& Set(const SphereHEALPix<T>& a);

// ------------- variables internes -----------------------

int_4 nSide_;
int_4 nPix_;
double omeg_;

NDataBlock<T> pixels_;
NDataBlock<int> sliceBeginIndex_;
NDataBlock<int> sliceLenght_;

};


//////////////////////////////////////////////////////////////////////////
//
// ------------- Classe PIXELS_XY -----------------------
//
class PIXELS_XY
{

public :

static PIXELS_XY& instance();

NDataBlock<int_4> pix2x_;
NDataBlock<int_4> pix2y_;
NDataBlock<int_4> x2pix_;
NDataBlock<int_4> y2pix_;

private :

PIXELS_XY();
void mk_pix2xy();
void mk_xy2pix();
};

} // Fin du namespace

#endif