Ignore:
Timestamp:
Oct 15, 1999, 5:43:30 PM (26 years ago)
Author:
ansari
Message:

versions templatees, NdataBlocks etc. 15-OCT-99-GLM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/Samba/spherethetaphi.cc

    r276 r470  
    1 //
    2 //
    31#include "spherethetaphi.h"
    42#include "nbmath.h"
     3#include <complex>
     4#include "piocmplx.h"
    55#include <iostream.h>
     6
     7#ifdef __CXX_PRAGMA_TEMPLATES__
     8#pragma define_template SphereThetaPhi<double>
     9#pragma define_template SphereThetaPhi<float>
     10#pragma define_template SphereThetaPhi< complex<float> >
     11#pragma define_template SphereThetaPhi< complex<double> >
     12#pragma define_template FIO_SphereThetaPhi<double>
     13#pragma define_template FIO_SphereThetaPhi<float>
     14#pragma define_template FIO_SphereThetaPhi< complex<float> >
     15#pragma define_template FIO_SphereThetaPhi< complex<double> >
     16#endif
     17#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
     18template class SphereThetaPhi<double>;
     19template class SphereThetaPhi<float>;
     20template class SphereThetaPhi< complex<float> >;
     21template class SphereThetaPhi< complex<double> >;
     22template class FIO_SphereThetaPhi<double>;
     23template class FIO_SphereThetaPhi<float>;
     24template class FIO_SphereThetaPhi< complex<float> >;
     25template class FIO_SphereThetaPhi< complex<double> >;
     26#endif
     27
    628//***************************************************************
    729//++
    … …  
    3961//++
    4062
    41 
    42 SphereThetaPhi::SphereThetaPhi()
     63template <class T>
     64SphereThetaPhi<T>::SphereThetaPhi()
    4365
    4466//--
    … …  
    4971/* --Methode-- */
    5072//++
    51 SphereThetaPhi::SphereThetaPhi(char* flnm)
     73template <class T>
     74SphereThetaPhi<T>::SphereThetaPhi(char* flnm)
    5275
    5376//    Constructeur : charge une image à partir d'un fichier
    … …  
    5578{
    5679  InitNul();
    57   PInPersist s(flnm);
    58   Read(s);
    59 }
    60 
    61 /* --Methode-- */
    62 
    63 //++
    64 SphereThetaPhi::SphereThetaPhi(int_4 m, int_4 pet)
     80  cout << " ShereThetaPhi:: constructeur lecture sur fichier A FAIRE " << endl;
     81  //PInPersist s(flnm);
     82  //Read(s);
     83}
     84
     85/* --Methode-- */
     86
     87//++
     88template <class T>
     89SphereThetaPhi<T>::SphereThetaPhi(int_4 m)
    6590
    6691//    Constructeur : m est le nombre de bandes en theta sur un hémisphère
    … …  
    7196{
    7297  InitNul();
    73   Pixelize(m,pet);
    74 }
    75 
     98  Pixelize(m);
     99}
     100
     101template <class T>
     102SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)
     103{
     104  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
     105  NTheta_= s.NTheta_;
     106  NPix_  = s.NPix_;
     107  NPhi_  = new int_4[NTheta_];
     108  Theta_ = new r_4[NTheta_+1];
     109  TNphi_ = new int_4[NTheta_+1];
     110  for(int k = 0; k < NTheta_; k++)
     111    {
     112      NPhi_[k] = s.NPhi_[k];
     113      Theta_[k]= s.Theta_[k];
     114      TNphi_[k]= s.TNphi_[k];
     115    }
     116  Theta_[NTheta_]= s.Theta_[NTheta_];
     117  TNphi_[NTheta_]= s.TNphi_[NTheta_];
     118  Omega_ = s.Omega_;
     119  pixels_= s.pixels_;
     120}
    76121
    77122//++
    … …  
    79124//--
    80125//++
    81 SphereThetaPhi::~SphereThetaPhi()
     126template <class T>
     127SphereThetaPhi<T>::~SphereThetaPhi()
    82128
    83129//--
    … …  
    89135// Titre        Méthodes
    90136//--
    91 void SphereThetaPhi::InitNul()
     137template <class T>
     138void SphereThetaPhi<T>::InitNul()
    92139//
    93140//    initialise à zéro les variables de classe pointeurs
    94141{
    95   mTheta = NULL;
    96   mNphi = NULL;
    97   mTNphi = NULL;
    98   // mPix = NULL;
    99   _pixel=NULL;
    100   mNTheta = mNPet = 0;
    101   mNPix = 0;
    102 }
    103 
    104 /* --Methode-- */
    105 void SphereThetaPhi::Clear()
    106 
    107 {
    108   if (mTheta)  delete[] mTheta;
    109   if (mNphi)  delete[] mNphi;
    110   if (mTNphi) delete[] mTNphi;
    111   //if (mPix)   delete[] mPix;
    112   if (_pixel) _pixel->Detach();
    113   mTheta = NULL;
    114   mNphi = NULL;
    115   mTNphi = NULL;
    116   //mPix = NULL;
    117   _pixel=NULL;
    118   mNTheta = mNPet = 0;
    119   mNPix = 0;
    120 }
    121 
    122 /* --Methode-- */
    123 //++
    124 void SphereThetaPhi::WriteSelf(POutPersist& s) const
    125 
    126 //    créer un fichier image
    127 //--
    128 {
    129   char strg[256];
    130   if (mInfo_) {sprintf(strg, "SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo", (int_4)mNTheta, (int_4)mNPix);
    131   }
    132   else { sprintf(strg, "SphereThetaPhi: NSlices=%6d  NPix=%9d ",
    133                (int_4)mNTheta, (int_4)mNPix);
    134   }
    135   s.PutLine(strg);
    136   if (mInfo_)  s << (*mInfo_);
    137   s.PutI4(mNTheta);
    138   s.PutI4(mNPet);
    139   s.PutI4(mNPix);
    140   s.PutR8(mOmeg);
    141   s.PutR4s(mTheta, mNTheta+1);
    142   s.PutI4s(mNphi, mNTheta);
    143   s.PutI4s(mTNphi, mNTheta+1);
    144   s.PutR8s(_pixel->Data(), mNPix);
    145   //s.PutR8s(mPix, mNPix);
    146   return;
    147 }
    148 
    149 /* --Methode-- */
    150 //++
    151 void SphereThetaPhi::ReadSelf(PInPersist& s)
    152 
    153 //    relit un fichier d'image
     142  NTheta_= 0;
     143  NPix_  = 0;
     144  Theta_ = NULL;
     145  NPhi_  = NULL;
     146  TNphi_ = NULL;
     147  pixels_.Reset();
     148}
     149
     150/* --Methode-- */
     151template <class T>
     152void SphereThetaPhi<T>::Clear()
     153
     154{
     155  if(Theta_)  delete[] Theta_;
     156  if(NPhi_ )  delete[] NPhi_;
     157  if(TNphi_)  delete[] TNphi_;
     158  InitNul();
     159}
     160
     161//++
     162template <class T>
     163void SphereThetaPhi<T>::Resize(int_4 m)
     164//   re-pixelize the sphere
    154165//--
    155166{
    156167  Clear();
    157   char strg[256];
    158   s.GetLine(strg, 255);
    159 // Pour savoir s'il y avait un DVList Info associe
    160   bool hadinfo = false;
    161   if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
    162   if (hadinfo) {    // Lecture eventuelle du DVList Info
    163     if (mInfo_ == NULL)  mInfo_ = new DVList;
    164     s >> (*mInfo_);
    165   }
    166   s.GetI4(mNTheta);
    167   s.GetI4(mNPet);
    168   s.GetI4(mNPix);
    169   s.GetR8(mOmeg);
    170   mTheta = new r_4[mNTheta+1];
    171   mNphi = new int_4[mNTheta];
    172   mTNphi = new int_4[mNTheta+1];
    173   //mPix = new r_8[mNPix+1];
    174   _pixel=new PDataArray<r_8>(mNPix,true);
    175   s.GetR4s(mTheta, mNTheta+1);
    176   s.GetI4s(mNphi, mNTheta);
    177   s.GetI4s(mTNphi, mNTheta+1);
    178   s.GetR8s(_pixel->Data(), mNPix);
    179   //s.GetR8s(mPix, mNPix);
    180   return;
    181 }
    182 
    183 
    184 /* --Methode-- */
    185 //++
    186 int_4 SphereThetaPhi::NbPixels() const
     168  Pixelize(m);
     169}
     170
     171/* --Methode-- */
     172//++
     173template <class T>
     174int_4 SphereThetaPhi<T>::NbPixels() const
    187175
    188176//    Retourne le nombre de pixels du découpage
    189177//--
    190178{
    191   return(mNPix);
    192 }
    193 
    194 
    195 static r_8 dummy_pixel = 0;
    196 
    197 /* --Methode-- */
    198 //++
    199 r_8& SphereThetaPhi::PixVal(int_4 k)
     179  return(NPix_);
     180}
     181
     182/* --Methode-- */
     183//++
     184template <class T>
     185T& SphereThetaPhi<T>::PixVal(int_4 k)
    200186
    201187//    Retourne la valeur du contenu du pixel d'indice k
    202188//--
    203189{
    204   if ( (k<0) || (k >= mNPix) ) {
    205     //  THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    206     cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
    207   THROW(rangeCheckErr);
    208 
    209     //throw "SphereThetaPhi::PIxVal Pixel index out of range";
    210   }
    211   return (*(_pixel->Data()+k));
    212   //return(mPix[k]);
    213 }
    214 //++
    215 r_8 const& SphereThetaPhi::PixVal(int_4 k) const
     190  if((k < 0) || (k >= NPix_))
     191    {
     192      //  THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
     193      cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
     194      THROW(rangeCheckErr);
     195    }
     196  return pixels_(k);
     197}
     198
     199//++
     200template <class T>
     201T const& SphereThetaPhi<T>::PixVal(int_4 k) const
    216202
    217203//    Retourne la valeur du contenu du pixel d'indice k
    218204//--
    219205{
    220   if ( (k<0) || (k >= mNPix) ) {
    221     cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
    222   //    THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    223 
    224     throw "SphereThetaPhi::PIxVal Pixel index out of range";
    225   }
    226   return *(_pixel->Data()+k);
    227   //return(mPix[k]);
    228 }
    229 
    230 
    231 /* --Methode-- */
    232 //++
    233 int_4 SphereThetaPhi::PixIndexSph(r_4 theta, r_4 phi) const
     206  if((k < 0) || (k >= NPix_))
     207    {
     208      cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
     209      //    THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
     210      throw "SphereThetaPhi::PIxVal Pixel index out of range";
     211    }
     212  return *(pixels_.Data()+k);
     213}
     214
     215/* --Methode-- */
     216//++
     217template <class T>
     218int_4 SphereThetaPhi<T>::PixIndexSph(r_4 theta, r_4 phi) const
    234219
    235220//    Retourne l'indice du pixel vers lequel pointe une direction définie par
    … …  
    246231 
    247232  // La bande d'indice kt est limitée par les valeurs de theta contenues dans
    248   // mTheta[kt] et mTheta[kt+1]
    249   for( i=1; i< mNTheta; i++ )
    250     if( theta < mTheta[i] ) break;
    251  
    252   dphi= (r_4)DeuxPi/(r_4)mNphi[i-1];
    253  
    254   if (fgzn)  k= mNPix-mTNphi[i]+(int_4)(phi/dphi);
    255   else k= mTNphi[i-1]+(int_4)(phi/dphi);
    256  
     233  // Theta_[kt] et Theta_[kt+1]
     234  for( i=1; i< NTheta_; i++ )
     235    if( theta < Theta_[i] ) break;
     236 
     237  dphi= (r_4)DeuxPi/(r_4)NPhi_[i-1];
     238 
     239  if (fgzn)  k= NPix_-TNphi_[i]+(int_4)(phi/dphi);
     240  else k= TNphi_[i-1]+(int_4)(phi/dphi);
    257241  return(k);
    258242}
    259243
    260 
    261 /* --Methode-- */
    262 //++
    263 void SphereThetaPhi::PixThetaPhi(int_4 k, r_4& theta, r_4& phi) const
     244/* --Methode-- */
     245//++
     246template <class T>
     247void SphereThetaPhi<T>::PixThetaPhi(int_4 k, r_4& theta, r_4& phi) const
    264248
    265249//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
    … …  
    269253  bool fgzn = false;
    270254 
    271   if( (k < 0) || (k >= mNPix) ) {theta = -99999.; phi = -99999.;  return; }
    272   if( k >= mNPix/2)  {fgzn = true;  k = mNPix-1-k; }
     255  if( (k < 0) || (k >= NPix_) ) {theta = -99999.; phi = -99999.;  return; }
     256  if( k >= NPix_/2)  {fgzn = true;  k = NPix_-1-k; }
    273257
    274258  // recupère l'indice i de la tranche qui contient le pixel k
    275   for( i=0; i< mNTheta; i++ )
    276     if( k < mTNphi[i+1] ) break;
     259  for( i=0; i< NTheta_; i++ )
     260    if( k < TNphi_[i+1] ) break;
    277261
    278262  // angle theta
    279   theta= 0.5*(mTheta[i]+mTheta[i+1]);
     263  theta= 0.5*(Theta_[i]+Theta_[i+1]);
    280264  if (fgzn) theta= Pi-theta;
    281265 
    282266  // angle phi
    283   k -= mTNphi[i];
    284   phi= (r_4)DeuxPi/(r_4)mNphi[i]*(r_4)(k+.5);
     267  k -= TNphi_[i];
     268  phi= (r_4)DeuxPi/(r_4)NPhi_[i]*(r_4)(k+.5);
    285269  if (fgzn) phi= DeuxPi-phi;
    286 
    287   return;
    288 }
    289 
    290 
    291 //++
    292 r_8         SphereThetaPhi::PixSolAngle(int_4 dummy) const
     270}
     271
     272//++
     273template <class T>
     274r_8  SphereThetaPhi<T>::PixSolAngle(int_4 dummy) const
     275
    293276//    Pixel Solid angle  (steradians)
    294277//    All the pixels have the same solid angle. The dummy argument is
    … …  
    297280//--
    298281{
    299   return mOmeg;
    300 }
    301 
    302 
    303 /* --Methode-- */
    304 //++
    305 void SphereThetaPhi::Limits(int_4 k, r_4& tetMin, r_4& tetMax,
    306                            r_4& phiMin, r_4& phiMax)
     282  return Omega_;
     283}
     284
     285/* --Methode-- */
     286//++
     287template <class T>
     288void SphereThetaPhi<T>::Limits(int_4 k, r_4& tetMin, r_4& tetMax, r_4& phiMin, r_4& phiMax)
    307289
    308290//    Retourne les valeurs de theta et phi limitant le pixel d'indice k
    … …  
    313295  bool fgzn= false;
    314296 
    315   if( (k< 0) || (k >= mNPix) ) {
     297  if( (k< 0) || (k >= NPix_) ) {
    316298    tetMin= -99999.;
    317299    phiMin= -99999.; 
    … …  
    322304 
    323305  // si on se trouve dans l'hémisphère Sud
    324   if( k >= mNPix/2 ) {
     306  if( k >= NPix_/2 ) {
    325307    fgzn= true;
    326     k= mNPix-1-k;
     308    k= NPix_-1-k;
    327309  }
    328310 
    329311  // on recupere l'indice i de la tranche qui contient le pixel k
    330312  int i;
    331   for( i=0; i< mNTheta; i++ )
    332     if( k< mTNphi[i+1] ) break;
     313  for( i=0; i< NTheta_; i++ )
     314    if( k< TNphi_[i+1] ) break;
    333315 
    334316  // valeurs limites de theta dans l'hemisphere Nord
    335   tetMin= mTheta[i];
    336   tetMax= mTheta[i+1];
     317  tetMin= Theta_[i];
     318  tetMax= Theta_[i+1];
    337319  // valeurs limites de theta dans l'hemisphere Sud
    338320  if (fgzn) {
    339     tetMin= Pi-mTheta[i+1];
    340     tetMax= Pi-mTheta[i];
     321    tetMin= Pi-Theta_[i+1];
     322    tetMax= Pi-Theta_[i];
    341323  }
    342324 
    343325  // pixel correspondant dans l'hemisphere Nord
    344   if (fgzn) k= mTNphi[i+1]-k+mTNphi[i]-1;
     326  if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1;
    345327 
    346328  // indice j de discretisation ( phi= j*dphi )
    347   j= k-mTNphi[i];
    348   dphi= (r_4)DeuxPi/(r_4)mNphi[i];
     329  j= k-TNphi_[i];
     330  dphi= (r_4)DeuxPi/(r_4)NPhi_[i];
    349331 
    350332  // valeurs limites de phi
    … …  
    356338/* --Methode-- */
    357339//++
    358 int_4 SphereThetaPhi::NbThetaSlices() const
     340template <class T>
     341int_4 SphereThetaPhi<T>::NbThetaSlices() const
    359342
    360343//    Retourne le nombre de tranches en theta sur la sphere
    … …  
    362345{
    363346  int_4 nbslices;
    364   nbslices= 2*mNTheta;
     347  nbslices= 2*NTheta_;
    365348  return(nbslices);
    366349}
    … …  
    368351/* --Methode-- */
    369352//++
    370 int_4 SphereThetaPhi::NPhi(int_4 kt) const
     353template <class T>
     354int_4 SphereThetaPhi<T>::NPhi(int_4 kt) const
    371355
    372356//    Retourne le nombre de pixels en phi de la tranche kt
    373357//--
    374358{
    375  
    376   int_4 nbpix;
    377  
     359  int_4 nbpix; 
    378360  // verification
    379   if( (kt< 0) || (kt>= 2*mNTheta) ) return(-1);
     361  if( (kt< 0) || (kt>= 2*NTheta_) ) return(-1);
    380362 
    381363  // si on se trouve dans l'hemisphere Sud
    382   if( kt >= mNTheta ) {
    383     kt= 2*mNTheta-1-kt;
     364  if( kt >= NTheta_ ) {
     365    kt= 2*NTheta_-1-kt;
    384366  }
    385367 
    386368  // nombre de pixels
    387   nbpix= mNphi[kt];
     369  nbpix= NPhi_[kt];
    388370  return(nbpix);
    389371}
    … …  
    392374/* --Methode-- */
    393375//++
    394 void SphereThetaPhi::Theta(int_4 kt,r_4& tetMin,r_4& tetMax)
     376template <class T>
     377void SphereThetaPhi<T>::Theta(int_4 kt,r_4& tetMin,r_4& tetMax)
    395378
    396379//    Retourne les valeurs de theta limitant la tranche kt
    397380//--
    398381{
    399  
    400382  bool fgzn= false;
    401  
    402383  // verification
    403   if( (kt< 0) || (kt>= 2*mNTheta) ) {
     384  if( (kt< 0) || (kt>= 2*NTheta_) ) {
    404385    tetMin= -99999.;
    405386    tetMax= -99999.;
    … …  
    408389
    409390  // si on se trouve dans l'hemisphere Sud
    410   if( kt >= mNTheta ) {
     391  if( kt >= NTheta_ ) {
    411392    fgzn= true;
    412     kt= 2*mNTheta-1-kt;
     393    kt= 2*NTheta_-1-kt;
    413394  }
    414395 
    415396  // valeurs limites de theta dans l'hemisphere Nord
    416   tetMin= mTheta[kt];
    417   tetMax= mTheta[kt+1];
     397  tetMin= Theta_[kt];
     398  tetMax= Theta_[kt+1];
    418399  // valeurs limites de theta dans l'hemisphere Sud
    419400  if (fgzn) {
    420     tetMin= Pi-mTheta[kt+1];
    421     tetMax= Pi-mTheta[kt];
     401    tetMin= Pi-Theta_[kt+1];
     402    tetMax= Pi-Theta_[kt];
    422403  } 
    423404}
    … …  
    425406/* --Methode-- */
    426407//++
    427 void SphereThetaPhi::Phi(int_4 kt,int_4 jp,r_4& phiMin,r_4& phiMax)
     408template <class T>
     409void SphereThetaPhi<T>::Phi(int_4 kt,int_4 jp,r_4& phiMin,r_4& phiMax)
    428410
    429411//    Retourne les valeurs de phi limitant le pixel jp de la tranche kt
    430412//--
    431413{
    432  
    433   r_4 dphi;
    434  
    435414  // verification
    436   if( (kt< 0) || (kt>= 2*mNTheta) ) {
     415  if( (kt< 0) || (kt>= 2*NTheta_) ) {
    437416    phiMin= -99999.;
    438417    phiMax= -99999.;
    … …  
    441420 
    442421  // si on se trouve dans l'hemisphere Sud
    443   if( kt >= mNTheta ) kt= 2*mNTheta-1-kt;
     422  if( kt >= NTheta_ ) kt= 2*NTheta_-1-kt;
    444423 
    445424  // verifie si la tranche kt contient au moins jp pixels
    446   if( (jp< 0) || (jp >= mNphi[kt]) ) {
     425  if( (jp< 0) || (jp >= NPhi_[kt]) ) {
    447426    phiMin= -88888.;
    448427    phiMax= -88888.;
    … …  
    450429  }
    451430 
    452   dphi= (r_4)DeuxPi/(r_4)mNphi[kt];
     431  r_4 dphi= (r_4)DeuxPi/(r_4)NPhi_[kt];
    453432  phiMin= dphi*(r_4)(jp);
    454433  phiMax= dphi*(r_4)(jp+1);
    … …  
    458437/* --Methode-- */
    459438//++
    460 int_4 SphereThetaPhi::Index(int_4 kt,int_4 jp) const
     439template <class T>
     440int_4 SphereThetaPhi<T>::Index(int_4 kt,int_4 jp) const
    461441
    462442//    Retourne l'indice du pixel d'indice jp dans la tranche kt
    463443//--
    464444{
    465  
    466445  int_4 k;
    467446  bool fgzn= false;
    468447 
    469448  // si on se trouve dans l'hemisphere Sud
    470   if( kt >= mNTheta ) {
     449  if( kt >= NTheta_ ) {
    471450    fgzn= true;
    472     kt= 2*mNTheta-1-kt;
     451    kt= 2*NTheta_-1-kt;
    473452  }
    474453 
    475454  // si la tranche kt contient au moins i pixels
    476   if( (jp>=0) && (jp<mNphi[kt]) ) {
    477     // dans l'hemisphere Sud
    478     if (fgzn) k= mNPix-mTNphi[kt+1]+jp;
    479     // dans l'hemisphere Nord
    480     else k= mTNphi[kt]+jp;
    481   }
    482   else{
    483     k= 9999;
    484     printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
    485   }
     455  if( (jp>=0) && (jp<NPhi_[kt]) )
     456    {
     457      // dans l'hemisphere Sud
     458      if (fgzn) k= NPix_-TNphi_[kt+1]+jp;
     459      // dans l'hemisphere Nord
     460      else k= TNphi_[kt]+jp;
     461    }
     462  else
     463    {
     464      k= 9999;
     465      printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
     466    }
    486467  return(k);
    487468}
    … …  
    489470/* --Methode-- */
    490471//++
    491 void SphereThetaPhi::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp)
     472template <class T>
     473void SphereThetaPhi<T>::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp)
    492474
    493475//    Retourne les indices kt et jp du pixel d'indice k
    494476//--
    495477{
    496  
    497   bool fgzn= false;
    498  
     478  bool fgzn= false; 
    499479  // si on se trouve dans l'hemisphere Sud
    500   if( k >= mNPix/2 ) {
     480  if( k >= NPix_/2 ) {
    501481    fgzn= true;
    502     k= mNPix-1-k;
     482    k= NPix_-1-k;
    503483  }
    504484 
    505485  // on recupere l'indice kt de la tranche qui contient le pixel k
    506486  int i;
    507   for( i=0; i< mNTheta; i++ )
    508     if( k< mTNphi[i+1] ) break;
     487  for( i=0; i< NTheta_; i++ )
     488    if( k< TNphi_[i+1] ) break;
    509489 
    510490  // indice  kt de tranche
    511   if (fgzn) kt= 2*mNTheta-1-i;
     491  if (fgzn) kt= 2*NTheta_-1-i;
    512492  else kt= i;
    513493 
    514494  // indice jp de pixel
    515   if (fgzn) jp= mTNphi[i+1]-k-1;
    516   else jp= k-mTNphi[i];
    517  
    518 }
    519 //++
    520 void  SphereThetaPhi::Pixelize( int_4 m, int_4 pet)
     495  if (fgzn) jp= TNphi_[i+1]-k-1;
     496  else jp= k-TNphi_[i];
     497 
     498}
     499//++
     500template <class T>
     501void  SphereThetaPhi<T>::Pixelize( int_4 m)
    521502
    522503//    effectue le découpage en pixels (m et pet ont la même signification
    … …  
    531512{
    532513  int_4 ntotpix,i,j;
    533   //r_8 x, omeg, omeg2pi, teti, tetm, tet, nphi, costeti;
    534   r_8 omeg2pi;
     514 
    535515  // Decodage et controle des arguments d'appel :
    536516  // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
    … …  
    538518  if (m > 16384) m = 16384;
    539519 
    540   // Au moins 4 et au plus 256 pixels dans la premiere bande decoupee en phi
    541   if (pet < 4) pet = 4;
    542   if (pet > 256) pet = 256;
    543  
    544520  // On memorise les arguments d'appel
    545   mNTheta = m;  mNPet = pet;
     521  NTheta_ = m; 
    546522 
    547523  // On commence par decouper l'hemisphere z>0.
    548524  // Creation des vecteurs contenant :
    549525  // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
    550   mTheta = new r_4[m+1];
     526  Theta_= new r_4[m+1];
    551527 
    552528  // Le nombre de pixels en phi de chacune des bandes en theta
    553   mNphi = new int_4[m];
     529  NPhi_ = new int_4[m];
    554530 
    555531  // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
    556532  // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
    557   mTNphi = new int_4[m+1];
     533  TNphi_= new int_4[m+1];
    558534 
    559535  // Calcul du nombre total de pixels dans chaque bande pour optimiser
    … …  
    561537 
    562538  //calotte polaire
    563   mTNphi[0]=0;
    564   mNphi[0] = 1;
     539  TNphi_[0]= 0;
     540  NPhi_[0] = 1;
    565541 
    566542  //bandes jusqu'a l'equateur
    567   for(j=1; j < m; j++) {
    568     mTNphi[j] = mTNphi[j-1]+mNphi[j-1];
    569     mNphi[j]  = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
    570   }
     543  for(j = 1; j < m; j++)
     544    {
     545      TNphi_[j]= TNphi_[j-1]+NPhi_[j-1];
     546      NPhi_[j] = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
     547    }
    571548 
    572549  // Nombre total de pixels sur l'hemisphere
    573   ntotpix = mTNphi[m-1]+mNphi[m-1];
    574   mTNphi[m] = ntotpix;
     550  ntotpix  = TNphi_[m-1]+NPhi_[m-1];
     551  TNphi_[m]= ntotpix;
    575552  // et sur la sphere entiere
    576   mNPix=2*ntotpix;
     553  NPix_= 2*ntotpix;
    577554 
    578555  // Creation et initialisation du vecteur des contenus des pixels
    579   //mPix = new r_8[mNPix];
    580   //for(i=0; i<mNPix; i++)  mPix[i] = 0.;
    581   _pixel=new PDataArray<r_8>(mNPix,true);
    582   for(i=0; i<mNPix; i++)  *(_pixel->Data()+i)=0.;
     556  pixels_.ReSize(NPix_);
     557  pixels_.Reset();
     558
    583559  // Determination des limites des bandes en theta :
    584560  // omeg est l'angle solide couvert par chaque pixel,
    585   // une bande donnee kt couvre un angle solide mNphi[kt]*omeg
    586   // egal a 2* Pi*(cos mTheta[kt]-cos mTheta[kt+1]). De meme, l'angle solide
     561  // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg
     562  // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide
    587563  //de la
    588564  // calotte allant du pole a la limite haute de la bande kt vaut
    589   // 2* Pi*(1.-cos mTheta[kt+1])= mTNphi[kt]*omeg...
    590  
    591   omeg2pi = 1./(r_8)ntotpix;
    592   mOmeg = omeg2pi*DeuxPi;
    593  
    594   for(j=0; j <= m; j++) {
    595     mTheta[j] = acos(1.-(double)mTNphi[j]*omeg2pi);
    596   }
    597 }
    598 
    599 
    600 //++
    601 void  SphereThetaPhi::GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const
     565  // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...
     566 
     567  double omeg2pi= 1./(double)ntotpix;
     568  Omega_ = omeg2pi*DeuxPi;
     569 
     570  for(j=0; j <= m; j++)
     571    {
     572      Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi);
     573    }
     574}
     575
     576//++
     577template <class T>
     578void  SphereThetaPhi<T>::GetThetaSlice(int_4 index, r_4& theta, TVector<float>& phi, TVector<T>& value) const
    602579
    603580//    Retourne, pour la tranche en theta d'indice 'index' le theta
    … …  
    609586{
    610587  cout << "entree GetThetaSlice, couche no " << index << endl;
    611   if (index<0 || index > NbThetaSlices()) {
    612     // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    613     cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
    614   THROW(rangeCheckErr);
    615   }
    616   int_4 iring=Index(index,0);
    617   int_4 bid=this->NPhi(index);
    618   int lring=bid;
     588
     589  if(index < 0 || index > NbThetaSlices())
     590    {
     591      // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
     592      cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
     593      THROW(rangeCheckErr);
     594    }
     595
     596  int_4 iring= Index(index,0);
     597  int_4 bid  = this->NPhi(index);
     598  int lring  = bid;
    619599  cout << " iring= " << iring << " lring= " << lring << endl;
    620   phi.Realloc(lring);
    621   value.Realloc(lring);
    622   float T=0.;
    623   float F=0.;
    624   for (int kk=0; kk<lring;kk++) {
    625     PixThetaPhi(kk+iring,T,F);
    626     phi(kk)=F;
    627     value(kk)=PixVal(kk+iring);
    628   }
    629   theta=T;
    630 
    631 }
    632 
    633 
    634 /* --Methode-- */
    635 
     600
     601  phi.ReSize(lring);
     602  value.ReSize(lring);
     603  float Te= 0.;
     604  float Fi= 0.;
     605  for(int kk = 0; kk < lring; kk++)
     606    {
     607      PixThetaPhi(kk+iring,Te,Fi);
     608      phi(kk)= Fi;
     609      value(kk)= PixVal(kk+iring);
     610    }
     611  theta= Te;
     612}
     613
     614template <class T>
     615void SphereThetaPhi<T>::setmNPhi(int_4* array, int_4 m)
     616  //remplit  le tableau contenant le nombre de pixels en phi de chacune des bandes en theta
     617  //--
     618{
     619  NPhi_= new int_4[m];
     620  for(int k = 0; k < m; k++) NPhi_[k]= array[k];
     621}
     622
     623template <class T>
     624void SphereThetaPhi<T>::setmTNphi(int_4* array, int_4 m)
     625  //remplit  le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes
     626  //--
     627{
     628  TNphi_= new int_4[m];
     629  for(int k = 0; k < m; k++) TNphi_[k]= array[k];
     630}
     631
     632template <class T>
     633void SphereThetaPhi<T>::setmTheta(r_4* array, int_4 m)
     634  //remplit  le tableau contenant les valeurs limites de theta
     635  //--
     636{
     637  Theta_= new r_4[m];
     638  for(int k = 0; k < m; k++) Theta_[k]= array[k];
     639}
     640
     641template <class T>
     642void SphereThetaPhi<T>::setDataBlock(T* data, int_4 m)
     643  // remplit  le vecteur des contenus des pixels
     644{
     645  pixels_.FillFrom(m,data);
     646}
     647
     648template <class T>
     649void SphereThetaPhi<T>::print(ostream& os) const
     650{
     651  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
     652  //
     653  os << " NTheta_= " << NTheta_ << endl;
     654  os << " NPix_    = " << NPix_   << endl;
     655  os << " Omega_  =  " << Omega_   << endl;
     656
     657  os << " contenu de NPhi_ : ";
     658  for(int i=0; i < NTheta_; i++)
     659    {
     660      if(i%5 == 0) os << endl;
     661      os << NPhi_[i] <<", ";
     662    }
     663  os << endl;
     664
     665  os << " contenu de Theta_ : ";
     666  for(int i=0; i < NTheta_+1; i++)
     667    {
     668      if(i%5 == 0) os << endl;
     669      os << Theta_[i] <<", ";
     670    }
     671  os << endl;
     672
     673  os << " contenu de TNphi_ : ";
     674  for(int i=0; i < NTheta_+1; i++)
     675    {
     676      if(i%5 == 0) os << endl;
     677      os << TNphi_[i] <<", ";
     678    }
     679  os << endl;
     680
     681  os << " contenu de pixels : ";
     682  for(int i=0; i < NPix_; i++)
     683    {
     684      if(i%5 == 0) os << endl;
     685      os <<  pixels_(i) <<", ";
     686    }
     687  os << endl;
     688}
     689
     690///////////////////////////////////////////////////////////
     691// --------------------------------------------------------
     692//   Les objets delegues pour la gestion de persistance
     693// --------------------------------------------------------
     694//////////////////////////////////////////////////////////
     695
     696template <class T>
     697FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi()
     698{
     699  dobj= new SphereThetaPhi<T>;
     700  ownobj= true;
     701}
     702
     703template <class T>
     704FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename)
     705{
     706  dobj= new SphereThetaPhi<T>;
     707  ownobj= true;
     708  Read(filename);
     709}
     710
     711template <class T>
     712FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
     713{
     714  dobj= new SphereThetaPhi<T>(obj);
     715  ownobj= true;
     716}
     717
     718template <class T>
     719FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj)
     720{
     721  dobj= obj;
     722  ownobj= false;
     723}
     724
     725template <class T>
     726FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi()
     727{
     728  if (ownobj && dobj) delete dobj;
     729}
     730
     731template <class T>
     732AnyDataObj* FIO_SphereThetaPhi<T>::DataObj()
     733{
     734  return(dobj);
     735}
     736
     737template <class T>
     738void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is)
     739{
     740  cout << " FIO_SphereThetaPhi:: ReadSelf " << endl;
     741
     742  if(dobj == NULL)
     743    {
     744      dobj= new SphereThetaPhi<T>;
     745    }
     746
     747  // Pour savoir s'il y avait un DVList Info associe
     748  char strg[256];
     749  is.GetLine(strg, 255);
     750  bool hadinfo= false;
     751  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
     752  if(hadinfo)
     753    {    // Lecture eventuelle du DVList Info
     754      is >> dobj->Info();
     755    }
     756
     757  int_4 mNTheta;
     758  is.GetI4(mNTheta); 
     759  dobj->setSizeIndex(mNTheta);
     760
     761  int_4 mNPix;
     762  is.GetI4(mNPix);
     763  dobj->setNbPixels(mNPix);
     764
     765  r_8 mOmeg;
     766  is.GetR8(mOmeg);
     767  dobj->setPixSolAngle(mOmeg);
     768
     769  int_4* mNphi= new int_4[mNTheta];
     770  is.GetI4s(mNphi, mNTheta);
     771  dobj->setmNPhi(mNphi, mNTheta);
     772  delete [] mNphi;
     773
     774  int_4* mTNphi= new int_4[mNTheta+1];
     775  is.GetI4s(mTNphi, mNTheta+1);
     776  dobj->setmTNphi(mTNphi, mNTheta+1);
     777  delete [] mTNphi;
     778
     779  r_4* mTheta= new r_4[mNTheta+1];
     780  is.GetR4s(mTheta, mNTheta+1);
     781  dobj->setmTheta(mTheta, mNTheta+1);
     782  delete [] mTheta;
     783
     784  T* mPixels= new T[mNPix];
     785  //is.Get(mPixels, mNPix);
     786  PIOSReadArray(is, mPixels, mNPix);
     787  dobj->setDataBlock(mPixels, mNPix);
     788  delete [] mPixels;
     789}
     790
     791template <class T>
     792void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const
     793{
     794  cout << " FIO_SphereThetaPhi:: WriteSelf " << endl;
     795
     796  if(dobj == NULL)
     797    {
     798      cout << " WriteSelf:: dobj= null " << endl;
     799      return;
     800    }
     801
     802  char strg[256];
     803  int_4 mNTheta= dobj->SizeIndex();
     804  int_4 mNPix  = dobj->NbPixels();
     805 
     806  if(dobj->ptrInfo())
     807    {
     808      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo",mNTheta,mNPix);
     809      os.PutLine(strg);
     810      os << dobj->Info();
     811    }
     812  else
     813    {
     814      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d ",mNTheta,mNPix);
     815      os.PutLine(strg); 
     816    }
     817
     818  os.PutI4(mNTheta); 
     819  os.PutI4(mNPix);
     820  os.PutR8(dobj->PixSolAngle(0));
     821  os.PutI4s(dobj->getmNPhi() , mNTheta);
     822  os.PutI4s(dobj->getmTNphi(), mNTheta+1);
     823  os.PutR4s(dobj->getmTheta(), mNTheta+1);
     824  //os.Put((dobj->getDataBlock())->Data(), mNPix);
     825  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix);
     826}
Note: See TracChangeset for help on using the changeset viewer.