source: Sophya/trunk/SophyaLib/SkyMap/spherethetaphi.h@ 2782

#ifndef SPHERETHETAPHI_SEEN
#define SPHERETHETAPHI_SEEN

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

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

namespace SOPHYA {



template <class T>
class FIO_SphereThetaPhi;  

template<class T>
class FITS_SphereThetaPhi;


// ***************** Class SphereThetaPhi *****************************
/*!    sphere splitted with respect to theta, phi : each hemisphere is 
    splitted into (m-1) parallels (equator does not enter into account).
   This operation defines m slices, each of which is splitted into 
    equidistant meridians. This splitting is realized in such a way that 
    all pixels have the same area and are as square as possible.

    One begins with the hemisphere with positive z, starting from the pole
    toward the equator. The first pixel is the polar cap ; it is circular
    and centered on theta=0. 
*/
template <class T>
class SphereThetaPhi : public SphericalMap<T>
{

public :

SphereThetaPhi();
/*!    m is the number of slices in theta on an hemisphere (the polar cap
   forms the first slice).
    pet is a dummy parameter at the moment. 
*/
SphereThetaPhi(int_4 m);
SphereThetaPhi(const SphereThetaPhi<T>& s, bool share);
SphereThetaPhi(const SphereThetaPhi<T>& s);
virtual ~SphereThetaPhi();

  // Temporaire?
inline virtual bool IsTemp(void) const {

    if (  NPhi_.IsTemp()  != pixels_.IsTemp() || 
          TNphi_.IsTemp() != pixels_.IsTemp()|| 
          Theta_.IsTemp() != pixels_.IsTemp()    )
      throw PException(" l'etat 'temporaire' de la spherethetaphi est incoherent"); 
    return pixels_.IsTemp();
}
/*! Setting blockdata to temporary (see ndatablock documentation) */
inline virtual void SetTemp(bool temp=false) const 
  {
    NPhi_.SetTemp(temp);
    TNphi_.SetTemp(temp);
    Theta_.SetTemp(temp);
    pixels_.SetTemp(temp);
  };

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

/* retourne le nombre de pixels */ 
/*!    Return total number of pixels  */
virtual int_4 NbPixels() const;

/* retourne la valeur du pixel d'indice k */ 
/*!    Return value of pixel with index k */
virtual T&       PixVal(int_4 k);
virtual T const& PixVal(int_4 k) const;

/* Return true if teta,phi in map  */
virtual bool ContainsSph(double theta, double phi) const;
/* retourne l'indice du pixel a (theta,phi) */ 
/*    Return index of the pixel corresponding to direction (theta, phi). */
virtual int_4 PixIndexSph(double theta, double phi) const;

/* retourne les coordonnees Spheriques du centre du pixel d'indice k */ 
/*!   Return (theta,phi) coordinates of middle of  pixel with  index k */
virtual void PixThetaPhi(int_4 k, double& theta, double& phi) const;

/*! Setting pixel values to a constant */
virtual T SetPixels(T v);

/* retourne/fixe l'angle Solide de Pixel   (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;
 
/* retourne/fixe la valeur du parametre de decoupage m */ 
inline virtual int_4 SizeIndex() const { return( NTheta_); }

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

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

/*!   re-pixelize the sphere */
virtual void Resize(int_4 m);

inline virtual string TypeOfMap() const {return string("TETAFI");};

/* Valeurs de theta des paralleles et phi des meridiens limitant le pixel d'indice k */
/*   Return values of theta,phi which limit the pixel with  index k */
virtual void Limits(int_4 k,double& th1,double& th2,double& phi1,double& phi2);

/* Nombre de tranches en theta */ 
/*!    Return number of theta-slices on the sphere */
uint_4 NbThetaSlices() const;

/* Nombre de pixels en phi de la tranche d'indice kt */ 
int_4 NPhi(int_4 kt) const;

/* Renvoie dans t1,t2 les valeurs respectives de theta min et theta max  */
/* de la tranche d'indice kt  */ 
/*!    Return  theta values which limit the slice kt */
void Theta(int_4 kt, double& t1, double& t2);

/* Renvoie dans p1,p2 les valeurs phimin et phimax du pixel d'indice jp  */
/* dans la tranche d'indice kt  */ 
/*!   Return values of phi which limit the jp-th pixel of the kt-th slice */
void Phi(int_4 kt, int_4 jp, double& p1, double& p2);

/* Renvoie l'indice k du pixel d'indice jp dans la tranche d'indice kt   */
/*!   Return pixel index  with sequence index jp in the slice kt */
int_4 Index(int_4 kt, int_4 jp) const;

/* Indice kt de la tranche et indice jp du pixel d'indice k  */ 
/*!    Return indices kt (theta) and jp (phi) of  pixel with index k */
void ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp);

/*!    achieve the splitting into pixels (m has the same signification 
    as for the constructor)

    Each theta-slice of the north hemisphere will be spitted starting f
    from  phi=0 ...

    South hemisphere is scanned in the same direction according to phi
    and from equator to the pole (the pixel following the last one of
    the slice closest to the equator with z>0, is the pixel with lowest 
    phi of the slice closest of the equator with z<0).
*/
void Pixelize(int_4); 

/*!   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 
*/
virtual 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 
*/
virtual void GetThetaSlice(int_4 index, r_8& theta, r_8& phi0,TVector<int_4>& pixelIndices, TVector<T>& value) const ;


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



// Operations diverses  = , +=, ...


SphereThetaPhi<T>& Set(const SphereThetaPhi<T>& a);
inline  SphereThetaPhi<T>& operator = (const SphereThetaPhi<T>& a) 
                                                        {return Set(a);}

// A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)

