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

#ifndef SPHEREGORSKI_SEEN
#define SPHEREGORSKI_SEEN

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

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

// les attributs de classe pix2x_ et pix2y_ sont relatifs a la traduction des 
// indices  RING en indices NESTED (ou l inverse). Ce sont des tableaux 
// d'entiers initialises et remplis par le constructeur. Dans la version 
// FORTRAN de healpix, il s'agissait de tableaux mis en COMMON. Ils etaient 
// initialises au premier appel d'une conversion d indice. Je ne peux pas 
// garder cette procedure car, par exemple, la fonction PixValNest() const 
// n'autorisera pas la constitution de ces tableaux (a cause du const). Ainsi, 
// des la creation d un objet SphereGorski ces tableaux sont calcules, ce qui 
// est, certes, inutile si on ne se sert pas des indices NESTED. Ca ne me 
// rejouit pas, mais c est le seul moyen que j ai trouve pour temir compte de 
// toutes les demandes et des options prises.
// 
//                                                G. Le Meur 

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

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

public :

SphereGorski();
SphereGorski(int m);
SphereGorski(const SphereGorski<T>& s);
virtual ~SphereGorski();

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

/* Nombre de pixels du decoupage */
virtual int NbPixels() const;
inline void setNbPixels(int n) {nPix_= n;}

/* Valeur du contenu du pixel d'indice "RING" k  */
virtual T& PixVal(int k);
virtual T const& PixVal(int k) const;

/* Nombre de tranches en theta */ 
int NbThetaSlices() const;
void GetThetaSlice(int index,double& theta,TVector<double>& phi,TVector<T>& value) const;

/* Indice "RING" du pixel vers lequel pointe une direction definie par 
ses  coordonnees spheriques */                                    
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;

/* Pixel Solid angle  (steradians) */
virtual double PixSolAngle(int dummy) const;
inline void setPixSolAngle(double x) {omeg_= x;}

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

virtual void Resize(int m);

inline virtual char* TypeOfMap() const {return "RING";};

/* Valeur du contenu du pixel d'indice "NEST" k */
virtual T& PixValNest(int k);
virtual T const& PixValNest(int k) const;

/* Indice "NEST" du pixel vers lequel pointe une direction definie par 
ses  coordonnees spheriques */                                    
virtual int PixIndexSphNest(double theta,double phi) const;

/* Coordonnees spheriques du milieu du pixel d'indice "NEST" k       */
virtual void PixThetaPhiNest(int k,double& theta,double& phi) const;

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

/* convertit index nested en ring  */
int NestToRing(int) const;

/* convertit index ring en nested" */
int RingToNest(int) const;

/* analyse en harmoniques spheriques des valeurs des pixels de la 
   sphere : appel du module anafast (Gorski-Hivon) */
//void anharm(int, float, float*);
 
/*synthese  des valeurs des pixels de la sphere par l'intermediaire 
  des coefficients en harmoniques spheriques reconstitues apartir d'un 
  spectre en puissance : appel du module synfast (Gorski-Hivon) */
//void synharm(int, int, float, float*); 

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

/* retourne/fixe les parametres des modules anafast et synfast */
void getParafast(int& nlmax,int& nmmax,int& iseed,float& fwhm,float& quadr,float& cut) const;
void setParafast(int nlmax,int nmmax,int iseed,float fwhm,float quadr,float cut,char* filename);

/* retourne/fixe le nom du fichier qui contient le spectre de puissance */
void powfile(char filename[]) const;

/* 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_;

int nlmax_;
int nmmax_;
int iseed_;
float fwhm_;
float quadrupole_;
float sym_cut_deg_;
char powFile_[128];
};

//
// ------------- 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