source: Sophya/trunk/SophyaLib/Samba/spherethetaphi.cc@ 738

#include "spherethetaphi.h"
#include "nbmath.h"
#include <complex>
#include "piocmplx.h"
#include <iostream.h>


//***************************************************************
//++ 
// Class        SphereThetaPhi
// 
//
// include      spherethetaphi.h nbmath.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()

//--
{
  InitNul();
  pixels_.Reset();
}


/* --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)
  : pixels_(s.pixels_ , share)
{
  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
  NTheta_= s.NTheta_;
  NPix_  = s.NPix_;
  NPhi_.ReSize(NTheta_);
  Theta_.ReSize(NTheta_+1);
  TNphi_.ReSize(NTheta_+1);

  NPhi_  = s.NPhi_;
  TNphi_ = s.TNphi_;
  Theta_ = s.Theta_;
  Omega_ = s.Omega_;
}

//++
// 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);
}

/* --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];
  NPhi_.ReSize(m);
  
  // 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,double& theta, TVector<double>& 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,double& theta, double& phi0,TVector<int>& 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;
}

///////////////////////////////////////////////////////////
// --------------------------------------------------------
//   Les objets delegues pour la gestion de persistance 
// --------------------------------------------------------
//////////////////////////////////////////////////////////

template <class T>
FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi()
{
  dobj= new SphereThetaPhi<T>;
  ownobj= true;
}

template <class T>
FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename)
{
  dobj= new SphereThetaPhi<T>;
  dobj->pixels_.SetTemp(true);
  ownobj= true;
  Read(filename);
}

template <class T>
FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
{
  dobj= new SphereThetaPhi<T>(obj, true);
  dobj->pixels_.SetTemp(true);
  ownobj= true;
}

template <class T>
FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj)
{
  dobj= obj;
  ownobj= false;
}

template <class T>
FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi()
{
  if (ownobj && dobj) delete dobj;
}

template <class T>
AnyDataObj* FIO_SphereThetaPhi<T>::DataObj()
{
  return(dobj);
}

template <class T>
void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is)
{

  if(dobj == NULL) 
    {
      dobj= new SphereThetaPhi<T>;
      dobj->pixels_.SetTemp(true);
      ownobj= true;      
    }

  // On lit les 3 premiers uint_8
  uint_8 itab[3];
  is.Get(itab, 3);
// Let's Read the SphereCoordSys object  -- ATTENTIOn - $CHECK$
  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
  fio_scs.Read(is);

  // Pour savoir s'il y avait un DVList Info associe
  char strg[256];
  is.GetLine(strg, 255);
  bool hadinfo= false;
  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
  if(hadinfo) 
    {    // Lecture eventuelle du DVList Info
      is >> dobj->Info(); 
    }

  int_4 mNTheta;
  is.GetI4(mNTheta);  
  int_4 mNPix;
  is.GetI4(mNPix);
  double mOmeg;
  is.GetR8(mOmeg);
  dobj->setParameters(mNPix, mOmeg, mNTheta); 
  FIO_NDataBlock<int_4> fio_nd_nphi(&dobj->NPhi_);
  fio_nd_nphi.Read(is);
  FIO_NDataBlock<int_4> fio_nd_Tnphi(&dobj->TNphi_);
  fio_nd_Tnphi.Read(is);
  FIO_NDataBlock<r_8> fio_nd_Theta(&dobj->Theta_);
  fio_nd_Theta.Read(is);

FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
fio_nd.Read(is);
}

template <class T>
void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const
{

  if(dobj == NULL) 
    {
      cout << " WriteSelf:: dobj= null " << endl;
      return;
    }

//  On ecrit 3 uint_8 
//  0 : Numero de version,  1 : Size index,  2  reserve a l
uint_8 itab[3];
itab[0] = 1;
itab[1] = dobj->SizeIndex();
itab[2] = 0;
os.Put(itab, 3);
// Let's write the SphereCoordSys object
  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
  fio_scs.Write(os);
  
  char strg[256];
  int_4 mNTheta= dobj->SizeIndex();
  int_4 mNPix  = dobj->NbPixels();
 
  if(dobj->ptrInfo()) 
    {
      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo",mNTheta,mNPix);
      os.PutLine(strg);
      os << dobj->Info();
    }
  else 
    { 
      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d ",mNTheta,mNPix);
      os.PutLine(strg);  
    }

  os.PutI4(mNTheta);  
  os.PutI4(mNPix);
  os.PutR8(dobj->PixSolAngle());
  FIO_NDataBlock<int_4> fio_nd_nphi(&dobj->NPhi_);
  fio_nd_nphi.Write(os);
  FIO_NDataBlock<int_4> fio_nd_Tnphi(&dobj->TNphi_);
  fio_nd_Tnphi.Write(os);
  FIO_NDataBlock<r_8> fio_nd_Theta(&dobj->Theta_);
  fio_nd_Theta.Write(os);
FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
fio_nd.Write(os);
}

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