Changeset 515 in Sophya for trunk/SophyaLib


Ignore:
Timestamp:
Oct 26, 1999, 2:44:44 PM (26 years ago)
Author:
ansari
Message:

Portage SGI-CC Reza 26/10/99

Location:
trunk/SophyaLib
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/SophyaLib/NTools/dvlist.cc

    r514 r515  
    6464
    6565/* --Methode-- */
    66 DVList::DVList(DVList& dvl)
     66DVList::DVList(const DVList& dvl)
    6767{
    6868Merge(dvl);
  • trunk/SophyaLib/NTools/dvlist.h

    r490 r515  
    6868
    6969                    DVList();
    70                     DVList(DVList&);
     70                    DVList(const DVList&);
    7171                    DVList(char *flnm);
    7272
  • trunk/SophyaLib/NTools/fftserver.cc

    r514 r515  
    9292  checkint_cfft(l);
    9393  float* foo = new float[2*l];
    94   for (int i=0;i<l;i++){
     94  int i;
     95  for (i=0;i<l;i++){
    9596    foo[2*i]=inout[i].real();
    9697    foo[2*i+1]=inout[i].imag();
     
    9899  cfftf_(&l, foo, ws_cfft);
    99100  inout[0]=complex<float> (foo[0],foo[1]);
    100   for (int i=1;i<l;i++) inout[l-i]= complex<float> (foo[2*i], foo[2*i+1]);
     101  for (i=1;i<l;i++) inout[l-i]= complex<float> (foo[2*i], foo[2*i+1]);
    101102  delete[] foo;
    102103}
     
    106107  checkint_cdfft(l);
    107108  double* foo=new double[2*l];
    108   for (int i=0;i<l;i++){
     109  int i;
     110  for (i=0;i<l;i++){
    109111    foo[2*i]=inout[i].real();
    110112    foo[2*i+1]=inout[i].imag();
     
    112114  cdfftf_(&l, foo, ws_cdfft);
    113115  inout[0]=complex<double> (foo[0],foo[1]);
    114   for (int i=1;i<l;i++) {
     116  for (i=1;i<l;i++) {
    115117    inout[l-i]= complex<double> (foo[2*i],foo[2*i+1]);
    116118  }
     
    134136  checkint_cfft(l);
    135137  float* foo = new float[2*l];
    136   for (int i=0;i<l;i++){
     138  int i;
     139  for (i=0;i<l;i++){
    137140    foo[2*i]=inout[i].real();
    138141    foo[2*i+1]=inout[i].imag();
    139142  }
    140143  cfftf_(&l, foo, ws_cfft);
    141   for (int i=0;i<l;i++) inout[i]=complex<float> (foo[2*i],foo[2*i+1]);
     144  for (i=0;i<l;i++) inout[i]=complex<float> (foo[2*i],foo[2*i+1]);
    142145  delete[] foo;
    143146}
     
    147150  checkint_cdfft(l);
    148151  double* foo = new double[2*l];
    149   for (int i=0;i<l;i++){
     152  int i;
     153  for (i=0;i<l;i++){
    150154    foo[2*i]=inout[i].real();
    151155    foo[2*i+1]=inout[i].imag();
    152156  }
    153157  cdfftf_(&l, foo, ws_cdfft);
    154   for (int i=0;i<l;i++) inout[i]=complex<double> (foo[2*i],foo[2*i+1]);
     158  for (i=0;i<l;i++) inout[i]=complex<double> (foo[2*i],foo[2*i+1]);
    155159  delete[] foo;
    156160}
  • trunk/SophyaLib/Samba/scan.cc

    r276 r515  
    116116  PhiZero_ = s. PhiZero_;
    117117  sPix_=new r_8[ NmaxPts_];
    118   for (int_4 k=0; k<NmaxPts_; k++)  sPix_[k]=s.sPix_[k];
    119   for (int k=0; k<9; k++) Rota_[k]=s. Rota_[k];
     118  int_4 k;
     119  for (k=0; k<NmaxPts_; k++)  sPix_[k]=s.sPix_[k];
     120  for (k=0; k<9; k++) Rota_[k]=s. Rota_[k];
    120121}
    121122//++
  • trunk/SophyaLib/Samba/spheregorski.cc

    r487 r515  
    11#include "spheregorski.h"
    22#include "strutil.h"
     3#include <math.h>
    34#include <complex>
    45#include "piocmplx.h"
     
    886887  iy_hi  =     iy/128;
    887888  ipf= (PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
    888   ipf = ipf / pow(ns_max/nside,2);//  ! in {0, nside**2 - 1}
    889   return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
     889  //  ipf = ipf / pow(ns_max/nside,2.);//  ! in {0, nside**2 - 1}
     890  //  return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
     891  // $CHECK$  Reza 25/10/99 , pow remplace par *
     892  //  ipf = ipf / ((ns_max/nside)*(ns_max/nside)); //  ! in {0, nside**2 - 1}
     893  //  return ( ipf + face_num*(nside*nside);//    ! in {0, 12*nside**2 - 1}
    890894}
    891895
  • trunk/SophyaLib/Samba/spherethetaphi.cc

    r487 r515  
    44#include "piocmplx.h"
    55#include <iostream.h>
     6
     7
     8//***************************************************************
     9//++
     10// Class        SphereThetaPhi
     11//
     12//
     13// include      spherethetaphi.h nbmath.h
     14//
     15//    Découpage de la sphère selon theta et phi, chaque
     16//    hémisphère étant découpé en (m-1) parallèles (l'équateur compte pour du
     17//    beurre), chacune des m bandes de theta ainsi définies étant découpée par
     18//    des méridiens équirepartis, ce découpage étant fait de sorte que tous
     19//    les pixels aient la même surface et soient le plus carré possible.
     20//    On commence par découper l'hémisphère de z positif en partant du pôle et
     21//    en allant vers l'équateur. Le premier pixel est la calotte polaire,
     22//    il est circulaire et centré sur theta=0.
     23//--
     24//++
     25//
     26// Links        Parents
     27//
     28//    SphericalMap
     29//--
     30//++
     31//
     32// Links        Descendants
     33//
     34//     
     35//--
     36
     37/* --Methode-- */
     38//++
     39// Titre        Constructeurs
     40//--
     41//++
     42
     43template <class T>
     44SphereThetaPhi<T>::SphereThetaPhi()
     45
     46//--
     47{
     48  InitNul();
     49}
     50
     51
     52/* --Methode-- */
     53
     54//++
     55template <class T>
     56SphereThetaPhi<T>::SphereThetaPhi(int m)
     57
     58//    Constructeur : m est le nombre de bandes en theta sur un hémisphère
     59//    (la calotte constituant la premiere bande).
     60//    pet est le nombre de pixels (pétales) de la bande en contact avec la
     61//    calotte polaire. Pour l'instant pet est inopérant!
     62//--
     63{
     64  InitNul();
     65  Pixelize(m);
     66}
     67
     68template <class T>
     69SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)
     70{
     71  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
     72  NTheta_= s.NTheta_;
     73  NPix_  = s.NPix_;
     74  NPhi_  = new int[NTheta_];
     75  Theta_ = new double[NTheta_+1];
     76  TNphi_ = new int[NTheta_+1];
     77  for(int k = 0; k < NTheta_; k++)
     78    {
     79      NPhi_[k] = s.NPhi_[k];
     80      Theta_[k]= s.Theta_[k];
     81      TNphi_[k]= s.TNphi_[k];
     82    }
     83  Theta_[NTheta_]= s.Theta_[NTheta_];
     84  TNphi_[NTheta_]= s.TNphi_[NTheta_];
     85  Omega_ = s.Omega_;
     86  pixels_= s.pixels_;
     87}
     88
     89//++
     90// Titre        Destructeur
     91//--
     92//++
     93template <class T>
     94SphereThetaPhi<T>::~SphereThetaPhi()
     95
     96//--
     97{
     98  Clear();
     99}
     100
     101//++
     102// Titre        Méthodes
     103//--
     104template <class T>
     105void SphereThetaPhi<T>::InitNul()
     106//
     107//    initialise à zéro les variables de classe pointeurs
     108{
     109  NTheta_= 0;
     110  NPix_  = 0;
     111  Theta_ = NULL;
     112  NPhi_  = NULL;
     113  TNphi_ = NULL;
     114  pixels_.Reset();
     115}
     116
     117/* --Methode-- */
     118template <class T>
     119void SphereThetaPhi<T>::Clear()
     120
     121{
     122  if(Theta_)  delete[] Theta_;
     123  if(NPhi_ )  delete[] NPhi_;
     124  if(TNphi_)  delete[] TNphi_;
     125  InitNul();
     126}
     127
     128//++
     129template <class T>
     130void SphereThetaPhi<T>::Resize(int m)
     131//   re-pixelize the sphere
     132//--
     133{
     134  Clear();
     135  Pixelize(m);
     136}
     137
     138/* --Methode-- */
     139//++
     140template <class T>
     141int SphereThetaPhi<T>::NbPixels() const
     142
     143//    Retourne le nombre de pixels du découpage
     144//--
     145{
     146  return(NPix_);
     147}
     148
     149/* --Methode-- */
     150//++
     151template <class T>
     152T& SphereThetaPhi<T>::PixVal(int k)
     153
     154//    Retourne la valeur du contenu du pixel d'indice k
     155//--
     156{
     157  if((k < 0) || (k >= NPix_))
     158    {
     159      //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
     160      cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl;
     161      THROW(rangeCheckErr);
     162    }
     163  return pixels_(k);
     164}
     165
     166//++
     167template <class T>
     168T const& SphereThetaPhi<T>::PixVal(int k) const
     169
     170//    Retourne la valeur du contenu du pixel d'indice k
     171//--
     172{
     173  if((k < 0) || (k >= NPix_))
     174    {
     175      cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
     176      //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
     177      throw "SphereThetaPhi::PIxVal Pixel index out of range";
     178    }
     179  return *(pixels_.Data()+k);
     180}
     181
     182/* --Methode-- */
     183//++
     184template <class T>
     185int SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const
     186
     187//  Retourne l'indice du pixel vers lequel pointe une direction définie par
     188//    ses coordonnées sphériques
     189//--
     190{
     191  double dphi;
     192  int i,j,k;
     193  bool fgzn = false;
     194
     195  if((theta > Pi) || (theta < 0.)) return(-1);
     196  if((phi < 0.) || (phi > DeuxPi)) return(-1);
     197  if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;}
     198 
     199  // La bande d'indice kt est limitée par les valeurs de theta contenues dans
     200  // Theta_[kt] et Theta_[kt+1]
     201  for( i=1; i< NTheta_; i++ )
     202    if( theta < Theta_[i] ) break;
     203 
     204  dphi= DeuxPi/(double)NPhi_[i-1];
     205 
     206  if (fgzn)  k= NPix_-TNphi_[i]+(int)(phi/dphi);
     207  else k= TNphi_[i-1]+(int)(phi/dphi);
     208  return(k);
     209}
     210
     211/* --Methode-- */
     212//++
     213template <class T>
     214void SphereThetaPhi<T>::PixThetaPhi(int k,double& theta,double& phi) const
     215
     216//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
     217//--
     218{
     219  int i;
     220  bool fgzn = false;
     221 
     222  if((k < 0) || (k >= NPix_)) {theta = -99999.; phi = -99999.; return; }
     223  if( k >= NPix_/2)  {fgzn = true; k = NPix_-1-k;}
     224
     225  // recupère l'indice i de la tranche qui contient le pixel k
     226  for( i=0; i< NTheta_; i++ )
     227    if( k < TNphi_[i+1] ) break;
     228
     229  // angle theta
     230  theta= 0.5*(Theta_[i]+Theta_[i+1]);
     231  if (fgzn) theta= Pi-theta;
     232 
     233  // angle phi
     234  k -= TNphi_[i];
     235  phi= DeuxPi/(double)NPhi_[i]*(double)(k+.5);
     236  if (fgzn) phi= DeuxPi-phi;
     237}
     238
     239//++
     240template <class T>
     241double SphereThetaPhi<T>::PixSolAngle(int dummy) const
     242
     243//    Pixel Solid angle  (steradians)
     244//    All the pixels have the same solid angle. The dummy argument is
     245//    for compatibility with eventual pixelizations which would not
     246//    fulfil this requirement.
     247//--
     248{
     249  return Omega_;
     250}
     251
     252/* --Methode-- */
     253//++
     254template <class T>
     255void SphereThetaPhi<T>::Limits(int k,double& tetMin,double& tetMax,double& phiMin,double& phiMax)
     256
     257//    Retourne les valeurs de theta et phi limitant le pixel d'indice k
     258//--
     259  {
     260  int j;
     261  double dphi;
     262  bool fgzn= false;
     263 
     264  if((k < 0) || (k >= NPix_)) {
     265    tetMin= -99999.;
     266    phiMin= -99999.; 
     267    tetMax= -99999.;
     268    phiMax= -99999.; 
     269    return;
     270  }
     271 
     272  // si on se trouve dans l'hémisphère Sud
     273  if(k >= NPix_/2) {
     274    fgzn= true;
     275    k= NPix_-1-k;
     276  }
     277 
     278  // on recupere l'indice i de la tranche qui contient le pixel k
     279  int i;
     280  for( i=0; i< NTheta_; i++ )
     281    if(k < TNphi_[i+1]) break;
     282 
     283  // valeurs limites de theta dans l'hemisphere Nord
     284  tetMin= Theta_[i];
     285  tetMax= Theta_[i+1];
     286  // valeurs limites de theta dans l'hemisphere Sud
     287  if (fgzn) {
     288    tetMin= Pi-Theta_[i+1];
     289    tetMax= Pi-Theta_[i];
     290  }
     291 
     292  // pixel correspondant dans l'hemisphere Nord
     293  if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1;
     294 
     295  // indice j de discretisation ( phi= j*dphi )
     296  j= k-TNphi_[i];
     297  dphi= DeuxPi/(double)NPhi_[i];
     298 
     299  // valeurs limites de phi
     300  phiMin= dphi*(double)(j);
     301  phiMax= dphi*(double)(j+1);
     302  return;
     303}
     304
     305/* --Methode-- */
     306//++
     307template <class T>
     308int SphereThetaPhi<T>::NbThetaSlices() const
     309
     310//    Retourne le nombre de tranches en theta sur la sphere
     311//--
     312{
     313  int nbslices;
     314  nbslices= 2*NTheta_;
     315  return(nbslices);
     316}
     317
     318/* --Methode-- */
     319//++
     320template <class T>
     321int SphereThetaPhi<T>::NPhi(int kt) const
     322
     323//    Retourne le nombre de pixels en phi de la tranche kt
     324//--
     325{
     326  int nbpix; 
     327  // verification
     328  if((kt < 0) || (kt >= 2*NTheta_)) return(-1);
     329 
     330  // si on se trouve dans l'hemisphere Sud
     331  if(kt >= NTheta_) {
     332    kt= 2*NTheta_-1-kt;
     333  }
     334 
     335  // nombre de pixels
     336  nbpix= NPhi_[kt];
     337  return(nbpix);
     338}
     339
     340
     341/* --Methode-- */
     342//++
     343template <class T>
     344void SphereThetaPhi<T>::Theta(int kt,double& tetMin,double& tetMax)
     345
     346//    Retourne les valeurs de theta limitant la tranche kt
     347//--
     348{
     349  bool fgzn= false;
     350  // verification
     351  if( (kt< 0) || (kt>= 2*NTheta_) ) {
     352    tetMin= -99999.;
     353    tetMax= -99999.;
     354    return;
     355  }
     356
     357  // si on se trouve dans l'hemisphere Sud
     358  if( kt >= NTheta_ ) {
     359    fgzn= true;
     360    kt= 2*NTheta_-1-kt;
     361  }
     362 
     363  // valeurs limites de theta dans l'hemisphere Nord
     364  tetMin= Theta_[kt];
     365  tetMax= Theta_[kt+1];
     366  // valeurs limites de theta dans l'hemisphere Sud
     367  if (fgzn) {
     368    tetMin= Pi-Theta_[kt+1];
     369    tetMax= Pi-Theta_[kt];
     370  } 
     371}
     372
     373/* --Methode-- */
     374//++
     375template <class T>
     376void SphereThetaPhi<T>::Phi(int kt,int jp,double& phiMin,double& phiMax)
     377
     378//    Retourne les valeurs de phi limitant le pixel jp de la tranche kt
     379//--
     380{
     381  // verification
     382  if((kt < 0) || (kt >= 2*NTheta_)) {
     383    phiMin= -99999.;
     384    phiMax= -99999.;
     385    return;
     386  }
     387 
     388  // si on se trouve dans l'hemisphere Sud
     389  if(kt >= NTheta_) kt= 2*NTheta_-1-kt;
     390 
     391  // verifie si la tranche kt contient au moins jp pixels
     392  if( (jp< 0) || (jp >= NPhi_[kt]) ) {
     393    phiMin= -88888.;
     394    phiMax= -88888.;
     395    return;
     396  }
     397 
     398  double dphi= DeuxPi/(double)NPhi_[kt];
     399  phiMin= dphi*(double)(jp);
     400  phiMax= dphi*(double)(jp+1);
     401  return;
     402}
     403
     404/* --Methode-- */
     405//++
     406template <class T>
     407int SphereThetaPhi<T>::Index(int kt,int jp) const
     408
     409//    Retourne l'indice du pixel d'indice jp dans la tranche kt
     410//--
     411{
     412  int k;
     413  bool fgzn= false;
     414 
     415  // si on se trouve dans l'hemisphere Sud
     416  if(kt >= NTheta_) {
     417    fgzn= true;
     418    kt= 2*NTheta_-1-kt;
     419  }
     420 
     421  // si la tranche kt contient au moins i pixels
     422  if( (jp>=0) && (jp<NPhi_[kt]) )
     423    {
     424      // dans l'hemisphere Sud
     425      if (fgzn) k= NPix_-TNphi_[kt+1]+jp;
     426      // dans l'hemisphere Nord
     427      else k= TNphi_[kt]+jp;
     428    }
     429  else
     430    {
     431      k= 9999;
     432      printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
     433    }
     434  return(k);
     435}
     436
     437/* --Methode-- */
     438//++
     439template <class T>
     440void SphereThetaPhi<T>::ThetaPhiIndex(int k,int& kt,int& jp)
     441
     442//    Retourne les indices kt et jp du pixel d'indice k
     443//--
     444{
     445  bool fgzn= false; 
     446  // si on se trouve dans l'hemisphere Sud
     447  if(k >= NPix_/2)
     448    {
     449      fgzn= true;
     450      k= NPix_-1-k;
     451    }
     452 
     453  // on recupere l'indice kt de la tranche qui contient le pixel k
     454  int i;
     455  for(i = 0; i < NTheta_; i++)
     456    if(k < TNphi_[i+1]) break;
     457 
     458  // indice  kt de tranche
     459  if (fgzn) kt= 2*NTheta_-1-i;
     460  else kt= i;
     461 
     462  // indice jp de pixel
     463  if (fgzn) jp= TNphi_[i+1]-k-1;
     464  else jp= k-TNphi_[i]; 
     465}
     466//++
     467template <class T>
     468void  SphereThetaPhi<T>::Pixelize(int m)
     469
     470//    effectue le découpage en pixels (m et pet ont la même signification
     471//    que pour le constructeur)
     472//
     473//    Chaque bande de theta sera découpée en partant de phi=0 ...
     474//    L'autre hémisphère est parcourue dans le même sens en phi et de
     475//    l'équateur vers le pôle (le pixel qui suit le dernier de la bande la plus
     476//    proche de l'équateur a z>0 est celui de plus petit phi de la bande la
     477//    plus proche de l'equateur a z<0).
     478//--
     479{
     480  int ntotpix,i,j;
     481 
     482  // Decodage et controle des arguments d'appel :
     483  // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
     484  if (m < 2) m = 2;
     485  if (m > 16384) m = 16384;
     486 
     487  // On memorise les arguments d'appel
     488  NTheta_ = m; 
     489 
     490  // On commence par decouper l'hemisphere z>0.
     491  // Creation des vecteurs contenant :
     492  // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
     493  Theta_= new double[m+1];
     494 
     495  // Le nombre de pixels en phi de chacune des bandes en theta
     496  NPhi_ = new int[m];
     497 
     498  // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
     499  // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
     500  TNphi_= new int[m+1];
     501 
     502  // Calcul du nombre total de pixels dans chaque bande pour optimiser
     503  // le rapport largeur/hauteur des pixels
     504 
     505  //calotte polaire
     506  TNphi_[0]= 0;
     507  NPhi_[0] = 1;
     508 
     509  //bandes jusqu'a l'equateur
     510  for(j = 1; j < m; j++)
     511    {
     512      TNphi_[j]= TNphi_[j-1]+NPhi_[j-1];
     513      NPhi_[j] = (int)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
     514    }
     515 
     516  // Nombre total de pixels sur l'hemisphere
     517  ntotpix  = TNphi_[m-1]+NPhi_[m-1];
     518  TNphi_[m]= ntotpix;
     519  // et sur la sphere entiere
     520  NPix_= 2*ntotpix;
     521 
     522  // Creation et initialisation du vecteur des contenus des pixels
     523  pixels_.ReSize(NPix_);
     524  pixels_.Reset();
     525
     526  // Determination des limites des bandes en theta :
     527  // omeg est l'angle solide couvert par chaque pixel,
     528  // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg
     529  // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide
     530  //de la
     531  // calotte allant du pole a la limite haute de la bande kt vaut
     532  // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...
     533 
     534  double omeg2pi= 1./(double)ntotpix;
     535  Omega_ = omeg2pi*DeuxPi;
     536 
     537  for(j=0; j <= m; j++)
     538    {
     539      Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi);
     540    }
     541}
     542
     543//++
     544template <class T>
     545void SphereThetaPhi<T>::GetThetaSlice(int index,double& theta, TVector<double>& phi, TVector<T>& value) const
     546
     547//    Retourne, pour la tranche en theta d'indice 'index' le theta
     548//    correspondant, un vecteur (Peida) contenant les phi des pixels de
     549//    la tranche, un vecteur (Peida) contenant les valeurs de pixel
     550//    correspondantes
     551//--
     552
     553{
     554  cout << "entree GetThetaSlice, couche no " << index << endl;
     555
     556  if(index < 0 || index > NbThetaSlices())
     557    {
     558      // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
     559      cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
     560      THROW(rangeCheckErr);
     561    }
     562
     563  int iring= Index(index,0);
     564  int bid  = this->NPhi(index);
     565  int lring  = bid;
     566  cout << " iring= " << iring << " lring= " << lring << endl;
     567
     568  phi.ReSize(lring);
     569  value.ReSize(lring);
     570  double Te= 0.;
     571  double Fi= 0.;
     572  for(int kk = 0; kk < lring; kk++)
     573    {
     574      PixThetaPhi(kk+iring,Te,Fi);
     575      phi(kk)= Fi;
     576      value(kk)= PixVal(kk+iring);
     577    }
     578  theta= Te;
     579}
     580
     581template <class T>
     582void SphereThetaPhi<T>::setmNPhi(int* array, int m)
     583  //remplit le tableau contenant le nombre de pixels en phi de chacune des bandes en theta
     584  //--
     585{
     586  NPhi_= new int[m];
     587  for(int k = 0; k < m; k++) NPhi_[k]= array[k];
     588}
     589
     590template <class T>
     591void SphereThetaPhi<T>::setmTNphi(int* array, int m)
     592  //remplit  le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes
     593  //--
     594{
     595  TNphi_= new int[m];
     596  for(int k = 0; k < m; k++) TNphi_[k]= array[k];
     597}
     598
     599template <class T>
     600void SphereThetaPhi<T>::setmTheta(double* array, int m)
     601  //remplit  le tableau contenant les valeurs limites de theta
     602  //--
     603{
     604  Theta_= new double[m];
     605  for(int k = 0; k < m; k++) Theta_[k]= array[k];
     606}
     607
     608template <class T>
     609void SphereThetaPhi<T>::setDataBlock(T* data, int m)
     610  // remplit  le vecteur des contenus des pixels
     611{
     612  pixels_.FillFrom(m,data);
     613}
     614
     615template <class T>
     616void SphereThetaPhi<T>::print(ostream& os) const
     617{
     618  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
     619  //
     620  os << " NTheta_= " << NTheta_ << endl;
     621  os << " NPix_    = " << NPix_   << endl;
     622  os << " Omega_  =  " << Omega_   << endl;
     623
     624  os << " contenu de NPhi_ : ";
     625  int i;
     626  for(i=0; i < NTheta_; i++)
     627    {
     628      if(i%5 == 0) os << endl;
     629      os << NPhi_[i] <<", ";
     630    }
     631  os << endl;
     632
     633  os << " contenu de Theta_ : ";
     634  for(i=0; i < NTheta_+1; i++)
     635    {
     636      if(i%5 == 0) os << endl;
     637      os << Theta_[i] <<", ";
     638    }
     639  os << endl;
     640
     641  os << " contenu de TNphi_ : ";
     642  for(i=0; i < NTheta_+1; i++)
     643    {
     644      if(i%5 == 0) os << endl;
     645      os << TNphi_[i] <<", ";
     646    }
     647  os << endl;
     648
     649  os << " contenu de pixels : ";
     650  for(i=0; i < NPix_; i++)
     651    {
     652      if(i%5 == 0) os << endl;
     653      os <<  pixels_(i) <<", ";
     654    }
     655  os << endl;
     656}
     657
     658///////////////////////////////////////////////////////////
     659// --------------------------------------------------------
     660//   Les objets delegues pour la gestion de persistance
     661// --------------------------------------------------------
     662//////////////////////////////////////////////////////////
     663
     664template <class T>
     665FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi()
     666{
     667  dobj= new SphereThetaPhi<T>;
     668  ownobj= true;
     669}
     670
     671template <class T>
     672FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename)
     673{
     674  dobj= new SphereThetaPhi<T>;
     675  ownobj= true;
     676  Read(filename);
     677}
     678
     679template <class T>
     680FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
     681{
     682  dobj= new SphereThetaPhi<T>(obj);
     683  ownobj= true;
     684}
     685
     686template <class T>
     687FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj)
     688{
     689  dobj= obj;
     690  ownobj= false;
     691}
     692
     693template <class T>
     694FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi()
     695{
     696  if (ownobj && dobj) delete dobj;
     697}
     698
     699template <class T>
     700AnyDataObj* FIO_SphereThetaPhi<T>::DataObj()
     701{
     702  return(dobj);
     703}
     704
     705template <class T>
     706void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is)
     707{
     708  cout << " FIO_SphereThetaPhi:: ReadSelf " << endl;
     709
     710  if(dobj == NULL)
     711    {
     712      dobj= new SphereThetaPhi<T>;
     713    }
     714
     715  // Pour savoir s'il y avait un DVList Info associe
     716  char strg[256];
     717  is.GetLine(strg, 255);
     718  bool hadinfo= false;
     719  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
     720  if(hadinfo)
     721    {    // Lecture eventuelle du DVList Info
     722      is >> dobj->Info();
     723    }
     724
     725  int mNTheta;
     726  is.GetI4(mNTheta); 
     727  dobj->setSizeIndex(mNTheta);
     728
     729  int mNPix;
     730  is.GetI4(mNPix);
     731  dobj->setNbPixels(mNPix);
     732
     733  double mOmeg;
     734  is.GetR8(mOmeg);
     735  dobj->setPixSolAngle(mOmeg);
     736
     737  int* mNphi= new int[mNTheta];
     738  is.GetI4s(mNphi, mNTheta);
     739  dobj->setmNPhi(mNphi, mNTheta);
     740  delete [] mNphi;
     741
     742  int* mTNphi= new int[mNTheta+1];
     743  is.GetI4s(mTNphi, mNTheta+1);
     744  dobj->setmTNphi(mTNphi, mNTheta+1);
     745  delete [] mTNphi;
     746
     747  double* mTheta= new double[mNTheta+1];
     748  is.GetR8s(mTheta, mNTheta+1);
     749  dobj->setmTheta(mTheta, mNTheta+1);
     750  delete [] mTheta;
     751
     752  T* mPixels= new T[mNPix];
     753  PIOSReadArray(is, mPixels, mNPix);
     754  dobj->setDataBlock(mPixels, mNPix);
     755  delete [] mPixels;
     756}
     757
     758template <class T>
     759void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const
     760{
     761  cout << " FIO_SphereThetaPhi:: WriteSelf " << endl;
     762
     763  if(dobj == NULL)
     764    {
     765      cout << " WriteSelf:: dobj= null " << endl;
     766      return;
     767    }
     768
     769  char strg[256];
     770  int mNTheta= dobj->SizeIndex();
     771  int mNPix  = dobj->NbPixels();
     772 
     773  if(dobj->ptrInfo())
     774    {
     775      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo",mNTheta,mNPix);
     776      os.PutLine(strg);
     777      os << dobj->Info();
     778    }
     779  else
     780    {
     781      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d ",mNTheta,mNPix);
     782      os.PutLine(strg); 
     783    }
     784
     785  os.PutI4(mNTheta); 
     786  os.PutI4(mNPix);
     787  os.PutR8(dobj->PixSolAngle(0));
     788  os.PutI4s(dobj->getmNPhi() , mNTheta);
     789  os.PutI4s(dobj->getmTNphi(), mNTheta+1);
     790  os.PutR8s(dobj->getmTheta(), mNTheta+1);
     791  //os.Put((dobj->getDataBlock())->Data(), mNPix);
     792  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix);
     793}
    6794
    7795#ifdef __CXX_PRAGMA_TEMPLATES__
     
    25813template class FIO_SphereThetaPhi< complex<double> >;
    26814#endif
    27 
    28 //***************************************************************
    29 //++
    30 // Class        SphereThetaPhi
    31 //
    32 //
    33 // include      spherethetaphi.h nbmath.h
    34 //
    35 //    Découpage de la sphère selon theta et phi, chaque
    36 //    hémisphère étant découpé en (m-1) parallèles (l'équateur compte pour du
    37 //    beurre), chacune des m bandes de theta ainsi définies étant découpée par
    38 //    des méridiens équirepartis, ce découpage étant fait de sorte que tous
    39 //    les pixels aient la même surface et soient le plus carré possible.
    40 //    On commence par découper l'hémisphère de z positif en partant du pôle et
    41 //    en allant vers l'équateur. Le premier pixel est la calotte polaire,
    42 //    il est circulaire et centré sur theta=0.
    43 //--
    44 //++
    45 //
    46 // Links        Parents
    47 //
    48 //    SphericalMap
    49 //--
    50 //++
    51 //
    52 // Links        Descendants
    53 //
    54 //     
    55 //--
    56 
    57 /* --Methode-- */
    58 //++
    59 // Titre        Constructeurs
    60 //--
    61 //++
    62 
    63 template <class T>
    64 SphereThetaPhi<T>::SphereThetaPhi()
    65 
    66 //--
    67 {
    68   InitNul();
    69 }
    70 
    71 
    72 /* --Methode-- */
    73 
    74 //++
    75 template <class T>
    76 SphereThetaPhi<T>::SphereThetaPhi(int m)
    77 
    78 //    Constructeur : m est le nombre de bandes en theta sur un hémisphère
    79 //    (la calotte constituant la premiere bande).
    80 //    pet est le nombre de pixels (pétales) de la bande en contact avec la
    81 //    calotte polaire. Pour l'instant pet est inopérant!
    82 //--
    83 {
    84   InitNul();
    85   Pixelize(m);
    86 }
    87 
    88 template <class T>
    89 SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)
    90 {
    91   if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
    92   NTheta_= s.NTheta_;
    93   NPix_  = s.NPix_;
    94   NPhi_  = new int[NTheta_];
    95   Theta_ = new double[NTheta_+1];
    96   TNphi_ = new int[NTheta_+1];
    97   for(int k = 0; k < NTheta_; k++)
    98     {
    99       NPhi_[k] = s.NPhi_[k];
    100       Theta_[k]= s.Theta_[k];
    101       TNphi_[k]= s.TNphi_[k];
    102     }
    103   Theta_[NTheta_]= s.Theta_[NTheta_];
    104   TNphi_[NTheta_]= s.TNphi_[NTheta_];
    105   Omega_ = s.Omega_;
    106   pixels_= s.pixels_;
    107 }
    108 
    109 //++
    110 // Titre        Destructeur
    111 //--
    112 //++
    113 template <class T>
    114 SphereThetaPhi<T>::~SphereThetaPhi()
    115 
    116 //--
    117 {
    118   Clear();
    119 }
    120 
    121 //++
    122 // Titre        Méthodes
    123 //--
    124 template <class T>
    125 void SphereThetaPhi<T>::InitNul()
    126 //
    127 //    initialise à zéro les variables de classe pointeurs
    128 {
    129   NTheta_= 0;
    130   NPix_  = 0;
    131   Theta_ = NULL;
    132   NPhi_  = NULL;
    133   TNphi_ = NULL;
    134   pixels_.Reset();
    135 }
    136 
    137 /* --Methode-- */
    138 template <class T>
    139 void SphereThetaPhi<T>::Clear()
    140 
    141 {
    142   if(Theta_)  delete[] Theta_;
    143   if(NPhi_ )  delete[] NPhi_;
    144   if(TNphi_)  delete[] TNphi_;
    145   InitNul();
    146 }
    147 
    148 //++
    149 template <class T>
    150 void SphereThetaPhi<T>::Resize(int m)
    151 //   re-pixelize the sphere
    152 //--
    153 {
    154   Clear();
    155   Pixelize(m);
    156 }
    157 
    158 /* --Methode-- */
    159 //++
    160 template <class T>
    161 int SphereThetaPhi<T>::NbPixels() const
    162 
    163 //    Retourne le nombre de pixels du découpage
    164 //--
    165 {
    166   return(NPix_);
    167 }
    168 
    169 /* --Methode-- */
    170 //++
    171 template <class T>
    172 T& SphereThetaPhi<T>::PixVal(int k)
    173 
    174 //    Retourne la valeur du contenu du pixel d'indice k
    175 //--
    176 {
    177   if((k < 0) || (k >= NPix_))
    178     {
    179       //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    180       cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl;
    181       THROW(rangeCheckErr);
    182     }
    183   return pixels_(k);
    184 }
    185 
    186 //++
    187 template <class T>
    188 T const& SphereThetaPhi<T>::PixVal(int k) const
    189 
    190 //    Retourne la valeur du contenu du pixel d'indice k
    191 //--
    192 {
    193   if((k < 0) || (k >= NPix_))
    194     {
    195       cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
    196       //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    197       throw "SphereThetaPhi::PIxVal Pixel index out of range";
    198     }
    199   return *(pixels_.Data()+k);
    200 }
    201 
    202 /* --Methode-- */
    203 //++
    204 template <class T>
    205 int SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const
    206 
    207 //  Retourne l'indice du pixel vers lequel pointe une direction définie par
    208 //    ses coordonnées sphériques
    209 //--
    210 {
    211   double dphi;
    212   int i,j,k;
    213   bool fgzn = false;
    214 
    215   if((theta > Pi) || (theta < 0.)) return(-1);
    216   if((phi < 0.) || (phi > DeuxPi)) return(-1);
    217   if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;}
    218  
    219   // La bande d'indice kt est limitée par les valeurs de theta contenues dans
    220   // Theta_[kt] et Theta_[kt+1]
    221   for( i=1; i< NTheta_; i++ )
    222     if( theta < Theta_[i] ) break;
    223  
    224   dphi= DeuxPi/(double)NPhi_[i-1];
    225  
    226   if (fgzn)  k= NPix_-TNphi_[i]+(int)(phi/dphi);
    227   else k= TNphi_[i-1]+(int)(phi/dphi);
    228   return(k);
    229 }
    230 
    231 /* --Methode-- */
    232 //++
    233 template <class T>
    234 void SphereThetaPhi<T>::PixThetaPhi(int k,double& theta,double& phi) const
    235 
    236 //    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
    237 //--
    238 {
    239   int i;
    240   bool fgzn = false;
    241  
    242   if((k < 0) || (k >= NPix_)) {theta = -99999.; phi = -99999.; return; }
    243   if( k >= NPix_/2)  {fgzn = true; k = NPix_-1-k;}
    244 
    245   // recupère l'indice i de la tranche qui contient le pixel k
    246   for( i=0; i< NTheta_; i++ )
    247     if( k < TNphi_[i+1] ) break;
    248 
    249   // angle theta
    250   theta= 0.5*(Theta_[i]+Theta_[i+1]);
    251   if (fgzn) theta= Pi-theta;
    252  
    253   // angle phi
    254   k -= TNphi_[i];
    255   phi= DeuxPi/(double)NPhi_[i]*(double)(k+.5);
    256   if (fgzn) phi= DeuxPi-phi;
    257 }
    258 
    259 //++
    260 template <class T>
    261 double SphereThetaPhi<T>::PixSolAngle(int dummy) const
    262 
    263 //    Pixel Solid angle  (steradians)
    264 //    All the pixels have the same solid angle. The dummy argument is
    265 //    for compatibility with eventual pixelizations which would not
    266 //    fulfil this requirement.
    267 //--
    268 {
    269   return Omega_;
    270 }
    271 
    272 /* --Methode-- */
    273 //++
    274 template <class T>
    275 void SphereThetaPhi<T>::Limits(int k,double& tetMin,double& tetMax,double& phiMin,double& phiMax)
    276 
    277 //    Retourne les valeurs de theta et phi limitant le pixel d'indice k
    278 //--
    279   {
    280   int j;
    281   double dphi;
    282   bool fgzn= false;
    283  
    284   if((k < 0) || (k >= NPix_)) {
    285     tetMin= -99999.;
    286     phiMin= -99999.; 
    287     tetMax= -99999.;
    288     phiMax= -99999.; 
    289     return;
    290   }
    291  
    292   // si on se trouve dans l'hémisphère Sud
    293   if(k >= NPix_/2) {
    294     fgzn= true;
    295     k= NPix_-1-k;
    296   }
    297  
    298   // on recupere l'indice i de la tranche qui contient le pixel k
    299   int i;
    300   for( i=0; i< NTheta_; i++ )
    301     if(k < TNphi_[i+1]) break;
    302  
    303   // valeurs limites de theta dans l'hemisphere Nord
    304   tetMin= Theta_[i];
    305   tetMax= Theta_[i+1];
    306   // valeurs limites de theta dans l'hemisphere Sud
    307   if (fgzn) {
    308     tetMin= Pi-Theta_[i+1];
    309     tetMax= Pi-Theta_[i];
    310   }
    311  
    312   // pixel correspondant dans l'hemisphere Nord
    313   if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1;
    314  
    315   // indice j de discretisation ( phi= j*dphi )
    316   j= k-TNphi_[i];
    317   dphi= DeuxPi/(double)NPhi_[i];
    318  
    319   // valeurs limites de phi
    320   phiMin= dphi*(double)(j);
    321   phiMax= dphi*(double)(j+1);
    322   return;
    323 }
    324 
    325 /* --Methode-- */
    326 //++
    327 template <class T>
    328 int SphereThetaPhi<T>::NbThetaSlices() const
    329 
    330 //    Retourne le nombre de tranches en theta sur la sphere
    331 //--
    332 {
    333   int nbslices;
    334   nbslices= 2*NTheta_;
    335   return(nbslices);
    336 }
    337 
    338 /* --Methode-- */
    339 //++
    340 template <class T>
    341 int SphereThetaPhi<T>::NPhi(int kt) const
    342 
    343 //    Retourne le nombre de pixels en phi de la tranche kt
    344 //--
    345 {
    346   int nbpix; 
    347   // verification
    348   if((kt < 0) || (kt >= 2*NTheta_)) return(-1);
    349  
    350   // si on se trouve dans l'hemisphere Sud
    351   if(kt >= NTheta_) {
    352     kt= 2*NTheta_-1-kt;
    353   }
    354  
    355   // nombre de pixels
    356   nbpix= NPhi_[kt];
    357   return(nbpix);
    358 }
    359 
    360 
    361 /* --Methode-- */
    362 //++
    363 template <class T>
    364 void SphereThetaPhi<T>::Theta(int kt,double& tetMin,double& tetMax)
    365 
    366 //    Retourne les valeurs de theta limitant la tranche kt
    367 //--
    368 {
    369   bool fgzn= false;
    370   // verification
    371   if( (kt< 0) || (kt>= 2*NTheta_) ) {
    372     tetMin= -99999.;
    373     tetMax= -99999.;
    374     return;
    375   }
    376 
    377   // si on se trouve dans l'hemisphere Sud
    378   if( kt >= NTheta_ ) {
    379     fgzn= true;
    380     kt= 2*NTheta_-1-kt;
    381   }
    382  
    383   // valeurs limites de theta dans l'hemisphere Nord
    384   tetMin= Theta_[kt];
    385   tetMax= Theta_[kt+1];
    386   // valeurs limites de theta dans l'hemisphere Sud
    387   if (fgzn) {
    388     tetMin= Pi-Theta_[kt+1];
    389     tetMax= Pi-Theta_[kt];
    390   } 
    391 }
    392 
    393 /* --Methode-- */
    394 //++
    395 template <class T>
    396 void SphereThetaPhi<T>::Phi(int kt,int jp,double& phiMin,double& phiMax)
    397 
    398 //    Retourne les valeurs de phi limitant le pixel jp de la tranche kt
    399 //--
    400 {
    401   // verification
    402   if((kt < 0) || (kt >= 2*NTheta_)) {
    403     phiMin= -99999.;
    404     phiMax= -99999.;
    405     return;
    406   }
    407  
    408   // si on se trouve dans l'hemisphere Sud
    409   if(kt >= NTheta_) kt= 2*NTheta_-1-kt;
    410  
    411   // verifie si la tranche kt contient au moins jp pixels
    412   if( (jp< 0) || (jp >= NPhi_[kt]) ) {
    413     phiMin= -88888.;
    414     phiMax= -88888.;
    415     return;
    416   }
    417  
    418   double dphi= DeuxPi/(double)NPhi_[kt];
    419   phiMin= dphi*(double)(jp);
    420   phiMax= dphi*(double)(jp+1);
    421   return;
    422 }
    423 
    424 /* --Methode-- */
    425 //++
    426 template <class T>
    427 int SphereThetaPhi<T>::Index(int kt,int jp) const
    428 
    429 //    Retourne l'indice du pixel d'indice jp dans la tranche kt
    430 //--
    431 {
    432   int k;
    433   bool fgzn= false;
    434  
    435   // si on se trouve dans l'hemisphere Sud
    436   if(kt >= NTheta_) {
    437     fgzn= true;
    438     kt= 2*NTheta_-1-kt;
    439   }
    440  
    441   // si la tranche kt contient au moins i pixels
    442   if( (jp>=0) && (jp<NPhi_[kt]) )
    443     {
    444       // dans l'hemisphere Sud
    445       if (fgzn) k= NPix_-TNphi_[kt+1]+jp;
    446       // dans l'hemisphere Nord
    447       else k= TNphi_[kt]+jp;
    448     }
    449   else
    450     {
    451       k= 9999;
    452       printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
    453     }
    454   return(k);
    455 }
    456 
    457 /* --Methode-- */
    458 //++
    459 template <class T>
    460 void SphereThetaPhi<T>::ThetaPhiIndex(int k,int& kt,int& jp)
    461 
    462 //    Retourne les indices kt et jp du pixel d'indice k
    463 //--
    464 {
    465   bool fgzn= false; 
    466   // si on se trouve dans l'hemisphere Sud
    467   if(k >= NPix_/2)
    468     {
    469       fgzn= true;
    470       k= NPix_-1-k;
    471     }
    472  
    473   // on recupere l'indice kt de la tranche qui contient le pixel k
    474   int i;
    475   for(i = 0; i < NTheta_; i++)
    476     if(k < TNphi_[i+1]) break;
    477  
    478   // indice  kt de tranche
    479   if (fgzn) kt= 2*NTheta_-1-i;
    480   else kt= i;
    481  
    482   // indice jp de pixel
    483   if (fgzn) jp= TNphi_[i+1]-k-1;
    484   else jp= k-TNphi_[i]; 
    485 }
    486 //++
    487 template <class T>
    488 void  SphereThetaPhi<T>::Pixelize(int m)
    489 
    490 //    effectue le découpage en pixels (m et pet ont la même signification
    491 //    que pour le constructeur)
    492 //
    493 //    Chaque bande de theta sera découpée en partant de phi=0 ...
    494 //    L'autre hémisphère est parcourue dans le même sens en phi et de
    495 //    l'équateur vers le pôle (le pixel qui suit le dernier de la bande la plus
    496 //    proche de l'équateur a z>0 est celui de plus petit phi de la bande la
    497 //    plus proche de l'equateur a z<0).
    498 //--
    499 {
    500   int ntotpix,i,j;
    501  
    502   // Decodage et controle des arguments d'appel :
    503   // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
    504   if (m < 2) m = 2;
    505   if (m > 16384) m = 16384;
    506  
    507   // On memorise les arguments d'appel
    508   NTheta_ = m; 
    509  
    510   // On commence par decouper l'hemisphere z>0.
    511   // Creation des vecteurs contenant :
    512   // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
    513   Theta_= new double[m+1];
    514  
    515   // Le nombre de pixels en phi de chacune des bandes en theta
    516   NPhi_ = new int[m];
    517  
    518   // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
    519   // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
    520   TNphi_= new int[m+1];
    521  
    522   // Calcul du nombre total de pixels dans chaque bande pour optimiser
    523   // le rapport largeur/hauteur des pixels
    524  
    525   //calotte polaire
    526   TNphi_[0]= 0;
    527   NPhi_[0] = 1;
    528  
    529   //bandes jusqu'a l'equateur
    530   for(j = 1; j < m; j++)
    531     {
    532       TNphi_[j]= TNphi_[j-1]+NPhi_[j-1];
    533       NPhi_[j] = (int)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
    534     }
    535  
    536   // Nombre total de pixels sur l'hemisphere
    537   ntotpix  = TNphi_[m-1]+NPhi_[m-1];
    538   TNphi_[m]= ntotpix;
    539   // et sur la sphere entiere
    540   NPix_= 2*ntotpix;
    541  
    542   // Creation et initialisation du vecteur des contenus des pixels
    543   pixels_.ReSize(NPix_);
    544   pixels_.Reset();
    545 
    546   // Determination des limites des bandes en theta :
    547   // omeg est l'angle solide couvert par chaque pixel,
    548   // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg
    549   // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide
    550   //de la
    551   // calotte allant du pole a la limite haute de la bande kt vaut
    552   // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...
    553  
    554   double omeg2pi= 1./(double)ntotpix;
    555   Omega_ = omeg2pi*DeuxPi;
    556  
    557   for(j=0; j <= m; j++)
    558     {
    559       Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi);
    560     }
    561 }
    562 
    563 //++
    564 template <class T>
    565 void SphereThetaPhi<T>::GetThetaSlice(int index,double& theta, TVector<double>& phi, TVector<T>& value) const
    566 
    567 //    Retourne, pour la tranche en theta d'indice 'index' le theta
    568 //    correspondant, un vecteur (Peida) contenant les phi des pixels de
    569 //    la tranche, un vecteur (Peida) contenant les valeurs de pixel
    570 //    correspondantes
    571 //--
    572 
    573 {
    574   cout << "entree GetThetaSlice, couche no " << index << endl;
    575 
    576   if(index < 0 || index > NbThetaSlices())
    577     {
    578       // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    579       cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
    580       THROW(rangeCheckErr);
    581     }
    582 
    583   int iring= Index(index,0);
    584   int bid  = this->NPhi(index);
    585   int lring  = bid;
    586   cout << " iring= " << iring << " lring= " << lring << endl;
    587 
    588   phi.ReSize(lring);
    589   value.ReSize(lring);
    590   double Te= 0.;
    591   double Fi= 0.;
    592   for(int kk = 0; kk < lring; kk++)
    593     {
    594       PixThetaPhi(kk+iring,Te,Fi);
    595       phi(kk)= Fi;
    596       value(kk)= PixVal(kk+iring);
    597     }
    598   theta= Te;
    599 }
    600 
    601 template <class T>
    602 void SphereThetaPhi<T>::setmNPhi(int* array, int m)
    603   //remplit le tableau contenant le nombre de pixels en phi de chacune des bandes en theta
    604   //--
    605 {
    606   NPhi_= new int[m];
    607   for(int k = 0; k < m; k++) NPhi_[k]= array[k];
    608 }
    609 
    610 template <class T>
    611 void SphereThetaPhi<T>::setmTNphi(int* array, int m)
    612   //remplit  le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes
    613   //--
    614 {
    615   TNphi_= new int[m];
    616   for(int k = 0; k < m; k++) TNphi_[k]= array[k];
    617 }
    618 
    619 template <class T>
    620 void SphereThetaPhi<T>::setmTheta(double* array, int m)
    621   //remplit  le tableau contenant les valeurs limites de theta
    622   //--
    623 {
    624   Theta_= new double[m];
    625   for(int k = 0; k < m; k++) Theta_[k]= array[k];
    626 }
    627 
    628 template <class T>
    629 void SphereThetaPhi<T>::setDataBlock(T* data, int m)
    630   // remplit  le vecteur des contenus des pixels
    631 {
    632   pixels_.FillFrom(m,data);
    633 }
    634 
    635 template <class T>
    636 void SphereThetaPhi<T>::print(ostream& os) const
    637 {
    638   if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
    639   //
    640   os << " NTheta_= " << NTheta_ << endl;
    641   os << " NPix_    = " << NPix_   << endl;
    642   os << " Omega_  =  " << Omega_   << endl;
    643 
    644   os << " contenu de NPhi_ : ";
    645   for(int i=0; i < NTheta_; i++)
    646     {
    647       if(i%5 == 0) os << endl;
    648       os << NPhi_[i] <<", ";
    649     }
    650   os << endl;
    651 
    652   os << " contenu de Theta_ : ";
    653   for(int i=0; i < NTheta_+1; i++)
    654     {
    655       if(i%5 == 0) os << endl;
    656       os << Theta_[i] <<", ";
    657     }
    658   os << endl;
    659 
    660   os << " contenu de TNphi_ : ";
    661   for(int i=0; i < NTheta_+1; i++)
    662     {
    663       if(i%5 == 0) os << endl;
    664       os << TNphi_[i] <<", ";
    665     }
    666   os << endl;
    667 
    668   os << " contenu de pixels : ";
    669   for(int i=0; i < NPix_; i++)
    670     {
    671       if(i%5 == 0) os << endl;
    672       os <<  pixels_(i) <<", ";
    673     }
    674   os << endl;
    675 }
    676 
    677 ///////////////////////////////////////////////////////////
    678 // --------------------------------------------------------
    679 //   Les objets delegues pour la gestion de persistance
    680 // --------------------------------------------------------
    681 //////////////////////////////////////////////////////////
    682 
    683 template <class T>
    684 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi()
    685 {
    686   dobj= new SphereThetaPhi<T>;
    687   ownobj= true;
    688 }
    689 
    690 template <class T>
    691 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename)
    692 {
    693   dobj= new SphereThetaPhi<T>;
    694   ownobj= true;
    695   Read(filename);
    696 }
    697 
    698 template <class T>
    699 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
    700 {
    701   dobj= new SphereThetaPhi<T>(obj);
    702   ownobj= true;
    703 }
    704 
    705 template <class T>
    706 FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj)
    707 {
    708   dobj= obj;
    709   ownobj= false;
    710 }
    711 
    712 template <class T>
    713 FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi()
    714 {
    715   if (ownobj && dobj) delete dobj;
    716 }
    717 
    718 template <class T>
    719 AnyDataObj* FIO_SphereThetaPhi<T>::DataObj()
    720 {
    721   return(dobj);
    722 }
    723 
    724 template <class T>
    725 void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is)
    726 {
    727   cout << " FIO_SphereThetaPhi:: ReadSelf " << endl;
    728 
    729   if(dobj == NULL)
    730     {
    731       dobj= new SphereThetaPhi<T>;
    732     }
    733 
    734   // Pour savoir s'il y avait un DVList Info associe
    735   char strg[256];
    736   is.GetLine(strg, 255);
    737   bool hadinfo= false;
    738   if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
    739   if(hadinfo)
    740     {    // Lecture eventuelle du DVList Info
    741       is >> dobj->Info();
    742     }
    743 
    744   int mNTheta;
    745   is.GetI4(mNTheta); 
    746   dobj->setSizeIndex(mNTheta);
    747 
    748   int mNPix;
    749   is.GetI4(mNPix);
    750   dobj->setNbPixels(mNPix);
    751 
    752   double mOmeg;
    753   is.GetR8(mOmeg);
    754   dobj->setPixSolAngle(mOmeg);
    755 
    756   int* mNphi= new int[mNTheta];
    757   is.GetI4s(mNphi, mNTheta);
    758   dobj->setmNPhi(mNphi, mNTheta);
    759   delete [] mNphi;
    760 
    761   int* mTNphi= new int[mNTheta+1];
    762   is.GetI4s(mTNphi, mNTheta+1);
    763   dobj->setmTNphi(mTNphi, mNTheta+1);
    764   delete [] mTNphi;
    765 
    766   double* mTheta= new double[mNTheta+1];
    767   is.GetR8s(mTheta, mNTheta+1);
    768   dobj->setmTheta(mTheta, mNTheta+1);
    769   delete [] mTheta;
    770 
    771   T* mPixels= new T[mNPix];
    772   PIOSReadArray(is, mPixels, mNPix);
    773   dobj->setDataBlock(mPixels, mNPix);
    774   delete [] mPixels;
    775 }
    776 
    777 template <class T>
    778 void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const
    779 {
    780   cout << " FIO_SphereThetaPhi:: WriteSelf " << endl;
    781 
    782   if(dobj == NULL)
    783     {
    784       cout << " WriteSelf:: dobj= null " << endl;
    785       return;
    786     }
    787 
    788   char strg[256];
    789   int mNTheta= dobj->SizeIndex();
    790   int mNPix  = dobj->NbPixels();
    791  
    792   if(dobj->ptrInfo())
    793     {
    794       sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo",mNTheta,mNPix);
    795       os.PutLine(strg);
    796       os << dobj->Info();
    797     }
    798   else
    799     {
    800       sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d ",mNTheta,mNPix);
    801       os.PutLine(strg); 
    802     }
    803 
    804   os.PutI4(mNTheta); 
    805   os.PutI4(mNPix);
    806   os.PutR8(dobj->PixSolAngle(0));
    807   os.PutI4s(dobj->getmNPhi() , mNTheta);
    808   os.PutI4s(dobj->getmTNphi(), mNTheta+1);
    809   os.PutR8s(dobj->getmTheta(), mNTheta+1);
    810   //os.Put((dobj->getDataBlock())->Data(), mNPix);
    811   PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix);
    812 }
Note: See TracChangeset for help on using the changeset viewer.