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

#ifndef SPHEREHEALPIX_SEEN
#define SPHEREHEALPIX_SEEN

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

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

#include "HEALPixUtils.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>
{
public :
// Static Methods to ease the use of HEALPix index <> angle conversion methods

  static inline int_4  nest2ring(int_4 nside,int_4 ipnest) 
         { return HEALPix::nest2ring(nside, ipnest); }
  static inline int_4  ring2nest(int_4 nside,int_4 ipring) 
         { return HEALPix::ring2nest(nside, ipring); }
  static inline int_4  ang2pix_ring(int_4 nside,double theta,double phi) 
         { return HEALPix::ang2pix_ring(nside, theta, phi); }
  static inline int_4  ang2pix_nest(int_4 nside,double theta,double phi) 
         { return HEALPix::ang2pix_nest(nside, theta, phi); }
  static inline void   pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi) 
         { HEALPix::pix2ang_ring(nside, ipix, theta, phi); }
  static inline void   pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi) 
         { HEALPix::pix2ang_nest(nside, ipix, theta, phi); }

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

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


 // friend declaration for classes which handle persistence and FITS IO
 friend class FIO_SphereHEALPix<T>;
 friend class FITS_SphereHEALPix<T>;
       
protected :

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

inline void setParameters(int_4 nside, int_4 nbpixels, double solangle) 
  {
    nSide_= nside;
    nPix_= nbpixels;
    omeg_= solangle;
  }

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

int_4 nSide_;
int_4 nPix_;
double omeg_;

NDataBlock<T> pixels_;
NDataBlock<int_4> sliceBeginIndex_;             // Rationalisation Mac.         D.Y.
NDataBlock<int_4> sliceLenght_;

};



} // Fin du namespace

#endif