Changeset 726 in Sophya for trunk/SophyaLib/Samba


Ignore:
Timestamp:
Feb 21, 2000, 5:38:11 PM (26 years ago)
Author:
ansari
Message:

nouvelle gestion des thetaslices

Location:
trunk/SophyaLib/Samba
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/Samba/spherethetaphi.cc

    r701 r726  
    6868  NTheta_= s.NTheta_;
    6969  NPix_  = s.NPix_;
    70   NPhi_  = new int_4[NTheta_];
    71   Theta_ = new double[NTheta_+1];
    72   TNphi_ = new int_4[NTheta_+1];
    73   for(int k = 0; k < NTheta_; k++)
    74     {
    75       NPhi_[k] = s.NPhi_[k];
    76       Theta_[k]= s.Theta_[k];
    77       TNphi_[k]= s.TNphi_[k];
    78     }
    79   Theta_[NTheta_]= s.Theta_[NTheta_];
    80   TNphi_[NTheta_]= s.TNphi_[NTheta_];
     70  NPhi_.ReSize(NTheta_);
     71  Theta_.ReSize(NTheta_+1);
     72  TNphi_.ReSize(NTheta_+1);
     73
     74  NPhi_  = s.NPhi_;
     75  TNphi_ = s.TNphi_;
     76  Theta_ = s.Theta_;
    8177  Omega_ = s.Omega_;
    8278}
    … …  
    9086
    9187//--
    92 {
    93   Clear();
    94 }
     88{}
    9589
    9690//++
    … …  
    10094void SphereThetaPhi<T>::InitNul()
    10195//
    102 //    initialise à zéro les variables de classe pointeurs
    10396{
    10497  NTheta_= 0;
    10598  NPix_  = 0;
    106   Theta_ = NULL;
    107   NPhi_  = NULL;
    108   TNphi_ = NULL;
    10999//  pixels_.Reset();  Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
    110100}
    111101
    112 /* --Methode-- */
    113 template <class T>
    114 void SphereThetaPhi<T>::Clear()
    115 
    116 {
    117   if(Theta_)  delete[] Theta_;
    118   if(NPhi_ )  delete[] NPhi_;
    119   if(TNphi_)  delete[] TNphi_;
    120   InitNul();
    121 }
    122102
    123103//++
    … …  
    127107//--
    128108{
    129   Clear();
     109  InitNul();
    130110  Pixelize(m);
    131111}
    … …  
    152132  if((k < 0) || (k >= NPix_))
    153133    {
    154       //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    155       cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl;
    156       THROW(rangeCheckErr);
     134      throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
    157135    }
    158136  return pixels_(k);
    … …  
    168146  if((k < 0) || (k >= NPix_))
    169147    {
    170       cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
    171       //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    172       throw "SphereThetaPhi::PIxVal Pixel index out of range";
     148       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
    173149    }
    174150  return *(pixels_.Data()+k);
    … …  
    204180  // Theta_[kt] et Theta_[kt+1]
    205181  for( i=1; i< NTheta_; i++ )
    206     if( theta < Theta_[i] ) break;
    207  
    208   dphi= DeuxPi/(double)NPhi_[i-1];
    209  
    210   if (fgzn)  k= NPix_-TNphi_[i]+(int_4)(phi/dphi);
    211   else k= TNphi_[i-1]+(int_4)(phi/dphi);
     182    if( theta < Theta_(i) ) break;
     183 
     184  dphi= DeuxPi/(double)NPhi_(i-1);
     185 
     186  if (fgzn)  k= NPix_-TNphi_(i)+(int_4)(phi/dphi);
     187  else k= TNphi_(i-1)+(int_4)(phi/dphi);
    212188  return(k);
    213189}
    … …  
    229205  // recupère l'indice i de la tranche qui contient le pixel k
    230206  for( i=0; i< NTheta_; i++ )
    231     if( k < TNphi_[i+1] ) break;
     207    if( k < TNphi_(i+1) ) break;
    232208
    233209  // angle theta
    234   theta= 0.5*(Theta_[i]+Theta_[i+1]);
     210  theta= 0.5*(Theta_(i)+Theta_(i+1));
    235211  if (fgzn) theta= Pi-theta;
    236212 
    237213  // angle phi
    238   k -= TNphi_[i];
    239   phi= DeuxPi/(double)NPhi_[i]*(double)(k+.5);
     214  k -= TNphi_(i);
     215  phi= DeuxPi/(double)NPhi_(i)*(double)(k+.5);
    240216  if (fgzn) phi= DeuxPi-phi;
    241217}
    … …  
    290266  int i;
    291267  for( i=0; i< NTheta_; i++ )
    292     if(k < TNphi_[i+1]) break;
     268    if(k < TNphi_(i+1)) break;
    293269 
    294270  // valeurs limites de theta dans l'hemisphere Nord
    295   tetMin= Theta_[i];
    296   tetMax= Theta_[i+1];
     271  tetMin= Theta_(i);
     272  tetMax= Theta_(i+1);
    297273  // valeurs limites de theta dans l'hemisphere Sud
    298274  if (fgzn) {
    299     tetMin= Pi-Theta_[i+1];
    300     tetMax= Pi-Theta_[i];
     275    tetMin= Pi - Theta_(i+1);
     276    tetMax= Pi - Theta_(i);
    301277  }
    302278 
    303279  // pixel correspondant dans l'hemisphere Nord
    304   if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1;
     280  if (fgzn) k= TNphi_(i+1)-k+TNphi_(i)-1;
    305281 
    306282  // indice j de discretisation ( phi= j*dphi )
    307   j= k-TNphi_[i];
    308   dphi= DeuxPi/(double)NPhi_[i];
     283  j= k-TNphi_(i);
     284  dphi= DeuxPi/(double)NPhi_(i);
    309285 
    310286  // valeurs limites de phi
    … …  
    317293//++
    318294template <class T>
    319 int_4 SphereThetaPhi<T>::NbThetaSlices() const
     295uint_4 SphereThetaPhi<T>::NbThetaSlices() const
    320296
    321297//    Return number of theta-slices on the sphere
    322298//--
    323299{
     300  if (NTheta_<=0)
     301    {
     302      throw PException(" sphere not pixelized, NbSlice=0 ");
     303    }
    324304  return( 2*NTheta_);
    325305}
    … …  
    343323 
    344324  // nombre de pixels
    345   nbpix= NPhi_[kt];
     325  nbpix= NPhi_(kt);
    346326  return(nbpix);
    347327}
    … …  
    371351 
    372352  // valeurs limites de theta dans l'hemisphere Nord
    373   tetMin= Theta_[kt];
    374   tetMax= Theta_[kt+1];
     353  tetMin= Theta_(kt);
     354  tetMax= Theta_(kt+1);
    375355  // valeurs limites de theta dans l'hemisphere Sud
    376356  if (fgzn) {
    377     tetMin= Pi-Theta_[kt+1];
    378     tetMax= Pi-Theta_[kt];
     357    tetMin= Pi - Theta_(kt+1);
     358    tetMax= Pi - Theta_(kt);
    379359  } 
    380360}
    … …  
    399379 
    400380  // verifie si la tranche kt contient au moins jp pixels
    401   if( (jp< 0) || (jp >= NPhi_[kt]) ) {
     381  if( (jp< 0) || (jp >= NPhi_(kt)) ) {
    402382    phiMin= -88888.;
    403383    phiMax= -88888.;
    … …  
    405385  }
    406386 
    407   double dphi= DeuxPi/(double)NPhi_[kt];
     387  double dphi= DeuxPi/(double)NPhi_(kt);
    408388  phiMin= dphi*(double)(jp);
    409389  phiMax= dphi*(double)(jp+1);
    … …  
    429409 
    430410  // si la tranche kt contient au moins i pixels
    431   if( (jp>=0) && (jp<NPhi_[kt]) )
     411  if( (jp>=0) && (jp<NPhi_(kt)) )
    432412    {
    433413      // dans l'hemisphere Sud
    434       if (fgzn) k= NPix_-TNphi_[kt+1]+jp;
     414      if (fgzn) k= NPix_-TNphi_(kt+1)+jp;
    435415      // dans l'hemisphere Nord
    436       else k= TNphi_[kt]+jp;
     416      else k= TNphi_(kt)+jp;
    437417    }
    438418  else
    … …  
    463443  int i;
    464444  for(i = 0; i < NTheta_; i++)
    465     if(k < TNphi_[i+1]) break;
     445    if(k < TNphi_(i+1)) break;
    466446 
    467447  // indice  kt de tranche
    … …  
    470450 
    471451  // indice jp de pixel
    472   if (fgzn) jp= TNphi_[i+1]-k-1;
    473   else jp= k-TNphi_[i]; 
     452  if (fgzn) jp= TNphi_(i+1)-k-1;
     453  else jp= k-TNphi_(i); 
    474454}
    475455//++
    … …  
    502482  // Creation des vecteurs contenant :
    503483  // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
    504   Theta_= new double[m+1];
     484  //  Theta_= new double[m+1];
     485  Theta_.ReSize(m+1);
    505486 
    506487  // Le nombre de pixels en phi de chacune des bandes en theta
    507   NPhi_ = new int_4[m];
     488  //  NPhi_ = new int_4[m];
     489  NPhi_.ReSize(m);
    508490 
    509491  // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
    510492  // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
    511   TNphi_= new int_4[m+1];
     493  //  TNphi_= new int_4[m+1];
     494    TNphi_.ReSize(m+1);
    512495 
    513496  // Calcul du nombre total de pixels dans chaque bande pour optimiser
    … …  
    515498 
    516499  //calotte polaire
    517   TNphi_[0]= 0;
    518   NPhi_[0] = 1;
     500  TNphi_(0)= 0;
     501  NPhi_(0) = 1;
    519502 
    520503  //bandes jusqu'a l'equateur
    521504  for(j = 1; j < m; j++)
    522505    {
    523       TNphi_[j]= TNphi_[j-1]+NPhi_[j-1];
    524       NPhi_[j] = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
     506      TNphi_(j)= TNphi_(j-1)+NPhi_(j-1);
     507      NPhi_(j) = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
    525508    }
    526509 
    527510  // Nombre total de pixels sur l'hemisphere
    528   ntotpix  = TNphi_[m-1]+NPhi_[m-1];
    529   TNphi_[m]= ntotpix;
     511  ntotpix  = TNphi_(m-1)+NPhi_(m-1);
     512  TNphi_(m)= ntotpix;
    530513  // et sur la sphere entiere
    531514  NPix_= 2*ntotpix;
    … …  
    548531  for(j=0; j <= m; j++)
    549532    {
    550       Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi);
     533      Theta_(j)= acos(1.-(double)TNphi_(j)*omeg2pi);
    551534    }
    552535}
    … …  
    566549  if(index < 0 || index >= NbThetaSlices())
    567550    {
    568       // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    569       cout << " SphereThetaPhi::GetThetaSlice : Pixel index out of range" <<endl;
    570       THROW(rangeCheckErr);
     551       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
    571552    }
    572553
    … …  
    587568}
    588569
    589 template <class T>
    590 void SphereThetaPhi<T>::setmNPhi(int_4* array, int_4 m)
    591   //remplit le tableau contenant le nombre de pixels en phi de chacune des bandes en theta
    592   //
    593 {
    594   NPhi_= new int_4[m];
    595   for(int k = 0; k < m; k++) NPhi_[k]= array[k];
    596 }
    597 
    598 template <class T>
    599 void SphereThetaPhi<T>::setmTNphi(int_4* array, int_4 m)
    600   //remplit  le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes
    601   //
    602 {
    603   TNphi_= new int_4[m];
    604   for(int k = 0; k < m; k++) TNphi_[k]= array[k];
    605 }
    606 
    607 template <class T>
    608 void SphereThetaPhi<T>::setmTheta(double* array, int_4 m)
    609   //remplit  le tableau contenant les valeurs limites de theta
    610   //
    611 {
    612   Theta_= new double[m];
    613   for(int k = 0; k < m; k++) Theta_[k]= array[k];
    614 }
     570//++
     571template <class T>
     572void SphereThetaPhi<T>::GetThetaSlice(int_4 index,double& theta, double& phi0,TVector<int>& pixelIndices, TVector<T>& value) const
     573
     574//    For a theta-slice with index 'index', return :
     575//    the corresponding "theta"
     576//    the corresponding "phi" for first pixel of the slice
     577//    a vector containing the indices of the pixels of the slice
     578//    (equally distributed in phi)
     579//    a vector containing the corresponding values of pixels
     580//--
     581
     582{
     583
     584  if(index < 0 || index >= NbThetaSlices())
     585    {
     586       throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
     587 
     588    }
     589
     590  int iring= Index(index,0);
     591  int lring  = NPhi(index);
     592
     593  pixelIndices.ReSize(lring);
     594  value.ReSize(lring);
     595  double Te= 0.;
     596  double Fi= 0.;
     597  for(int kk = 0; kk < lring; kk++)
     598    {
     599      pixelIndices(kk)=kk+iring ;
     600      value(kk)= PixVal(kk+iring);
     601    }
     602      PixThetaPhi(iring,theta,phi0);
     603}
     604
     605
    615606
    616607
    … …  
    629620    {
    630621      if(i%5 == 0) os << endl;
    631       os << NPhi_[i] <<", ";
     622      os << NPhi_(i) <<", ";
    632623    }
    633624  os << endl;
    … …  
    637628    {
    638629      if(i%5 == 0) os << endl;
    639       os << Theta_[i] <<", ";
     630      os << Theta_(i) <<", ";
    640631    }
    641632  os << endl;
    … …  
    645636    {
    646637      if(i%5 == 0) os << endl;
    647       os << TNphi_[i] <<", ";
     638      os << TNphi_(i) <<", ";
    648639    }
    649640  os << endl;
    … …  
    675666{
    676667  dobj= new SphereThetaPhi<T>;
    677   //  dobj->DataBlock().SetTemp(true);
    678668  dobj->pixels_.SetTemp(true);
    679669  ownobj= true;
    … …  
    685675{
    686676  dobj= new SphereThetaPhi<T>(obj, true);
    687   //  dobj->DataBlock().SetTemp(true);
    688677  dobj->pixels_.SetTemp(true);
    689678  ownobj= true;
    … …  
    716705    {
    717706      dobj= new SphereThetaPhi<T>;
    718       //      dobj->DataBlock().SetTemp(true);
    719707      dobj->pixels_.SetTemp(true);
    720708      ownobj= true;     
    721709    }
    722710
     711  // On lit les 3 premiers uint_8
     712  uint_8 itab[3];
     713  is.Get(itab, 3);
    723714// Let's Read the SphereCoordSys object  -- ATTENTIOn - $CHECK$
    724715  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
    … …  
    742733  is.GetR8(mOmeg);
    743734  dobj->setParameters(mNPix, mOmeg, mNTheta);
    744 
    745   int_4* mNphi= new int_4[mNTheta];
    746   is.GetI4s(mNphi, mNTheta);
    747   dobj->setmNPhi(mNphi, mNTheta);
    748   delete [] mNphi;
    749 
    750   int_4* mTNphi= new int_4[mNTheta+1];
    751   is.GetI4s(mTNphi, mNTheta+1);
    752   dobj->setmTNphi(mTNphi, mNTheta+1);
    753   delete [] mTNphi;
    754 
    755   double* mTheta= new double[mNTheta+1];
    756   is.GetR8s(mTheta, mNTheta+1);
    757   dobj->setmTheta(mTheta, mNTheta+1);
    758   delete [] mTheta;
    759 
    760 // On lit le DataBlock;
    761   //FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
     735  FIO_NDataBlock<int_4> fio_nd_nphi(&dobj->NPhi_);
     736  fio_nd_nphi.Read(is);
     737  FIO_NDataBlock<int_4> fio_nd_Tnphi(&dobj->TNphi_);
     738  fio_nd_Tnphi.Read(is);
     739  FIO_NDataBlock<r_8> fio_nd_Theta(&dobj->Theta_);
     740  fio_nd_Theta.Read(is);
     741
    762742FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
    763743fio_nd.Read(is);
    764 //  is >> dobj->DataBlock();
    765744}
    766745
    … …  
    775754    }
    776755
     756//  On ecrit 3 uint_8
     757//  0 : Numero de version,  1 : Size index,  2  reserve a l
     758uint_8 itab[3];
     759itab[0] = 1;
     760itab[1] = dobj->SizeIndex();
     761itab[2] = 0;
     762os.Put(itab, 3);
    777763// Let's write the SphereCoordSys object
    778764  FIO_SphereCoordSys fio_scs( dobj->GetCoordSys());
    … …  
    798784  os.PutI4(mNPix);
    799785  os.PutR8(dobj->PixSolAngle());
    800   os.PutI4s(dobj->getmNPhi() , mNTheta);
    801   os.PutI4s(dobj->getmTNphi(), mNTheta+1);
    802   os.PutR8s(dobj->getmTheta(), mNTheta+1);
    803 // On ecrit le datablock
    804   //FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
     786  FIO_NDataBlock<int_4> fio_nd_nphi(&dobj->NPhi_);
     787  fio_nd_nphi.Write(os);
     788  FIO_NDataBlock<int_4> fio_nd_Tnphi(&dobj->TNphi_);
     789  fio_nd_Tnphi.Write(os);
     790  FIO_NDataBlock<r_8> fio_nd_Theta(&dobj->Theta_);
     791  fio_nd_Theta.Write(os);
    805792FIO_NDataBlock<T> fio_nd(&dobj->pixels_);
    806793fio_nd.Write(os);
    807 //  os << dobj->DataBlock();
    808794}
    809795

  • trunk/SophyaLib/Samba/spherethetaphi.h

    r708 r726  
    4141virtual ~SphereThetaPhi();
    4242
     43/*! Setting blockdata to temporary (see ndatablock documentation) */
     44inline virtual void SetTemp(bool temp=false) const {pixels_.SetTemp(temp);};
     45
    4346// ------------ Definition of PixelMap abstract methods -
    4447
    … …  
    9497/* Nombre de tranches en theta */
    9598/*!    Return number of theta-slices on the sphere */
    96 int_4 NbThetaSlices() const;
     99uint_4 NbThetaSlices() const;
    97100
    98101/* Nombre de pixels en phi de la tranche d'indice kt */
    … …  
    140143void GetThetaSlice(int_4 index,double& theta,TVector<double>& phi,TVector<T>& value) const;
    141144
    142 
     145/*!   For a theta-slice with index 'index', return :
     146
     147   the corresponding "theta"
     148
     149   the corresponding "phi" for first pixel of the slice
     150
     151    a vector containing indices of the pixels of the slice
     152
     153   (equally distributed in phi)
     154
     155    a vector containing the corresponding values of pixels
     156*/
     157 void GetThetaSlice(int_4 index,double& theta, double& phi0,TVector<int>& pixelIndices, TVector<T>& value) const ;
    143158
    144159
    … …  
    150165// ------------- méthodes internes ----------------------         
    151166void InitNul();
    152 void Clear();
    153167inline void setParameters(int nbpix, double omega, int nbThetaIndex)
    154168  {
    … …  
    157171    NTheta_= nbThetaIndex;
    158172  }
    159 void setmNPhi(int* array, int m);
    160 void setmTNphi(int* array, int m);
    161 void setmTheta(double* array, int m);
    162  /*retourne l'adresse du tableau contenant le nombre de pixels en phi de chacune des bandes en theta */
    163 inline const int* getmNPhi() const { return(NPhi_); }
    164 
    165 /* retourne le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes */
    166 inline const int* getmTNphi() const { return(TNphi_); }
    167 
    168 /* retourne le tableau contenant les valeurs limites de theta */
    169 inline const double* getmTheta() const { return(Theta_); }
    170173
    171174// ------------- variables internes ---------------------
    … …  
    173176 int_4 NPix_;     // nombre total de pixels
    174177 double Omega_; // angle solide constant pour chaque pixel
    175  int_4* NPhi_;    // tableau donnant, pour chaque bande en theta, le nombre de
    176                 // pixels selon phi
    177 int_4* TNphi_;
    178 double* Theta_;
    179 NDataBlock<T> pixels_;
     178 NDataBlock<int_4>  NPhi_;  // tableau donnant, pour chaque bande en theta,
     179                         //le nombre de pixels selon phi
     180 NDataBlock<int_4> TNphi_;
     181 NDataBlock<r_8> Theta_;
     182 NDataBlock<T> pixels_;
    180183};
    181184
Note: See TracChangeset for help on using the changeset viewer.