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

#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();
}


/* --Methode-- */

//++
template <class T>
SphereThetaPhi<T>::SphereThetaPhi(int 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_  = new int[NTheta_];
  Theta_ = new double[NTheta_+1];
  TNphi_ = new int[NTheta_+1];
  for(int k = 0; k < NTheta_; k++) 
    {
      NPhi_[k] = s.NPhi_[k];
      Theta_[k]= s.Theta_[k];
      TNphi_[k]= s.TNphi_[k];
    }
  Theta_[NTheta_]= s.Theta_[NTheta_];
  TNphi_[NTheta_]= s.TNphi_[NTheta_];
  Omega_ = s.Omega_;
  pixels_= s.pixels_;
}

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

//--
{
  Clear();
}

//++
// Titre        Public Méthods
//--
template <class T>
void SphereThetaPhi<T>::InitNul()
//
//    initialise à zéro les variables de classe pointeurs
{
  NTheta_= 0;
  NPix_  = 0;
  Theta_ = NULL;
  NPhi_  = NULL;
  TNphi_ = NULL;
  pixels_.Reset();
}

/* --Methode-- */
template <class T>
void SphereThetaPhi<T>::Clear()

{
  if(Theta_)  delete[] Theta_;
  if(NPhi_ )  delete[] NPhi_;
  if(TNphi_)  delete[] TNphi_;
  InitNul();
}

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

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

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

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

//    Return value of pixel with index k
//--
{
  if((k < 0) || (k >= NPix_)) 
    {
      //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
      cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl;
      THROW(rangeCheckErr);
    }
  return pixels_(k);
}

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

//    Return value of pixel with index k
//--
{
  if((k < 0) || (k >= NPix_)) 
    {
      cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
      //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
      throw "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 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)(phi/dphi);
  else k= TNphi_[i-1]+(int)(phi/dphi); 
  return(k);
}

/* --Methode-- */
//++
template <class T>
void SphereThetaPhi<T>::PixThetaPhi(int 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 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 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>
int SphereThetaPhi<T>::NbThetaSlices() const

//    Return number of theta-slices on the sphere
//--
{
  int nbslices;
  nbslices= 2*NTheta_;
  return(nbslices);
}

/* --Methode-- */
//++
template <class T>
int SphereThetaPhi<T>::NPhi(int 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 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 kt,int 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 SphereThetaPhi<T>::Index(int kt,int 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 k,int& kt,int& 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 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];
  
  // Le nombre de pixels en phi de chacune des bandes en theta
  NPhi_ = new int[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[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)(.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 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 
//--

{
  cout << "entree GetThetaSlice, couche no " << index << endl;

  if(index < 0 || index > NbThetaSlices()) 
    {
      // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 
      cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
      THROW(rangeCheckErr);
    }

  int iring= Index(index,0);
  int bid  = this->NPhi(index);
  int lring  = bid;
  cout << " iring= " << iring << " lring= " << lring << endl;

  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>::setmNPhi(int* array, int m)
  //remplit le tableau contenant le nombre de pixels en phi de chacune des bandes en theta 
  //
{
  NPhi_= new int[m];
  for(int k = 0; k < m; k++) NPhi_[k]= array[k];
}

template <class T>
void SphereThetaPhi<T>::setmTNphi(int* array, int m)
  //remplit  le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes 
  //
{
  TNphi_= new int[m];
  for(int k = 0; k < m; k++) TNphi_[k]= array[k];
}

template <class T>
void SphereThetaPhi<T>::setmTheta(double* array, int m)
  //remplit  le tableau contenant les valeurs limites de theta
  //
{
  Theta_= new double[m];
  for(int k = 0; k < m; k++) Theta_[k]= array[k];
}

template <class T>
void SphereThetaPhi<T>::setDataBlock(T* data, int m)
  // remplit  le vecteur des contenus des pixels 
{
  pixels_.FillFrom(m,data);
} 

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>;
  ownobj= true;
  Read(filename);
}

template <class T>
FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
{
  dobj= new SphereThetaPhi<T>(obj, 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>;
    }

// Let's Read the SphereCoordSys object  -- ATTENTIOn - $CHECK$
  SphereCoordSys* cs = dynamic_cast<SphereCoordSys*>(is.ReadObject());
  dobj->SetCoordSys(cs);

  // 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 mNTheta;
  is.GetI4(mNTheta);  
  dobj->setSizeIndex(mNTheta);

  int mNPix;
  is.GetI4(mNPix);
  dobj->setNbPixels(mNPix);

  double mOmeg;
  is.GetR8(mOmeg);
  dobj->setPixSolAngle(mOmeg);

  int* mNphi= new int[mNTheta];
  is.GetI4s(mNphi, mNTheta);
  dobj->setmNPhi(mNphi, mNTheta);
  delete [] mNphi;

  int* mTNphi= new int[mNTheta+1];
  is.GetI4s(mTNphi, mNTheta+1);
  dobj->setmTNphi(mTNphi, mNTheta+1);
  delete [] mTNphi;

  double* mTheta= new double[mNTheta+1];
  is.GetR8s(mTheta, mNTheta+1);
  dobj->setmTheta(mTheta, mNTheta+1);
  delete [] mTheta;

  T* mPixels= new T[mNPix];
  PIOSReadArray(is, mPixels, mNPix);
  dobj->setDataBlock(mPixels, mNPix);
  delete [] mPixels;
}

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

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

// Let's write the SphereCoordSys object
  dobj->GetCoordSys()->Write(os);
  
  char strg[256];
  int mNTheta= dobj->SizeIndex();
  int 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(0));
  os.PutI4s(dobj->getmNPhi() , mNTheta);
  os.PutI4s(dobj->getmTNphi(), mNTheta+1);
  os.PutR8s(dobj->getmTheta(), mNTheta+1);
  //os.Put((dobj->getDataBlock())->Data(), mNPix);
  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix);
}

#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