source: Sophya/trunk/SophyaLib/SkyMap/spherethetaphi.cc@ 980

#include "spherethetaphi.h"
#include "smathconst.h"
#include <complex>
#include "piocmplx.h"
#include "fiondblock.h"
#include <iostream.h>


//***************************************************************
//++ 
// Class        SphereThetaPhi
// 
//
// include      spherethetaphi.h 
//
//    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. 
//--
//++
//
// Links        Parents
//
//    SphericalMap 
//-- 

/* --Methode-- */
//++
// Titre        Constructors
//--
//++

template <class T>
SphereThetaPhi<T>::SphereThetaPhi()
  : NPhi_(), TNphi_(), Theta_(), pixels_() 

//--
{
  InitNul();
}


/* --Methode-- */

//++
template <class T>
SphereThetaPhi<T>::SphereThetaPhi(int_4 m)

//    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. 
//--
{
  InitNul();
  Pixelize(m);
}

template <class T>
SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s, bool share)
  : NPhi_(s.NPhi_, share), TNphi_(s.TNphi_, share), Theta_(s.Theta_, share),
    pixels_(s.pixels_ , share) 
{

  NTheta_= s.NTheta_;
  NPix_  = s.NPix_;
  Omega_ = s.Omega_;
  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
}

template <class T>
SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)
  : NPhi_(s.NPhi_), TNphi_(s.TNphi_), Theta_(s.Theta_), pixels_(s.pixels_) 
{

  NTheta_= s.NTheta_;
  NPix_  = s.NPix_;
  Omega_ = s.Omega_;
  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
}

//++
// Titre        Destructor
//--
//++
template <class T>
SphereThetaPhi<T>::~SphereThetaPhi()

//--
{}

//++
// Titre        Public Méthods
//--
template <class T>
void SphereThetaPhi<T>::InitNul()
//
{
  NTheta_= 0;
  NPix_  = 0;
//  pixels_.Reset();  Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
}


//++
template <class T>
void SphereThetaPhi<T>::Resize(int_4 m)
//   re-pixelize the sphere
//--
{
  InitNul();
  Pixelize(m);
}

template<class T>
void  SphereThetaPhi<T>::CloneOrShare(const  SphereThetaPhi<T>& a)
{

  NTheta_= a.NTheta_;
  NPix_  = a.NPix_;
  Omega_ = a.Omega_;
  NPhi_.CloneOrShare(a.NPhi_);
  TNphi_.CloneOrShare(a.TNphi_);
  Theta_.CloneOrShare(a.Theta_);
  pixels_.CloneOrShare(a.pixels_);
}

////////////////////////// methodes de copie/share
template<class T>
SphereThetaPhi<T>& SphereThetaPhi<T>::Set(const SphereThetaPhi<T>& a)
{
  if (this != &a) 
    { 
      
      
      if (a.NbPixels() < 1) 
        throw RangeCheckError("SphereThetaPhi<T>::Set(a ) - Array a not allocated ! ");
      if (NbPixels() < 1) CloneOrShare(a);
      else CopyElt(a);
      if (mInfo_) delete mInfo_;
      mInfo_ = NULL;
      if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
    }
  return(*this);
}

template<class T>
SphereThetaPhi<T>& SphereThetaPhi<T>::CopyElt(const SphereThetaPhi<T>& a)
{
  if (NbPixels() < 1) 
    throw RangeCheckError("SphereThetaPhi<T>::CopyElt(const SphereThetaPhi<T>& )  - Not Allocated Array ! ");
  if (NbPixels() != a.NbPixels()) 
    throw(SzMismatchError("SphereThetaPhi<T>::CopyElt(const SphereThetaPhi<T>&) SizeMismatch")) ;

  NTheta_= a.NTheta_;
  NPix_  = a.NPix_;
  Omega_ = a.Omega_;
  int k;
  for (k=0; k< NPix_; k++) pixels_(k) = a.pixels_(k);
  for (k=0; k< a.NPhi_.Size(); k++) NPhi_(k) = a.NPhi_(k);
  for (k=0; k< a.TNphi_.Size(); k++) TNphi_(k) = a.TNphi_(k);
  for (k=0; k< a.Theta_.Size(); k++) Theta_(k) = a.Theta_(k);
  return(*this);

}

/* --Methode-- */
//++
template <class T>
int_4 SphereThetaPhi<T>::NbPixels() const

//    Return total number of pixels 
//--
{
  return(NPix_);
}

/* --Methode-- */
//++
template <class T>
T& SphereThetaPhi<T>::PixVal(int_4 k)

//    Return value of pixel with index k
//--
{
  if((k < 0) || (k >= NPix_)) 
    {
      throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
    }
  return pixels_(k);
}

