source: Sophya/trunk/SophyaLib/Samba/spheregorski.h@ 574

#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 SphereGorski : public SphericalMap<T>,  public AnyDataObj 
{

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

// ------------------ Definition of PixelMap abstract methods

/* Nombre de pixels du decoupage */
/*!    Return number of  pixels of the  splitting */
virtual int NbPixels() const;
inline void setNbPixels(int n) {nPix_= n;}

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

/* Nombre de tranches en theta */ 
/*!    Return number of slices in theta direction on the  sphere */
int 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 index,double& theta,TVector<double>& phi,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 PixIndexSph(double theta,double phi) const;

/* Coordonnees spheriques du milieu du pixel d'indice "RING" k   */
virtual void PixThetaPhi(int 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 dummy) const;
inline void setPixSolAngle(double x) {omeg_= x;}

// --------------- 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 m);

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 k);
/*!    Return value of  pixel with "NESTED" index k */
virtual T const& PixValNest(int 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 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 k,double& theta,double& phi) const;

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

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

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


/* retourne/fixe la valeur du parametre Gorski */
inline virtual int SizeIndex() const {return(nSide_);}
inline void setSizeIndex(int n) {nSide_= n;}

/* retourne les pointeurs /remplit les tableaux */
inline const NDataBlock<T>* getDataBlock() const { return (&pixels_); }
inline void setDataBlock(T* data,int m) { pixels_.FillFrom(m,data); }



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

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

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

int  ang2pix_ring(int nside,double theta,double phi) const;
int  ang2pix_nest(int nside,double theta,double phi) const;
void pix2ang_ring(int nside,int ipix,double& theta,double& phi) const;
void pix2ang_nest(int nside,int ipix,double& theta,double& phi) const;

// ------------- variables internes -----------------------
int nSide_;
int nPix_;
double omeg_;

NDataBlock<T> pixels_;

};

//
// ------------- 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();
inline operator SphereGorski<T>() { return(*dobj); }
//inline SphereGorski<T> getObj() { 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> pix2x_;
NDataBlock<int> pix2y_;
NDataBlock<int> x2pix_;
NDataBlock<int> y2pix_;

private :

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