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

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

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

/*! Setting blockdata to temporary (see ndatablock documentation) */
inline virtual void SetTemp(bool temp=false) const {pixels_.SetTemp(temp);};


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

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

/* return/set the number of pixels */
/*!    Return number of pixels */
virtual int_4 NbPixels() const;   // D.Y. int change en int_4 rationalisation Mac
  
/* return the value of pixel number k */
/*!    Return value of pixel with index k */
virtual T& PixVal(int_4 k);
/*!   const version of previous method */
virtual T const& PixVal(int_4 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_4 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_4 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_4 k) const; 

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

/*!    Resize storage area for pixels */
void ReSize(int_4 nx, int_4 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_4 x0,int_4 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_4 SizeIndex() const {return(nPix_);}
inline int_4 Size_x() const {return nSzX_;}
inline int_4 XSize() const {return nSzX_;}
inline int_4 Size_y() const {return nSzY_;}
inline int_4 YSize() const {return nSzY_;}

inline void Origin(double& theta,double& phi,int_4& x0,int_4& 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_;}


/*  Acces to the DataBlock  */
inline       NDataBlock<T>& DataBlock()       {return pixels_;}
inline const NDataBlock<T>& DataBlock() const {return pixels_;}

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

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

void InitNul();
/*!    Return 2 indices corresponding to the pixel number k */
void Getij(int_4 k,int_4& i,int_4& 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_4 nSzX_;
int_4 nSzY_;
int_4 nPix_;
bool originFlag_;
bool extensFlag_;
int_4 x0_;
int_4 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