//++
template <class T>
T const& SphereThetaPhi<T>::PixVal(int_4 k) const

//    Return value of pixel with index k
//--
{
  if((k < 0) || (k >= NPix_)) 
    {
       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
    } 
  return *(pixels_.Data()+k);
}

/* --Methode-- */
//++
template <class T>
bool SphereThetaPhi<T>::ContainsSph(double /*theta*/, double /*phi*/) const
//--
{
  return(true);
}

/* --Methode-- */
//++
template <class T>
int_4 SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const

//    Return index of the pixel corresponding to 
//    direction (theta, phi).
//--
{
  double dphi;
  int i,k;
  bool fgzn = false;

  if((theta > Pi) || (theta < 0.)) return(-1);
  if((phi < 0.) || (phi > DeuxPi)) return(-1);
  if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;}
  
  // La bande d'indice kt est limitée par les valeurs de theta contenues dans
  // Theta_[kt] et Theta_[kt+1]
  for( i=1; i< NTheta_; i++ ) 
    if( theta < Theta_(i) ) break;
  
  dphi= DeuxPi/(double)NPhi_(i-1);
  
  if (fgzn)  k= NPix_-TNphi_(i)+(int_4)(phi/dphi);
  else k= TNphi_(i-1)+(int_4)(phi/dphi); 
  return(k);
}

/* --Methode-- */
//++
template <class T>
void SphereThetaPhi<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const

//   Return (theta,phi) coordinates of middle of  pixel with  index k
//--
{
  int i; 
  bool fgzn = false;
  
  if((k < 0) || (k >= NPix_)) {theta = -99999.; phi = -99999.; return; }
  if( k >= NPix_/2)  {fgzn = true; k = NPix_-1-k;}

  // recupère l'indice i de la tranche qui contient le pixel k
  for( i=0; i< NTheta_; i++ ) 
    if( k < TNphi_(i+1) ) break;

  // angle theta
  theta= 0.5*(Theta_(i)+Theta_(i+1));
  if (fgzn) theta= Pi-theta;
  
  // angle phi
  k -= TNphi_(i);
  phi= DeuxPi/(double)NPhi_(i)*(double)(k+.5);
  if (fgzn) phi= DeuxPi-phi;
}

template <class T>
T SphereThetaPhi<T>::SetPixels(T v)
{
pixels_.Reset(v);
return(v);
}

//++
template <class T>
double SphereThetaPhi<T>::PixSolAngle(int_4 /*dummy*/) const

//    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.
//--
{
  return Omega_;
}

/* --Methode-- */
//++
template <class T>
void SphereThetaPhi<T>::Limits(int_4 k,double& tetMin,double& tetMax,double& phiMin,double& phiMax) 

//   Return values of theta,phi which limit the pixel with  index k
//--
  {
  int j;
  double dphi;
  bool fgzn= false;
  
  if((k < 0) || (k >= NPix_)) {
    tetMin= -99999.; 
    phiMin= -99999.;  
    tetMax= -99999.; 
    phiMax= -99999.;  
    return;
  }
  
  // si on se trouve dans l'hémisphère Sud
  if(k >= NPix_/2) { 
    fgzn= true; 
    k= NPix_-1-k; 
  }
  
  // on recupere l'indice i de la tranche qui contient le pixel k
  int i;
  for( i=0; i< NTheta_; i++ ) 
    if(k < TNphi_(i+1)) break;
  
  // valeurs limites de theta dans l'hemisphere Nord
  tetMin= Theta_(i);
  tetMax= Theta_(i+1);
  // valeurs limites de theta dans l'hemisphere Sud
  if (fgzn) {
    tetMin= Pi - Theta_(i+1);
    tetMax= Pi - Theta_(i);
  }
  
  // pixel correspondant dans l'hemisphere Nord
  if (fgzn) k= TNphi_(i+1)-k+TNphi_(i)-1;
  
  // indice j de discretisation ( phi= j*dphi )
  j= k-TNphi_(i);
  dphi= DeuxPi/(double)NPhi_(i);
  
  // valeurs limites de phi 
  phiMin= dphi*(double)(j);
  phiMax= dphi*(double)(j+1);
  return;
}

/* --Methode-- */
//++
template <class T>
uint_4 SphereThetaPhi<T>::NbThetaSlices() const

//    Return number of theta-slices on the sphere
//--
{
  if (NTheta_<=0) 
    {
      throw PException(" sphere not pixelized, NbSlice=0 ");
    }
  return( 2*NTheta_);
}

/* --Methode-- */
//++
template <class T>
int_4 SphereThetaPhi<T>::NPhi(int_4 kt) const

