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

// This may look like C code, but it is really -*- C++ -*-
#ifndef LOCALMAP_SEEN
#define LOCALMAP_SEEN

#include "pixelmap.h"
#include "sphericalmap.h"
#include "ndatablock.h"

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

//! A local map of a region of the sky, in cartesian coordinates.
/*! A local map has an origin in (theta0, phi0), mapped to pixel(x0, y0)
   (x0, y0 might be outside of this local map)
 default value of (x0, y0) is middle of the map, center of pixel(nx/2, ny/2)
    A local map is a 2 dimensional array, with i as column index and j
    as row index. The map is supposed to lie on a plan tangent to the 
    celestial sphere in a point whose coordinates are (x0,y0) on the local
   map and (theta0, phi0) on the sphere. The range of the map is defined
    by two values of angles covered respectively by all the pixels in 
    x direction and all the pixels in y direction (SetSize()).

    A "reference plane" is considered : this plane is tangent to the 
    celestial sphere in a point with angles theta=Pi/2 and phi=0. This 
    point is the origine of coordinates is of the reference plane. The 
    x-axis is the tangent parallel to the equatorial line and oriented 
    toward the increasing phi's ; the y-axis is parallel to the meridian 
    line and oriented toward the north pole.
    
    Internally, a map is first defined within this reference plane and
    tranported until the point (theta0, phi0) in such a way that both 
    axes are kept parallel to meridian and parallel lines of the sphere.
    The user can define its own map with axes rotated with respect to
    reference axes (this rotation is characterized by angle between 
    the local parallel line and the wanted x-axis-- see method 
    SetOrigin(...))
*/
//
//    la carte est consideree comme un tableau a deux indices i et j, i etant 
//    indice de colonne et j indice de ligne. La carte est supposee resider 
//    dans un plan tangent, dont le point de tangence est repere (x0,y0) dans 
//    la carte et (theta0, phi0) sur la sphere celeste. L extension de la 
//    carte est definie par les valeurs de deux angles couverts respectivement
//    par la totalite des pixels en x de la carte et la totalite des pixels 
//    en y. (SetSize()).
//    On considere un "plan de reference" : plan tangent a la sphere celeste 
//    aux angles theta=Pi/2 et phi=0. Dans ce plan L origine des coordonnees 
//    est le point de tangence. L axe Ox est la tangente parallele a 
//    lequateur, dirige vers les phi croissants, l axe Oy est parallele 
//    au meridien, dirige vers le pole nord.
//    De maniere interne a la classe une carte est definie dans ce plan de 
//    reference et transportee  jusqu au point (theta0, phi0) de sorte que les //    axes restent paralleles aux meridiens et paralleles. L utilisateur peut 
//    definir sa carte selon un repere en rotation par rapport au repere de 
//    reference (par l angle entre le parallele et l axe Ox souhaite -- 
//    methode SetOrigin(...))


// ***************** Class LocalMap *****************************

