source: Sophya/trunk/SophyaLib/Samba/alm.h@ 729

#ifndef ALM_SEEN
#define ALM_SEEN
#include "nbrandom.h"
#include "nbmath.h"
#include "tvector.h"
#include "pexceptions.h"
/*! Class for inferior triangular matrix (base class for the class Alm) */
template <class T>
class TriangularMatrix 
  {
    
  public :

TriangularMatrix()   {};
  /* instanciate a triangular matrix from the number of rows */
TriangularMatrix(int rowSize)  : long_diag_((uint_4)rowSize) {elem_.ReSize((uint_4) (rowSize*(rowSize+1)/2) ); };
TriangularMatrix(const TriangularMatrix<T>& a,  bool share=false)  : elem_(a.elem_, share),  long_diag_(a.long_diag_) {;}
/*! resize the matrix with a new number of rows */
inline void ReSizeRow(int rowSize) 
  {
    long_diag_=(uint_4)rowSize;
    elem_.ReSize(long_diag_*(long_diag_+1)/2);
  }
 inline void SetTemp(bool temp=false) const {elem_.SetTemp(temp);}

inline TriangularMatrix<T>& operator = (const TriangularMatrix<T>& a)
{
  elem_=a.elem_;
  long_diag_ = a.long_diag_;
  return *this;
}
inline T& operator()(int l, int m) 
  {
      return  elem_(adr_ij(l,m));
  }
inline T const& operator()(int l, int m) const 
  {
      return *(elem_.Begin()+ adr_ij(l,m));
  }
 inline  int_4  rowNumber() const {return (int_4)long_diag_;}
  private: 
 /*! compute the address of an element in the single array representing the matrix */
inline uint_4 adr_ij(int i,int j) const 
{
  int adr= i*(i+1)/2+j;
  if ( adr >= elem_.Size() || adr <0 )
{
 cout << " attention depassement dans triangularMatrix " << endl;
 cout << " l= " << i << " m= " << j << " tableau reserve longueur " << elem_.Size() << endl; 
}
  return(i*(i+1)/2+j);
}

 uint_4 long_diag_;
 NDataBlock<T> elem_;

  };


/*! class for the coefficients \f$a_{lm}\f$ of the development of a function defined on a sphere, in spherical harmonics */
template <class T>
class Alm : public TriangularMatrix<complex<T> >
  {
    public :
 
Alm() : TriangularMatrix<complex<T> >()  {};
    /* instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s */
Alm(int nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;}
Alm(const Alm<T>& a,  bool share=false)  : TriangularMatrix<complex<T> >(a, share)  {;}

Alm(const TVector<T>& clin, const r_8 fwhm) ;
/*! resize with a new lmax */
inline void ReSizeToLmax(int_4 nlmax) {ReSizeRow(nlmax+1);}
inline int_4 Lmax() const {return rowNumber()-1;}
TVector<T> powerSpectrum() const;


};
/*! class for a vector with an index running from \f$-m_{max}\f$ to \f$+m_{max}\f$ (then the size of the vector will be actually \f$2m_{max}+1)\f$ */
template <class T>
class Bm
  {
    public :
Bm(): nmmax_(0) {;};
    /* instanciate from the maximum absolute value of the index */
Bm(int mmax) : nmmax_((uint_4)mmax)   { bm_.ReSize( (uint_4)(2*mmax+1) );}
 Bm(const Bm<T>& b,  bool share=false) : nmmax_(b.nmmax_), bm_(b.bm_, share) {;}
 /*! resize with a new value of mmax */
inline void ReSizeToMmax(int_4 mmax) 
  {
    nmmax_=(uint_4)l;
    bm_.ReSize(2* nmmax_+1);
  }
//inline int_4 Size() const {return 2*nmmax_+1;};
inline T& operator()(int m) {return bm_( adr_i(m));};
inline T const& operator()(int m) const {return *(bm_.Begin()+ adr_i(m));};
/*! return the current value of the maximum absolute value of the index */
inline int_4 Mmax() const 
   {
     return (int_4)nmmax_;
   }
inline Bm<T>& operator = (const Bm<T>& b)
 {
   if (this != &b)
     {
       nmmax_= b.nmmax_;
       bm_=  b.bm_;
     }
   return *this;
 }

  private:
/*! return the address in the array representing the vector */
inline uint_4 adr_i(int i) const 
{
  return(uint_4)(i+nmmax_);
}



 uint_4 nmmax_;
 NDataBlock<T> bm_;

  };


#endif