//    Return number of pixels in phi-direction of the kt-th slice
//-- 
{
  int nbpix;  
  // verification
  if((kt < 0) || (kt >= 2*NTheta_)) return(-1);
  
  // si on se trouve dans l'hemisphere Sud
  if(kt >= NTheta_) {
    kt= 2*NTheta_-1-kt; 
  }
  
  // nombre de pixels
  nbpix= NPhi_(kt);
  return(nbpix);
}


/* --Methode-- */
//++
template <class T>
void SphereThetaPhi<T>::Theta(int_4 kt,double& tetMin,double& tetMax) 

//    Return  theta values which limit the slice kt
//--
{
  bool fgzn= false;
  // verification
  if( (kt< 0) || (kt>= 2*NTheta_) ) {
    tetMin= -99999.; 
    tetMax= -99999.; 
    return;
  }

  // si on se trouve dans l'hemisphere Sud
  if( kt >= NTheta_ ) {
    fgzn= true;
    kt= 2*NTheta_-1-kt; 
  }
  
  // valeurs limites de theta dans l'hemisphere Nord
  tetMin= Theta_(kt);
  tetMax= Theta_(kt+1);
  // valeurs limites de theta dans l'hemisphere Sud
  if (fgzn) {
    tetMin= Pi - Theta_(kt+1);
    tetMax= Pi - Theta_(kt);
  }  
}

/* --Methode-- */
//++
template <class T>
void SphereThetaPhi<T>::Phi(int_4 kt,int_4 jp,double& phiMin,double& phiMax) 

//   Return values of phi which limit the jp-th pixel of the kt-th slice
//--
{
  // verification
  if((kt < 0) || (kt >= 2*NTheta_)) {
    phiMin= -99999.; 
    phiMax= -99999.; 
    return;
  }
  
  // si on se trouve dans l'hemisphere Sud
  if(kt >= NTheta_) kt= 2*NTheta_-1-kt;
  
  // verifie si la tranche kt contient au moins jp pixels
  if( (jp< 0) || (jp >= NPhi_(kt)) ) {
    phiMin= -88888.; 
    phiMax= -88888.; 
    return;
  }
  
  double dphi= DeuxPi/(double)NPhi_(kt);
  phiMin= dphi*(double)(jp);
  phiMax= dphi*(double)(jp+1);
  return;
}

/* --Methode-- */
//++
template <class T>
int_4 SphereThetaPhi<T>::Index(int_4 kt,int_4 jp) const

//    Return pixel index  with sequence index jp in the slice kt 
//--
{
  int k;
  bool fgzn= false;
  
  // si on se trouve dans l'hemisphere Sud
  if(kt >= NTheta_) {
    fgzn= true; 
    kt= 2*NTheta_-1-kt;
  }
  
  // si la tranche kt contient au moins i pixels
  if( (jp>=0) && (jp<NPhi_(kt)) ) 
    {
      // dans l'hemisphere Sud
      if (fgzn) k= NPix_-TNphi_(kt+1)+jp;
      // dans l'hemisphere Nord
      else k= TNphi_(kt)+jp;
    }
  else
    {
      k= 9999;
      printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
    }
  return(k);
}

/* --Methode-- */
//++
template <class T>
void SphereThetaPhi<T>::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp) 

//    Return indices kt (theta) and jp (phi) of  pixel with index k
//--
{
  bool fgzn= false;  
  // si on se trouve dans l'hemisphere Sud
  if(k >= NPix_/2) 
    { 
      fgzn= true; 
      k= NPix_-1-k; 
    }
  
  // on recupere l'indice kt de la tranche qui contient le pixel k
  int i;
  for(i = 0; i < NTheta_; i++) 
    if(k < TNphi_(i+1)) break;
  
  // indice  kt de tranche
  if (fgzn) kt= 2*NTheta_-1-i;
  else kt= i;
  
  // indice jp de pixel
  if (fgzn) jp= TNphi_(i+1)-k-1;
  else jp= k-TNphi_(i);  
}
//++
template <class T>
void  SphereThetaPhi<T>::Pixelize(int_4 m) 