template<class T>
class LocalMap : public PixelMap<T>, public AnyDataObj 
{

public:

LocalMap();
LocalMap(int nx, int ny);
LocalMap(const LocalMap<T>& lm, bool share=false);
virtual ~LocalMap();

// ---------- Overloading of () to access pixel number k ----

inline T& operator()(int k) {return(PixVal(k));}
inline T const& operator()(int k) const {return(PixVal(k));}
inline T& operator()(int ix, int iy) {return PixVal(iy*nSzX_+ix);};
inline T const& operator()(int ix, int iy) const {return PixVal(iy*nSzX_+ix);};
   
// ---------- Definition of PixelMap abstract methods -------

/* return/set the number of pixels */
/*!    Return number of pixels */
virtual int NbPixels() const;
inline void setNbPixels(int n) {nPix_= n;}
  
/* return the value of pixel number k */
/*!    Return value of pixel with index k */
virtual T& PixVal(int k);
/*!   const version of previous method */
virtual T const& PixVal(int k) const;

/* Return true if teta,phi in map  */
virtual bool ContainsSph(double theta, double phi) const;
/* return the index of pixel at (theta,phi) */
/*!    Return index of the pixel with spherical coordinates (theta,phi) */
virtual int PixIndexSph(double theta,double phi) const;
   
/* return the spherical coordinates of center of pixel number k */
/*!    Return (theta, phi) coordinates of pixel with index k */
virtual void PixThetaPhi(int k,double& theta,double& phi) const;

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

/* return the Pixel Solid angle  (steradians) */
/*!    Pixel Solid angle  (steradians)

    All the pixels have not necessarly the same size in (theta, phi)
    because of the projection scheme which is not yet fixed. 
*/  
virtual double PixSolAngle(int k) const; 

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

/*!    Resize storage area for pixels */
void ReSize(int nx, int ny);

inline virtual char* TypeOfMap() const {return "LOCAL";};
 
/* Origin (with angle between x axis and phi axis, in degrees)  x0,y0  the default: middle of map*/
/*!    set the referential of the map (angles in degrees)

    (default x0=siz_x/2,  y0=siz_y/2) 
*/
virtual void SetOrigin(double theta=90.,double phi=0.,double angle=0.);
/*!    set the referential of the map (angles in degrees) */
virtual void SetOrigin(double theta,double phi,int x0,int y0,double angle=0.);

/* Pixel size (degres) */ 
/*!    angle range of tthe map (angles in degrees) */
virtual void SetSize(double angleX,double angleY);

/* Check to see if the local mapping is done */
inline bool LocalMap_isDone() const {return(originFlag_ && extensFlag_);};

/*! Projection to a spherical map */
virtual void Project(SphericalMap<T>& sphere) const;
  
/* There should be a more complex algorithm somewhere to combine *several* local maps to a full sphere.
      -> static method, or separate class */
  
/* provides a integer characterizing the pixelization refinement  (here : number of pixels) */
inline virtual int SizeIndex() const {return(nPix_);}
inline int Size_x() const {return nSzX_;}
inline void setSize_x(int n) {nSzX_= n;}
inline int Size_y() const {return nSzY_;}
inline void setSize_y(int n) {nSzY_= n;}

inline void Origin(double& theta,double& phi,int& x0,int& y0,double& angle) const {theta= theta0_; phi= phi0_; x0= x0_; y0= y0_;angle= angle_;}

inline void Aperture(double& anglex,double& angley) const {anglex= angleX_; angley= angleY_;}

/* retourne le pointeur vers/remplit  le vecteur des contenus des pixels */ 
inline const NDataBlock<T>* getDataBlock() const {return (&pixels_);}
inline void setDataBlock(T* data, int n) {pixels_.FillFrom(n,data);}

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

// ---------- Méthodes internes -----------------------------
          
private :

void InitNul();
/*!    Return 2 indices corresponding to the pixel number k */
void Getij(int k,int& i,int& j) const;
/*!    Transform a pair of coordinates (theta, phi) given in 
    reference coordinates into map coordinates
*/
void ReferenceToUser(double& theta,double& phi) const; 
/*!    Transform a pair of coordinates (theta, phi) given in 
   map coordinates into reference coordinates
*/
void UserToReference(double& theta,double& phi) const;
/*!   Given coordinates in pixel units in the REFERENCE PLANE, return
    (theta, phi) in "absolute" referential theta=pi/2 ,phi=0.   
*/
void PixProjToAngle(double x,double y,double& theta,double& phi) const; 
/*!    Given coordinates  (theta, phi) in "absolute" referential 
    theta=pi/2 ,phi=0  return pixel indices  (i,j) in the REFERENCE PLANE.
*/
void AngleProjToPix(double theta,double phi,double& x,double& y) const; 

// ---------- Variables internes ----------------------------

int nSzX_;
int nSzY_;
int nPix_;
bool originFlag_;
bool extensFlag_;
int x0_;
int y0_;
double theta0_;
double phi0_;
double angle_;
double cos_angle_;
double sin_angle_;
double angleX_;
double angleY_;
double tgAngleX_;
double tgAngleY_;
NDataBlock<T> pixels_;
};

// ------------- Classe pour la gestion de persistance --
//!    Delegated objects for persitance management
template <class T>
class FIO_LocalMap : public PPersist  
{

public:

FIO_LocalMap();
FIO_LocalMap(string const & filename); 
FIO_LocalMap(const LocalMap<T>& obj);
FIO_LocalMap(LocalMap<T>* obj);
virtual ~FIO_LocalMap();
virtual AnyDataObj* DataObj();
inline operator LocalMap<T>() { return(*dobj); }
//inline LocalMap<T> getObj() { return(*dobj); }

protected :

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

#endif