  //! Fill SphereThetaPhi with all elements equal to \b x
virtual SphereThetaPhi<T>& SetT(T a); 
inline  SphereThetaPhi<T>& operator = (T a) {return SetT(a);}

//! Add \b x to all elements
virtual SphereThetaPhi<T>& Add(T a); 
inline  SphereThetaPhi<T>&  operator += (T x)  { return Add(x); }
//! Substract \b x to all elements
virtual SphereThetaPhi<T>& Sub(T a,bool fginv=false); 
inline   SphereThetaPhi<T>&  operator -= (T x)  { return Sub(x); }
//! Multiply all elements by \b x
virtual SphereThetaPhi<T>& Mul(T a); 
inline  SphereThetaPhi<T>&  operator *= (T x)  { return Mul(x); }
//! Divide all elements by \b x
virtual SphereThetaPhi<T>& Div(T a); 
inline  SphereThetaPhi<T>&  operator /= (T x)  { return Div(x); }

// A += -=  (ajoute, soustrait element par element les deux spheres )
  //! Operator SphereThetaPhi += SphereThetaPhi
  virtual SphereThetaPhi<T>&  AddElt(const SphereThetaPhi<T>& a);
  inline  SphereThetaPhi<T>&  operator += (const SphereThetaPhi<T>& a)  { return AddElt(a); }



  virtual SphereThetaPhi<T>&  SubElt(const SphereThetaPhi<T>& a);
  //! Operator SphereThetaPhi -= SphereThetaPhi
  inline  SphereThetaPhi<T>&  operator -= (const SphereThetaPhi<T>& a)  { return SubElt(a); }
// Multiplication, division element par element les deux SphereThetaPhi 
  virtual SphereThetaPhi<T>&  MulElt(const SphereThetaPhi<T>& a);
  inline  SphereThetaPhi<T>&  operator *= (const SphereThetaPhi<T>& a)  { return MulElt(a); }
  virtual SphereThetaPhi<T>&  DivElt(const SphereThetaPhi<T>& a);
  inline  SphereThetaPhi<T>&  operator /= (const SphereThetaPhi<T>& a)  { return DivElt(a); }


  void CloneOrShare(const SphereThetaPhi<T>& a);
  void Share(const SphereThetaPhi<T>& a);

  SphereThetaPhi<T>& CopyElt(const SphereThetaPhi<T>& a);






 // friend declaration for classes which handle persistence and FITS IO
  friend class FIO_SphereThetaPhi<T>;
  friend class FITS_SphereThetaPhi<T>;

protected :

// ------------- méthodes internes ----------------------          
void InitNul();
inline void setParameters( int nbThetaIndex, int nbpix, double omega) 
  {
    NPix_= nbpix;
    Omega_= omega;
    NTheta_= nbThetaIndex;
  }

// ------------- variables internes ---------------------
 int_4 NTheta_;   // nombre de tranches en theta, pour une demi-sphere
 int_4 NPix_;     // nombre total de pixels
 double Omega_; // angle solide constant pour chaque pixel
 NDataBlock<int_4>  NPhi_;  // tableau donnant, pour chaque bande en theta, 
                         //le nombre de pixels selon phi
 NDataBlock<int_4> TNphi_; // tableau donnant ke nombre de pixels cumule,
                           // jusqu'au debut de chaque tranche
 NDataBlock<r_8> Theta_;
 NDataBlock<T> pixels_;
};

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs A (+,-,*,/) (T) x
/*! \ingroup SkyMap \fn operator+(const SphereThetaPhi<T>&,T)
  \brief Operator SphereThetaPhi = SphereThetaPhi + constant */
template <class T> inline SphereThetaPhi<T> operator + (const SphereThetaPhi<T>& a, T b)
    {SphereThetaPhi<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Add(b); return result;}
/*! \ingroup SkyMap \fn operator+(T,const SphereThetaPhi<T>&)
  \brief Operator SphereThetaPhi = constant + SphereThetaPhi */