//    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).
//--
{
  int ntotpix,j;
  
  // Decodage et controle des arguments d'appel :
  // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
  if (m < 2) m = 2;
  if (m > 16384) m = 16384;
  
  // On memorise les arguments d'appel
  NTheta_ = m;  
  
  // On commence par decouper l'hemisphere z>0.
  // Creation des vecteurs contenant :
  // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
  //  Theta_= new double[m+1];
  Theta_.ReSize(m+1);
  
  // Le nombre de pixels en phi de chacune des bandes en theta
  //  NPhi_ = new int_4[m];
  // une taille de m suffit, mais je mets m+1 pour que les 3 tableaux aient 
  // la meme taille pour une manipulation plus faciles par la librairie 
  // cfitsio -- GLM (13-04-00)
  NPhi_.ReSize(m+1);
  
  // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
  // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
  //  TNphi_= new int_4[m+1];
    TNphi_.ReSize(m+1);
  
  // Calcul du nombre total de pixels dans chaque bande pour optimiser
  // le rapport largeur/hauteur des pixels 
  
  //calotte polaire
  TNphi_(0)= 0;
  NPhi_(0) = 1;
  
  //bandes jusqu'a l'equateur
  for(j = 1; j < m; j++) 
    {
      TNphi_(j)= TNphi_(j-1)+NPhi_(j-1);
      NPhi_(j) = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
    }
  
  // Nombre total de pixels sur l'hemisphere
  ntotpix  = TNphi_(m-1)+NPhi_(m-1);
  TNphi_(m)= ntotpix;
  // et sur la sphere entiere
  NPix_= 2*ntotpix;
  
  // Creation et initialisation du vecteur des contenus des pixels 
  pixels_.ReSize(NPix_);
  pixels_.Reset();

  // Determination des limites des bandes en theta :
  // omeg est l'angle solide couvert par chaque pixel,
  // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg
  // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide 
  //de la
  // calotte allant du pole a la limite haute de la bande kt vaut
  // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...
  
  double omeg2pi= 1./(double)ntotpix;
  Omega_ = omeg2pi*DeuxPi;
  
  for(j=0; j <= m; j++) 
    {
      Theta_(j)= acos(1.-(double)TNphi_(j)*omeg2pi);
    }
}

//++
template <class T>
void SphereThetaPhi<T>::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" 
//    a vector containing the phi's of the pixels of the slice
//    a vector containing the corresponding values of pixels 
//--

{

  if(index < 0 || index >= NbThetaSlices()) 
    {
       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range"); 
    }

  int iring= Index(index,0);
  int lring  = NPhi(index);

  phi.ReSize(lring);
  value.ReSize(lring);
  double Te= 0.;
  double Fi= 0.;
  for(int kk = 0; kk < lring; kk++) 
    {
      PixThetaPhi(kk+iring,Te,Fi);
      phi(kk)= Fi;
      value(kk)= PixVal(kk+iring);
    }
  theta= Te;
}

//++
template <class T>
void SphereThetaPhi<T>::GetThetaSlice(int_4 index,r_8& theta, r_8& phi0,TVector<int_4>& pixelIndices, 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 the indices of the pixels of the slice
//    (equally distributed in phi)
//    a vector containing the corresponding values of pixels 
//--

{

  if(index < 0 || index >= NbThetaSlices()) 
    {
       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range"); 
 
    }

  int iring= Index(index,0);
  int lring  = NPhi(index);

  pixelIndices.ReSize(lring);
  value.ReSize(lring);
  double Te= 0.;
  double Fi= 0.;
  for(int kk = 0; kk < lring; kk++) 
    {
      pixelIndices(kk)=kk+iring ;
      value(kk)= PixVal(kk+iring);
    }
      PixThetaPhi(iring,theta,phi0);
}




template <class T>
void SphereThetaPhi<T>::print(ostream& os) const
{
  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
  //
  os << " NTheta_= " << NTheta_ << endl;
  os << " NPix_    = " << NPix_   << endl;
  os << " Omega_  =  " << Omega_   << endl;

  os << " contenu de NPhi_ : ";
  int i;
  for(i=0; i < NTheta_; i++)
    {
      if(i%5 == 0) os << endl;
      os << NPhi_(i) <<", ";
    }
  os << endl;

  os << " contenu de Theta_ : ";
  for(i=0; i < NTheta_+1; i++)
    {
      if(i%5 == 0) os << endl;
      os << Theta_(i) <<", ";
    }
  os << endl;

  os << " contenu de TNphi_ : ";
  for(i=0; i < NTheta_+1; i++)
    {
      if(i%5 == 0) os << endl;
      os << TNphi_(i) <<", ";
    }
  os << endl;

  os << " contenu de pixels : ";
  for(i=0; i < NPix_; i++)
    {
      if(i%5 == 0) os << endl;
      os <<  pixels_(i) <<", ";
    }
  os << endl;
}


#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template SphereThetaPhi<r_8>
#pragma define_template SphereThetaPhi<r_4>
#pragma define_template SphereThetaPhi< complex<r_4> >
#pragma define_template SphereThetaPhi< complex<r_8> >
#endif
#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class SphereThetaPhi<r_8>;
template class SphereThetaPhi<r_4>;
template class SphereThetaPhi< complex<r_4> >;
template class SphereThetaPhi< complex<r_8> >;
#endif