source: Sophya/trunk/SophyaLib/SkyMap/spheregorski.h@ 808

#ifndef SPHEREGORSKI_SEEN
#define SPHEREGORSKI_SEEN

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

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


// ***************** CLASSE SphereGorski *****************************

//!  class SphereGorski
/*!
   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_SphereGorski;

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

  friend class FIO_SphereGorski<T>;

public :

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

/*! Setting blockdata to temporary (see ndatablock documentation) */
inline virtual void SetTemp(bool temp=false) const {pixels_.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 
*/
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 
*/
 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_);}

/* retourne les pointeurs /remplit les tableaux */




/* impression */ 
void print(ostream& os) const;
         
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;
  }

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

int_4 nSide_;
int_4 nPix_;
double omeg_;

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

};

//
// ------------- Classe pour la gestion de persistance --
//
template <class T>
class FIO_SphereGorski : public PPersist  
{
public:

FIO_SphereGorski();
FIO_SphereGorski(string const & filename); 
FIO_SphereGorski(const SphereGorski<T>& obj);
FIO_SphereGorski(SphereGorski<T>* obj);
virtual ~FIO_SphereGorski();
virtual AnyDataObj* DataObj();
virtual void        SetDataObj(AnyDataObj & o);
inline operator SphereGorski<T>() { return(*dobj); }

protected :

virtual void ReadSelf(PInPersist&);           
virtual void WriteSelf(POutPersist&) const;  
SphereGorski<T>* dobj;
bool ownobj;
};

//
// ------------- 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();
};
#endif