template <class T> inline SphereThetaPhi<T> operator + (T b,const SphereThetaPhi<T>& a)
    {SphereThetaPhi<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Add(b); return result;}


/*! \ingroup SphereThetaPhi\fn operator-(const SphereThetaPhi<T>&,T)
  \brief Operator SphereThetaPhi = SphereThetaPhi - constant */
template <class T> inline SphereThetaPhi<T> operator - (const SphereThetaPhi<T>& a, T b)
    {SphereThetaPhi<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Sub(b); return result;}

/*! \ingroup  \fn operator-(T,const SphereThetaPhi<T>&)
  \brief Operator SphereThetaPhi = constant - SphereThetaPhi */
template <class T> inline SphereThetaPhi<T> operator - (T b,const SphereThetaPhi<T>& a)
    {SphereThetaPhi<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Sub(b,true); return result;}

/*! \ingroup SkyMap \fn operator*(const SphereThetaPhi<T>&,T)
  \brief Operator SphereThetaPhi = SphereThetaPhi * constant */
template <class T> inline SphereThetaPhi<T> operator * (const SphereThetaPhi<T>& a, T b)
    {SphereThetaPhi<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Mul(b); return result;}

/*! \ingroup SkyMap \fn operator*(T,const SphereThetaPhi<T>&)
  \brief Operator SphereThetaPhi = constant * SphereThetaPhi */
template <class T> inline SphereThetaPhi<T> operator * (T b,const SphereThetaPhi<T>& a)
    {SphereThetaPhi<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Mul(b); return result;}

/*! \ingroup SkyMap \fn operator/(const SphereThetaPhi<T>&,T)
  \brief Operator SphereThetaPhi = SphereThetaPhi / constant */
template <class T> inline SphereThetaPhi<T> operator / (const SphereThetaPhi<T>& a, T b)
    {SphereThetaPhi<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Div(b); return result;}

/*! \ingroup SkyMap \fn operator/(T,const SphereThetaPhi<T>&)
  \brief Operator SphereThetaPhi = constant / SphereThetaPhi  */
template <class T> inline SphereThetaPhi<T> operator / (T b, const SphereThetaPhi<T>& a)
    {SphereThetaPhi<T> result; result.CloneOrShare(a); result.SetTemp(true); 
    result.Div(b, true); return result;}

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs C = A (+,-) B

/*! \ingroup SkyMap \fn operator+(const SphereThetaPhi<T>&,const SphereThetaPhi<T>&)
  \brief Operator SphereThetaPhi = SphereThetaPhi + SphereThetaPhi */
template <class T>
inline SphereThetaPhi<T> operator + (const SphereThetaPhi<T>& a,const SphereThetaPhi<T>& b)
    { SphereThetaPhi<T> result; result.SetTemp(true); 
    if (b.IsTemp())  { result.Share(b); result.AddElt(a); }
    else { result.CloneOrShare(a); result.AddElt(b); }
    return result; }

/*! \ingroup SkyMap \fn operator-(const SphereThetaPhi<T>&,const SphereThetaPhi<T>&)
  \brief Operator SphereThetaPhi = SphereThetaPhi - SphereThetaPhi */
template <class T>
inline SphereThetaPhi<T> operator - (const SphereThetaPhi<T>& a,const SphereThetaPhi<T>& b)
    { SphereThetaPhi<T> result; result.SetTemp(true); 
    if (b.IsTemp())  { result.Share(b); result.SubElt(a, true); }
    else { result.CloneOrShare(a); result.SubElt(b); }
    return result; }

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs C = A (*,/) B

/*! \ingroup SkyMap \fn operator*(const SphereThetaPhi<T>&,const SphereThetaPhi<T>&)
  \brief Operator SphereThetaPhi = SphereThetaPhi * SphereThetaPhi (pixel by pixel multiply)*/
template <class T>
inline SphereThetaPhi<T> operator * (const SphereThetaPhi<T>& a,const SphereThetaPhi<T>& b)
    { SphereThetaPhi<T> result; result.SetTemp(true); 
    if (b.IsTemp())  { result.Share(b); result.MulElt(a); }
    else { result.CloneOrShare(a); result.MulElt(b); }
    return result; }

/*! \ingroup SkyMap \fn operator/(const SphereThetaPhi<T>&,const SphereThetaPhi<T>&)
  \brief Operator SphereThetaPhi = SphereThetaPhi / SphereThetaPhi (pixel by pixel divide) */
template <class T>
inline SphereThetaPhi<T> operator / (const SphereThetaPhi<T>& a,const SphereThetaPhi<T>& b)
    { SphereThetaPhi<T> result; result.SetTemp(true); 
    result.CloneOrShare(a); result.DivElt(b); 
    return result; }


} // Fin du namespace

#endif