Changeset 470 in Sophya for trunk/SophyaLib/Samba/spheregorski.cc


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/spheregorski.cc

    r276 r470  
    22#include "spheregorski.h"
    33#include "strutil.h"
    4 extern "C" {
     4#include <complex>
     5#include "piocmplx.h"
     6
     7extern "C"
     8{
    59#include <stdio.h>
    610#include <stdlib.h>
    … …  
    812}
    913     
    10 extern "C" {
    11   //void ang2pix_ring_(int&,double&,double&,int&);
    12   //void pix2ang_ring_(int&, int&, double&, double&);
    13 //void ang2pix_nest_(int&,double&,double&,int&);
    14   //void pix2ang_nest_(int&, int&, double&, double&);
    15 //void nest2ring_(int&, int&, int&);
    16 //void ring2nest_(int&, int&, int&);
    17 void anafast_(int&, int&, int&,double&,float*,float*,float*,float*,
    18               float*,float*,float*);
    19 void synfast_(int&, int&, int&,int&, float&,float*,float*,float*,
    20               double*, double*,double*,double*,double*,float*);
    21 void  input_map_(char*,float*,int&);
    22 void  ecrire_fits_(int&,int&,int&,float*,char*,char*,int&,float&,float&);
    23 
    24 }
     14extern "C"
     15{
     16void anafast_(int&, int&, int&,double&,float*,float*,float*,float*,float*,float*,float*);
     17void synfast_(int&, int&, int&,int&, float&,float*,float*,float*,double*, double*,double*,double*,double*,float*);
     18}
     19
     20//*******************************************************************
     21// Class PIXELS_XY
     22//  Construction des tableaux necessaires a la traduction des indices RING en
     23//  indices NESTED (ou l'inverse)
     24//*******************************************************************
     25
     26PIXELS_XY::PIXELS_XY()
     27{
     28  cout << " appel du constructeur PIXELS_XY " <<endl;
     29  pix2x_.ReSize(1024);
     30  pix2x_.Reset();
     31  pix2y_.ReSize(1024);
     32  pix2y_.Reset();
     33  x2pix_.ReSize(128);
     34  x2pix_.Reset();
     35  y2pix_.ReSize(128);
     36  y2pix_.Reset();
     37  mk_pix2xy();
     38  mk_xy2pix();
     39}
     40
     41PIXELS_XY& PIXELS_XY::instance()
     42{
     43  static PIXELS_XY single;
     44  return (single);
     45}
     46
     47void PIXELS_XY::mk_pix2xy()
     48{
     49  /*
     50    ==================================================
     51    subroutine mk_pix2xy
     52    ==================================================
     53    c     constructs the array giving x and y in the face from pixel number
     54    c     for the nested (quad-cube like) ordering of pixels
     55    c
     56    c     the bits corresponding to x and y are interleaved in the pixel number
     57    c     one breaks up the pixel number by even and odd bits
     58    ==================================================
     59    */
     60  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
     61  //  (16/12/98)
     62
     63  int kpix, jpix, IX, IY, IP, ID;
     64 
     65  for(kpix = 0; kpix < 1024; kpix++)
     66    {
     67      jpix = kpix;
     68      IX = 0;
     69      IY = 0;
     70      IP = 1 ;//  ! bit position (in x and y)
     71      while( jpix!=0 )
     72        { // ! go through all the bits
     73          ID=jpix%2;//  ! bit value (in kpix), goes in ix
     74          jpix = jpix/2;
     75          IX = ID*IP+IX;
     76     
     77          ID=jpix%2;//  ! bit value (in kpix), goes in iy
     78          jpix = jpix/2;
     79          IY = ID*IP+IY;
     80     
     81          IP = 2*IP;//        ! next bit (in x and y)
     82        }
     83      pix2x_(kpix) = IX;//     ! in 0,31
     84      pix2y_(kpix) = IY;//     ! in 0,31
     85    }
     86}
     87
     88void PIXELS_XY::mk_xy2pix()
     89{
     90  /*
     91    =================================================
     92    subroutine mk_xy2pix
     93    =================================================
     94    c     sets the array giving the number of the pixel lying in (x,y)
     95    c     x and y are in {1,128}
     96    c     the pixel number is in {0,128**2-1}
     97    c
     98    c     if  i-1 = sum_p=0  b_p * 2^p
     99    c     then ix = sum_p=0  b_p * 4^p
     100    c          iy = 2*ix
     101    c     ix + iy in {0, 128**2 -1}
     102    =================================================
     103    */
     104  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
     105  //  (16/12/98)
     106
     107  int K,IP,I,J,ID;
     108  for(I = 1; I <= 128; I++)
     109    {
     110      J  = I-1;//            !pixel numbers
     111      K  = 0;//
     112      IP = 1;//
     113      truc : if( J==0 )
     114        {
     115          x2pix_(I-1) = K;
     116          y2pix_(I-1) = 2*K;
     117        }
     118      else
     119        {
     120          ID = (int)fmod(J,2);
     121          J  = J/2;
     122          K  = IP*ID+K;
     123          IP = IP*4;
     124          goto truc;
     125        }
     126    }     
     127}
     128
    25129//*******************************************************************
    26130//++
    … …  
    81185//++
    82186
    83 SphereGorski::SphereGorski()
    84 
    85 //--
    86 {
     187template<class T>
     188SphereGorski<T>::SphereGorski()
     189
     190//--
     191{
     192  cout<<" appel du constructeur SphereGorski ()" <<endl;
    87193  InitNul();
    88194}
    89195
    90 /* --Methode-- */
    91 //++
    92 SphereGorski::SphereGorski(char* flnm)
    93 
    94 //    Constructeur : charge une image à partir d'un fichier
    95 //--
    96 {
    97   InitNul();
    98   PInPersist s(flnm);
    99   Read(s);
    100 }
    101 
    102 //++
    103 SphereGorski::SphereGorski(int_4 m)
     196//++
     197template<class T>
     198SphereGorski<T>::SphereGorski(int_4 m)
    104199
    105200//    Constructeur : m est la variable nside de l'algorithme de Gorski
    … …  
    108203//--
    109204{
    110   //printf(" initialisation par defaut SphereGorski \n");
     205  cout<<" appel du constructeur SphereGorski (m)" <<endl;
     206
     207  if(m <= 0 || m > 8192)
     208    {
     209      cout << "SphereGorski : m hors bornes [0,8192], m= " << m << endl;
     210      exit(1);
     211    }
     212  // verifier que m est une puissance de deux
     213  int_4 x=m;
     214  while(x%2 == 0) x/=2;
     215  if(x != 1)
     216    { 
     217      cout<<"SphereGorski: m doit etre une puissance de deux, m= "<<m<<endl;
     218      exit(1);
     219    }
     220  InitNul();
     221  Pixelize(m);
     222}
     223
     224template<class T>
     225SphereGorski<T>::SphereGorski(const SphereGorski<T>& s)
     226{
     227  cout << " constructeur de recopie " << endl;
     228  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
     229
     230  nSide_= s.nSide_;
     231  nPix_ = s.nPix_;
     232  omeg_ = s.omeg_;
     233
     234  pixels_= s.pixels_;
     235
     236  nlmax_= s.nlmax_;
     237  nmmax_= s.nmmax_;
     238  iseed_= s.iseed_;
     239  fwhm_ = s.fwhm_;
     240  quadrupole_ = s.quadrupole_;
     241  sym_cut_deg_= s.sym_cut_deg_;
     242  strcpy(powFile_,s.powFile_);
     243}
     244
     245//++
     246// Titre        Destructeur
     247//--
     248//++
     249template<class T>
     250SphereGorski<T>::~SphereGorski()
     251
     252//--
     253{
     254  InitNul();
     255}
     256
     257//++
     258// Titre        Méthodes
     259//--
     260
     261//++
     262template<class T>
     263void SphereGorski<T>::Resize(int_4 m)
     264
     265//    m est la variable nside de l'algorithme de Gorski
     266//    le nombre total de pixels sera Npix =  12*nside**2
     267//    m DOIT OBLIGATOIREMENT ETRE UNE PUISSANCE DE 2 (<= 8192)
     268//--
     269{
    111270  if (m<=0 || m> 8192) {
    112271    cout << "SphereGorski : m hors bornes [0,8192], m= " << m << endl;
    113272    exit(1);
    114273  }
    115 // verifier que m est une puissance de deux
     274  // verifier que m est une puissance de deux
    116275  int_4 x=m;
    117276  while (x%2==0) x/=2;
    118   if (x!=1) {
    119     cout << "SphereGorski : m doit etre une puissance de deux, m= " << m << endl;
    120     exit(1);
    121   }
    122 
     277  if(x != 1)
     278    { 
     279      cout<<"SphereGorski: m doit etre une puissance de deux, m= "<<m<<endl;
     280      exit(1);
     281    }
    123282  InitNul();
    124283  Pixelize(m);
    125   if (pix2x_==NULL) {
    126     pix2x_=new int[1024];
    127     pix2y_=new int[1024];
    128     mk_pix2xy(pix2x_,pix2y_);
    129   }
    130   if (x2pix_==NULL) {
    131     x2pix_=new int[128];
    132     y2pix_=new int[128];
    133     mk_xy2pix(x2pix_,y2pix_);
    134   }
    135 }
    136 
    137 
    138 
    139 
    140 
    141 
    142 //++
    143 // Titre        Destructeur
    144 //--
    145 //++
    146 SphereGorski::~SphereGorski()
    147 
    148 //--
    149 {
    150   Clear();
    151 }
    152 
    153 //++
    154 // Titre        Méthodes
    155 //--
    156 
    157 
    158 
    159 
    160 void  SphereGorski::Pixelize( int_4 m)
     284}
     285
     286template<class T>
     287void  SphereGorski<T>::Pixelize( int_4 m)
    161288
    162289//    prépare la pixelisation Gorski (m a la même signification
    … …  
    166293//--
    167294{
    168  
    169295  // On memorise les arguments d'appel
    170296  nSide_ = m; 
    171  
    172  
     297
    173298  // Nombre total de pixels sur la sphere entiere
    174299  nPix_=12*nSide_*nSide_;
    175300
    176 // pour le moment les tableaux qui suivent seront ranges dans l'ordre
    177 // de l'indexation GORSKY "RING"
    178 // on pourra ulterieurement changer de strategie et tirer profit
    179 // de la dualite d'indexation GORSKY (RING er NEST) : tout dependra
    180 // de pourquoi c'est faire
     301  // pour le moment les tableaux qui suivent seront ranges dans l'ordre
     302  // de l'indexation GORSKY "RING"
     303  // on pourra ulterieurement changer de strategie et tirer profit
     304  // de la dualite d'indexation GORSKY (RING et NEST) : tout dependra
     305  // de pourquoi c'est faire
    181306
    182307  // Creation et initialisation du vecteur des contenus des pixels
    183   mPix_ = new r_8[nPix_];
    184   for(int i=0; i<nPix_; i++)  mPix_[i] = 0.;
     308  pixels_.ReSize(nPix_);
     309  pixels_.Reset();
    185310
    186311  // solid angle per pixel   
    187   omeg_=4*Pi/nPix_;
    188 }
    189 
    190 void SphereGorski::InitNul()
     312  omeg_= 4*Pi/nPix_;
     313}
     314
     315template<class T>
     316void SphereGorski<T>::InitNul()
    191317//
    192318//    initialise à zéro les variables de classe
    193319{
    194   nlmax_=0;
    195   nmmax_=0;
    196   iseed_=0;
    197   fwhm_=0.;
    198   quadrupole_=0.;
    199   sym_cut_deg_=0.;
    200   for (int k=0; k<128;k++) powFile_[k]=' ';
    201   nSide_=0;
    202   nPix_ =0;
    203   mPix_ = NULL;
    204   pix2x_=NULL;
    205   pix2y_=NULL;
    206   x2pix_=NULL;
    207   y2pix_=NULL;
    208   pix2xy_=NULL;
    209 }
     320  nSide_= 0;
     321  nPix_ = 0;
     322  omeg_ = 0.;
     323  pixels_.Reset();
     324
     325  nlmax_= 0;
     326  nmmax_= 0;
     327  iseed_= 0;
     328  fwhm_ = 0.;
     329  quadrupole_ = 0.;
     330  sym_cut_deg_= 0.;
     331  for(int k = 0; k < 128; k++) powFile_[k]=' ';
     332}
     333
    210334/* --Methode-- */
    211 void SphereGorski::Clear()
    212 
    213 {
    214   if (mPix_)   delete[] mPix_;
    215   if (pix2x_)  delete[] pix2x_;
    216   if (pix2y_)  delete[] pix2y_;
    217   if (x2pix_)  delete[] x2pix_;
    218   if (y2pix_)  delete[] y2pix_;
    219   if (pix2xy_) delete pix2xy_;
    220 }
    221 //++
    222 void SphereGorski::WriteSelf(POutPersist& s) const
    223 
    224 //    créer un fichier image
    225 //--
    226 {
    227   char strg[256];
    228   if (mInfo_) sprintf(strg, "SphereGorski: NSlices=%6d  NPix=%9d HasInfo", 
    229                    (int_4)nSide_, (int_4)nPix_);
    230   else
    231     sprintf(strg, "SphereGorski: nSide=%6d  nPix=%9d ",
    232           (int_4)nSide_, (int_4)nPix_);
    233   s.PutLine(strg);
    234   if (mInfo_)  s << (*mInfo_);
    235   s.PutI4(nlmax_);
    236   s.PutI4(nmmax_);
    237   s.PutI4(iseed_);
    238   s.PutI4(nSide_);
    239   s.PutI4(nPix_);
    240   s.PutR4(fwhm_);
    241   s.PutR4(quadrupole_);
    242   s.PutR4(sym_cut_deg_);
    243   s.PutR8(omeg_);
    244   s.PutR8s(mPix_, nPix_);
    245   s.PutLine(powFile_);
    246   return;
    247 }
    248 
    249 /* --Methode-- */
    250 //++
    251 void SphereGorski::ReadSelf(PInPersist& s)
    252 
    253 //    relit un fichier d'image
    254 //--
    255 {
    256   Clear();
    257   char strg[256];
    258   s.GetLine(strg, 255);
    259 // Pour savoir s'il y avait un DVList Info associe
    260   bool hadinfo = false;
    261   if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
    262   if (hadinfo) {    // Lecture eventuelle du DVList Info
    263     if (mInfo_ == NULL)  mInfo_ = new DVList;
    264     s >> (*mInfo_);
    265   }
    266   s.GetI4(nlmax_);
    267   s.GetI4(nmmax_);
    268   s.GetI4(iseed_);
    269   s.GetI4(nSide_);
    270   s.GetI4(nPix_);
    271   s.GetR4(fwhm_);
    272   s.GetR4(quadrupole_);
    273   s.GetR4(sym_cut_deg_);
    274   mPix_=new r_8[nPix_];
    275   s.GetR8(omeg_);
    276   s.GetR8s(mPix_, nPix_);
    277   s.GetLine(powFile_, 127);
    278 
    279   return;
    280 }
    281 
    282 
    283 /* --Methode-- */
    284 //++
    285 void SphereGorski::ReadFits(char flnm[])
    286 
    287 //    remplit la sphere a partir d'un fichier FITS
    288 //--
    289 {
    290   strip(flnm,'B',' ');
    291   if (access(flnm,F_OK) != 0) {perror(flnm); exit(1);}
    292   if(!nPix_) {
    293     cout << " ReadFits : SphereGorski non pixelisee " << endl;
    294     exit(1);
    295   }
    296   int npixtot=nPix_;
    297 
    298 // quand map et mPix_ auront le meme type, map ne sera plus necessaire
    299 
    300   float* map=new float[12*nSide_*nSide_];
    301   input_map_(flnm,map,npixtot);
    302   // Remplissage de la sphère
    303 
    304   for(int  j=0; j<nPix_; j++ ) mPix_[j]=(r_8)map[j];
    305   delete [] map;
    306  
    307 
    308 }
    309 
    310 /* --Methode-- */
    311 //++
    312 void SphereGorski::WriteFits(char flnm[])
    313 
    314 //    ecrit la sphere sur un fichier FITS
    315 
    316 //--
    317 {
    318   strip(flnm,'B',' ');
    319   if (access(flnm,F_OK) == 0) {
    320     cout << " SphereGorski::WriteFits : le fichier existe deja" << endl;
    321     exit(1);
    322   }
    323   //else cout << "un fichier sera cree " << endl;
    324   if(!nPix_) {
    325     cout << " WriteFits : SphereGorski non pixelisee " << endl;
    326     exit(1);
    327   }
    328 
    329 
    330   char infile[128];
    331   for (int k=0; k< 128; k++) infile[k]=powFile_[k];
    332 
    333 // quand map et mPix_ auront le meme type, map ne sera plus necessaire
    334 
    335   float* map=new float[12*nSide_*nSide_];
    336 
    337   for(int  j=0; j<nPix_; j++ ) map[j]=(float)mPix_[j];
    338   int nlmax=nlmax_;
    339   int nsmax=nSide_;
    340   int nmmax=nmmax_;
    341   int iseed=iseed_;
    342   float fwhm=fwhm_;
    343   float quadrupole=quadrupole_;
    344   ecrire_fits_(nlmax,nsmax,nmmax,map,infile,flnm,iseed,fwhm,quadrupole);
    345   delete [] map;
    346  
    347 }
    348 
    349 /* --Methode-- */
    350 //++
    351 int_4 SphereGorski::NbPixels() const
     335//++
     336template<class T>
     337int_4 SphereGorski<T>::NbPixels() const
    352338
    353339//    Retourne le nombre de pixels du découpage
    354340//--
    355 {
     341{ 
    356342  return(nPix_);
    357343}
    358344
    359345//++
    360 int_4 SphereGorski::NbThetaSlices() const
     346template<class T>
     347int_4 SphereGorski<T>::NbThetaSlices() const
    361348
    362349//    Retourne le nombre de tranches en theta sur la sphere
    363350//--
    364 {return int_4(4*nSide_-1);}
    365 
    366 
    367 //++
    368 void  SphereGorski::GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const
     351{
     352  return int_4(4*nSide_-1);
     353}
     354
     355//++
     356template<class T>
     357void  SphereGorski<T>::GetThetaSlice(int_4 index, r_4& theta, TVector<float>& phi, TVector<T>& value) const
    369358
    370359//    Retourne, pour la tranche en theta d'indice 'index' le theta
    … …  
    374363//--
    375364{
    376 
    377 cout << "entree GetThetaSlice, couche no " << index << endl;
    378   if (index<0 || index > NbThetaSlices()) {
    379     // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range"));
    380     cout << " SphereGorski::GetThetaSlice : exceptions  a mettre en place" <<endl;
    381   THROW(rangeCheckErr);
    382   }
    383   int_4 iring=0;
    384   int lring=0;
    385   if (index<nSide_-1) {
    386     iring=2*index*(index+1);
    387     lring=4*(index+1);
    388   }else
    389     if (index<3*nSide_) {
    390       iring=2*nSide_*(2*index-nSide_+1);
    391       lring=4*nSide_;
    392     }else
    393       {
    394         int nc=4*nSide_-1-index;
    395         iring=nPix_-2*nc*(nc+1);
    396         lring=4*nc;
    397       }
    398   phi.Realloc(lring);
    399   value.Realloc(lring);
    400   float T=0.;
    401   float F=0.;
    402   for (int kk=0; kk<lring;kk++) {
    403     PixThetaPhi(kk+iring,T,F);
    404     phi(kk)=F;
    405     value(kk)=PixVal(kk+iring);
    406   }
    407   theta=T;
    408 
    409 }
    410 
    411 
     365  cout << "entree GetThetaSlice, couche no " << index << endl;
     366
     367  if (index<0 || index > NbThetaSlices())
     368    {
     369      // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range"));
     370      cout << " SphereGorski::GetThetaSlice : exceptions  a mettre en place" <<endl;
     371      THROW(rangeCheckErr);
     372    }
     373
     374  int_4 iring= 0;
     375  int lring  = 0;
     376  if(index < nSide_-1)
     377    {
     378      iring= 2*index*(index+1);
     379      lring= 4*(index+1);
     380    }
     381  else if(index < 3*nSide_)
     382    {
     383      iring= 2*nSide_*(2*index-nSide_+1);
     384      lring= 4*nSide_;
     385    }
     386  else
     387    {
     388      int nc= 4*nSide_-1-index;
     389      iring = nPix_-2*nc*(nc+1);
     390      lring = 4*nc;
     391    }
     392
     393  phi.ReSize(lring);
     394  value.ReSize(lring);
     395  float TH=0.;
     396  float F =0.;
     397  for(int kk = 0; kk < lring;kk++)
     398    {
     399      PixThetaPhi(kk+iring,TH,F);
     400      phi(kk)= F;
     401      value(kk)= PixVal(kk+iring);
     402    }
     403  theta= TH;
     404}
    412405
    413406/* --Methode-- */
    414407//++
    415 r_8& SphereGorski::PixVal(int_4 k)
     408template<class T>
     409T& SphereGorski<T>::PixVal(int_4 k)
    416410
    417411//    Retourne la valeur du contenu du pixel d'indice "RING" k
    418412//--
    419413{
    420   if ( (k<0) || (k >= nPix_) ) {
    421     // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range"));
    422     cout << " SphereGorski::PIxVal : exceptions  a mettre en place" <<endl;
    423   THROW(rangeCheckErr);
    424   }
    425   return(mPix_[k]);
     414  if((k < 0) || (k >= nPix_))
     415    {
     416      // THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range"));
     417      cout << " SphereGorski::PIxVal : exceptions  a mettre en place" <<endl;
     418      THROW(rangeCheckErr);
     419    }
     420  return pixels_(k);
    426421}
    427422
    428423/* --Methode-- */
    429424//++
    430 r_8 const& SphereGorski::PixVal(int_4 k) const
     425template<class T>
     426T const& SphereGorski<T>::PixVal(int_4 k) const
    431427
    432428//    Retourne la valeur du contenu du pixel d'indice "RING" k
    433429//--
    434430{
    435   if ( (k<0) || (k >= nPix_) ) {
    436     //THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range"));
    437     cout << " SphereGorski::PIxVal : exceptions  a mettre en place" <<endl;
    438   THROW(rangeCheckErr);
    439   }
    440   return(mPix_[k]);
    441 }
    442 
    443 
    444 //++
    445 r_8& SphereGorski::PixValNest(int_4 k)
     431  if((k < 0) || (k >= nPix_))
     432    {
     433      //THROW(out_of_range("SphereGorski::PIxVal Pixel index out of range"));
     434      cout << " SphereGorski::PIxVal : exceptions  a mettre en place" <<endl;
     435      THROW(rangeCheckErr);
     436    }
     437  return *(pixels_.Data()+k);
     438}
     439
     440//++
     441template<class T>
     442T& SphereGorski<T>::PixValNest(int_4 k)
    446443
    447444//    Retourne la valeur du contenu du pixel d'indice "NESTED" k
    448445//--
    449446{
    450   if ( (k<0) || (k >= nPix_) ) {
    451    //THROW(out_of_range("SphereGorski::PIxValNest Pixel index out of range"));
    452     cout << " SphereGorski::PIxValNest : exceptions  a mettre en place" <<endl;
    453   THROW(rangeCheckErr);
    454   }
    455   return mPix_[nest2ring(nSide_,k)];
    456 }
    457 //++
    458 
    459 r_8 const& SphereGorski::PixValNest(int_4 k) const
     447  if((k < 0) || (k >= nPix_))
     448    {
     449      //THROW(out_of_range("SphereGorski::PIxValNest Pixel index out of range"));
     450      cout<<" SphereGorski::PIxValNest : exceptions a mettre en place" <<endl;
     451      THROW(rangeCheckErr);
     452    }
     453  return pixels_(nest2ring(nSide_,k));
     454}
     455//++
     456
     457template<class T>
     458T const& SphereGorski<T>::PixValNest(int_4 k) const
    460459
    461460//    Retourne la valeur du contenu du pixel d'indice "NESTED" k
    462461//--
    463462{
    464   if ( (k<0) || (k >= nPix_) ) {
    465    //THROW(out_of_range("SphereGorski::PIxValNest Pixel index out of range"));
    466     cout << " SphereGorski::PIxValNest : exceptions  a mettre en place" <<endl;
    467   THROW(rangeCheckErr);
    468   }
    469   int_4 pix=nest2ring(nSide_,k);
    470   return mPix_[pix];
    471 }
    472 
    473 
     463  if((k < 0) || (k >= nPix_))
     464    {
     465      //THROW(out_of_range("SphereGorski::PIxValNest Pixel index out of range"));
     466      cout<<" SphereGorski::PIxValNest : exceptions a mettre en place" <<endl;
     467      THROW(rangeCheckErr);
     468  }
     469  int_4 pix= nest2ring(nSide_,k);
     470  return *(pixels_.Data()+pix);
     471}
    474472
    475473/* --Methode-- */
    476474//++
    477 int_4 SphereGorski::PixIndexSph(r_4 theta, r_4 phi) const
     475template<class T>
     476int_4 SphereGorski<T>::PixIndexSph(r_4 theta, r_4 phi) const
    478477
    479478//    Retourne l'indice "RING" du pixel vers lequel pointe une direction
    … …  
    481480//--
    482481{
    483   return ang2pix_ring(nSide_, double(theta), double(phi));
    484 }
    485 
    486 //++
    487 int_4 SphereGorski::PixIndexSphNest(r_4 theta, r_4 phi) const
     482  return ang2pix_ring(nSide_,double(theta),double(phi));
     483}
     484
     485//++
     486template<class T>
     487int_4 SphereGorski<T>::PixIndexSphNest(r_4 theta, r_4 phi) const
    488488
    489489//    Retourne l'indice NESTED" du pixel vers lequel pointe une direction
    … …  
    491491//--
    492492{
    493   return ang2pix_nest(nSide_, double(theta), double(phi));
     493  return ang2pix_nest(nSide_,double(theta),double(phi));
    494494}
    495495
    … …  
    497497/* --Methode-- */
    498498//++
    499 void SphereGorski::PixThetaPhi(int_4 k, r_4& teta, r_4& phi) const
     499template<class T>
     500void SphereGorski<T>::PixThetaPhi(int_4 k, r_4& teta, r_4& phi) const
    500501
    501502//    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice "RING" k
    … …  
    504505  double t;
    505506  double p;
    506   pix2ang_ring(nSide_,k, t,  p);
    507   teta=(r_4)t;
    508   phi=(r_4)p;
    509 }
    510 
    511 //++
    512 r_8       SphereGorski::PixSolAngle(int_4 dummy) const
     507  pix2ang_ring(nSide_,k, t, p);
     508  teta= (r_4)t;
     509  phi = (r_4)p;
     510}
     511
     512//++
     513template<class T>
     514r_8 SphereGorski<T>::PixSolAngle(int_4 dummy) const
    513515//    Pixel Solid angle  (steradians)
    514516//    All the pixels have the same solid angle. The dummy argument is
    … …  
    520522}
    521523
    522 
    523 //++
    524 void SphereGorski::PixThetaPhiNest(int_4 k, float& teta, float& phi)  const
     524//++
     525template<class T>
     526void SphereGorski<T>::PixThetaPhiNest(int_4 k, float& teta, float& phi)  const
    525527
    526528//    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice
    … …  
    531533  double p;
    532534  pix2ang_nest(nSide_, k, t, p);
    533   teta=(r_4)t;
    534   phi=(r_4)p;
    535 }
    536 
    537 //++
    538 int_4 SphereGorski::NestToRing(int_4 k)
     535  teta= (r_4)t;
     536  phi = (r_4)p;
     537}
     538
     539//++
     540template<class T>
     541int_4 SphereGorski<T>::NestToRing(int_4 k) const
    539542
    540543//    conversion d'index NESTD en un index RING
    … …  
    544547  return  nest2ring(nSide_,k);
    545548}
    546 //++
    547 int_4 SphereGorski::RingToNest(int_4 k)
     549
     550//++
     551template<class T>
     552int_4 SphereGorski<T>::RingToNest(int_4 k) const
    548553//
    549554//    conversion d'index RING en un index NESTED
    … …  
    553558  return  ring2nest(nSide_,k);
    554559}
    555 //++
    556 void SphereGorski::anharm(int nlmax, float sym_c,float* powspec)
     560
     561/*
     562//++
     563template<class T>
     564void SphereGorski<T>::anharm(int nlmax, float sym_c,float* powspec)
    557565//
    558566//    analyse en harmoniques spheriques des valeurs des pixels de la
    … …  
    573581//
    574582{
    575   if (nlmax > 2*nSide_) {
    576     cout << " anharm : nlmax= " << nlmax <<
    577             " doit etre <= 2*nsmax (cf. Nyquist), soit :" << 2*nSide_ << endl;
    578     exit(1);
     583    if (nlmax > 2*nSide_) {
     584          cout << " anharm : nlmax= " << nlmax <<
     585                          " doit etre <= 2*nsmax (cf. Nyquist), soit :" << 2*nSide_ << endl;
     586         exit(1);
    579587  }
    580588  else {
    … …  
    584592  sym_cut_deg_=sym_c;
    585593  float* map=new float[nPix_];
    586   for (int k=0; k<nPix_; k++) map[k]=(float)mPix_[k];
     594  for (int k=0; k<nPix_; k++) map[k]=(float)pixels_(k);
    587595  int nsmax=nSide_;
    588596  int nmmax=nmmax_;
    … …  
    607615  delete [] phas_s;
    608616  delete [] dataw;
    609   delete [] work;
    610  
    611 }
    612 
    613 
    614 //++
    615 void SphereGorski::synharm(int nlmax, int iseed,float fwhm, float* powspec)
     617  delete [] work; 
     618}
     619*/
     620
     621/*
     622//++
     623template<class T>
     624void SphereGorski<T>::synharm(int nlmax, int iseed,float fwhm, float* powspec)
    616625//
    617626//    synthese  des valeurs des pixels de la sphere par l'intermediaire
    … …  
    654663  int nsmax=nSide_;
    655664  int nmmax=nmmax_;
    656   float* alm_T=new float[2*(nlmax+1)*(nmmax+1)];
    657  
    658 
    659 //    tableaux de travail
     665  float* alm_T=new float[2*(nlmax+1)*(nmmax+1)];
     666
     667  //    tableaux de travail
    660668  double* b_north=new double[2*(2*nmmax+1)];
    661669  double* b_south=new double[2*(2*nmmax+1)];
    … …  
    666674  synfast_(nsmax,nlmax,nmmax,iseed,fwhm, map,alm_T, powspec,
    667675                  b_north,b_south,bw,data,work,lread);
    668   for (int k=0; k<nPix_; k++) mPix_[k]=(double)map[k];
     676  for (int k=0; k<nPix_; k++) pixels_(k) = (T)map[k];
    669677  delete [] map;
    670678  delete [] alm_T;
    … …  
    675683  delete [] work;
    676684  delete [] lread;
    677 
    678 }
    679 
    680 int SphereGorski::nest2ring(int nside, int ipnest) const {
     685}
     686*/
     687
     688template<class T>
     689int SphereGorski<T>::nest2ring(int nside, int ipnest) const {
    681690  /*
    682     c=======================================================================
     691    ====================================================
    683692    subroutine nest2ring(nside, ipnest, ipring)
    684     c=======================================================================
     693    ====================================================
    685694    c     conversion from NESTED to RING pixel number
    686     c=======================================================================
    687   */
    688   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    689   //  (16/12/98)
    690 
    691       int npix, npface, face_num, ncap, n_before;
    692       int ipf, ip_low, ip_trunc, ip_med, ip_hi;
    693       int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
    694       int ns_max=8192;
    695       int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
    696       int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
    697 
    698       if(  nside<1 ||  nside>ns_max ) {
    699         cout << "nside out of range" << endl;
    700         exit(0);
    701       }
    702       npix = 12 *  nside* nside;
    703       if( ipnest<0 || ipnest>npix-1 ) {
    704         cout << "ipnest out of range" << endl;
    705         exit(0);
    706       }
    707 
    708 
    709       ncap  = 2* nside*( nside-1);// ! number of points in the North Polar cap
    710       nl4   = 4* nside;
    711 
    712       //c     finds the face, and the number in the face
    713       npface =  nside* nside;
    714       //cccccc      ip = ipnest - 1         ! in {0,npix-1}
    715 
    716       face_num = ipnest/npface;//  ! face number in {0,11}
    717       ipf =ipnest%npface;//  ! pixel number in the face {0,npface-1}
    718       //c     finds the x,y on the face (starting from the lowest corner)
    719       //c     from the pixel number
    720       ip_low=ipf%1024;                //   ! content of the last 10 bits
    721       ip_trunc =   ipf/1024;         //    ! truncation of the last 10 bits
    722       ip_med=ip_trunc%1024;         //     ! content of the next 10 bits
    723       ip_hi  =     ip_trunc/1024;//   ! content of the high weight 10 bits
    724 
    725       ix = 1024*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
    726       iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + pix2y_[ip_low];
    727 
    728       //c     transforms this in (horizontal, vertical) coordinates
    729       jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
    730       jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
    731 
    732       //c     computes the z coordinate on the sphere
    733       //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
    734       jr =  jrll[face_num]*nside - jrt - 1;
    735       nr = nside;//                  ! equatorial region (the most frequent)
    736       n_before = ncap + nl4 * (jr - nside);
    737       kshift=(jr - nside)%2;
    738       if( jr<nside ) {//then     ! north pole region
    739          nr = jr;
    740          n_before = 2 * nr * (nr - 1);
    741          kshift = 0;
    742       }
    743       else if( jr>3*nside ) {//then ! south pole region
    744          nr = nl4 - jr;
    745          n_before = npix - 2 * (nr + 1) * nr;
    746          kshift = 0;
    747       }
    748 
    749       //c     computes the phi coordinate on the sphere, in [0,2Pi]
    750       jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
    751 
    752       if( jp>nl4 ) jp = jp - nl4;
    753       if( jp<1 )   jp = jp + nl4;
    754 
    755       int aux=n_before + jp - 1;
    756       return (n_before + jp - 1);// ! in {0, npix-1}
    757 
    758 }
    759 
    760 void SphereGorski::mk_pix2xy(int *pix2x,int *pix2y) {
    761   /*
    762     c=======================================================================
    763     subroutine mk_pix2xy
    764     c=======================================================================
    765     c     constructs the array giving x and y in the face from pixel number
    766     c     for the nested (quad-cube like) ordering of pixels
    767     c
    768     c     the bits corresponding to x and y are interleaved in the pixel number
    769     c     one breaks up the pixel number by even and odd bits
    770     c=======================================================================
    771   */
    772   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    773   //  (16/12/98)
    774 
    775       int kpix, jpix, IX, IY, IP, ID;
    776       for (int i=0;i<1023;i++ ) pix2x[i]=0;
    777 
    778       for( kpix=0;kpix<1024;kpix++ ) {
    779         jpix = kpix;
    780         IX = 0;
    781         IY = 0;
    782         IP = 1 ;//              ! bit position (in x and y)
    783         while( jpix!=0 ){// ! go through all the bits
    784           ID=jpix%2;//  ! bit value (in kpix), goes in ix
    785           jpix = jpix/2;
    786           IX = ID*IP+IX;
    787          
    788           ID=jpix%2;//  ! bit value (in kpix), goes in iy
    789           jpix = jpix/2;
    790           IY = ID*IP+IY;
    791          
    792           IP = 2*IP;//         ! next bit (in x and y)
    793         }
    794        
    795         pix2x[kpix] = IX;//     ! in 0,31
    796         pix2y[kpix] = IY;//     ! in 0,31
    797       }
    798 }
    799 
    800 int SphereGorski::ring2nest(int nside, int ipring) const {
    801   /*
    802     c=======================================================================
    803     subroutine ring2nest(nside, ipring, ipnest)
    804     c=======================================================================
    805     c     conversion from RING to NESTED pixel number
    806     c=======================================================================
     695    ====================================================
    807696    */
    808697  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    809698  //  (16/12/98)
     699 
     700  const PIXELS_XY& PXY= PIXELS_XY::instance();
     701
     702  int npix, npface, face_num, ncap, n_before;
     703  int ipf, ip_low, ip_trunc, ip_med, ip_hi;
     704  int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
     705  int ns_max=8192;
     706  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
     707  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
     708 
     709  if(  nside<1 ||  nside>ns_max ) {
     710    cout << "nside out of range" << endl;
     711    exit(0);
     712  }
     713  npix = 12 *  nside* nside;
     714  if( ipnest<0 || ipnest>npix-1 ) {
     715    cout << "ipnest out of range" << endl;
     716    exit(0);
     717  }
     718
     719  ncap  = 2* nside*( nside-1);// ! number of points in the North Polar cap
     720  nl4   = 4* nside;
     721 
     722  //c     finds the face, and the number in the face
     723  npface =  nside* nside;
     724  //cccccc      ip = ipnest - 1         ! in {0,npix-1}
     725 
     726  face_num = ipnest/npface;//  ! face number in {0,11}
     727  ipf =ipnest%npface;//  ! pixel number in the face {0,npface-1}
     728  //c     finds the x,y on the face (starting from the lowest corner)
     729  //c     from the pixel number
     730  ip_low=ipf%1024;                //   ! content of the last 10 bits
     731  ip_trunc =   ipf/1024;         //    ! truncation of the last 10 bits
     732  ip_med=ip_trunc%1024;         //     ! content of the next 10 bits
     733  ip_hi  =     ip_trunc/1024;//   ! content of the high weight 10 bits
     734 
     735  ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
     736  iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
     737 
     738  //c     transforms this in (horizontal, vertical) coordinates
     739  jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
     740  jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
     741 
     742  //c     computes the z coordinate on the sphere
     743  //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
     744  jr =  jrll[face_num]*nside - jrt - 1;
     745  nr = nside;//                  ! equatorial region (the most frequent)
     746  n_before = ncap + nl4 * (jr - nside);
     747  kshift=(jr - nside)%2;
     748  if( jr<nside ) {//then     ! north pole region
     749    nr = jr;
     750    n_before = 2 * nr * (nr - 1);
     751    kshift = 0;
     752  }
     753  else if( jr>3*nside ) {//then ! south pole region
     754    nr = nl4 - jr;
     755    n_before = npix - 2 * (nr + 1) * nr;
     756    kshift = 0;
     757  }
     758 
     759  //c     computes the phi coordinate on the sphere, in [0,2Pi]
     760  jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
     761 
     762  if( jp>nl4 ) jp = jp - nl4;
     763  if( jp<1 )   jp = jp + nl4;
     764 
     765  int aux=n_before + jp - 1;
     766  return (n_before + jp - 1);// ! in {0, npix-1}
     767}
     768
     769template<class T>
     770int SphereGorski<T>::ring2nest(int nside, int ipring) const
     771{
     772  /*
     773    ==================================================
     774    subroutine ring2nest(nside, ipring, ipnest)
     775    ==================================================
     776    c     conversion from RING to NESTED pixel number
     777    ==================================================
     778    */
     779  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
     780  //  (16/12/98)
     781
     782  const PIXELS_XY& PXY= PIXELS_XY::instance();
    810783
    811784  double fihip, hip;
    … …  
    898871  iy_low = (int)fmod(iy,128);
    899872  iy_hi  = iy/128;
    900   ipf =  (x2pix_[ix_hi]+y2pix_[iy_hi]) * (128 * 128)
    901     + (x2pix_[ix_low]+y2pix_[iy_low]);
     873  ipf=(PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
    902874
    903875  return (ipf + face_num* nside *nside);//   ! in {0, 12*nside**2 - 1}
    904876}
    905877
    906 void SphereGorski::mk_xy2pix(int *x2pix, int *y2pix) {
     878template<class T>
     879int SphereGorski<T>::ang2pix_ring(int nside, double theta, double phi) const
     880{
    907881  /*
    908     c=======================================================================
    909     subroutine mk_xy2pix
    910     c=======================================================================
    911     c     sets the array giving the number of the pixel lying in (x,y)
    912     c     x and y are in {1,128}
    913     c     the pixel number is in {0,128**2-1}
    914     c
    915     c     if  i-1 = sum_p=0  b_p * 2^p
    916     c     then ix = sum_p=0  b_p * 4^p
    917     c          iy = 2*ix
    918     c     ix + iy in {0, 128**2 -1}
    919     c=======================================================================
    920   */
    921   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    922   //  (16/12/98)
    923 
    924   int K,IP,I,J,ID;
    925 
    926   for( int i(0);i<127;i++ ) x2pix[i]=0;
    927   for( I=1;I<=128;I++ ) {
    928     J  = I-1;//            !pixel numbers
    929     K  = 0;//
    930     IP = 1;//
    931     truc : if( J==0 ) {
    932       x2pix[I-1] = K;
    933       y2pix[I-1] = 2*K;
    934     }
    935     else {
    936       ID = (int)fmod(J,2);
    937       J  = J/2;
    938       K  = IP*ID+K;
    939       IP = IP*4;
    940       goto truc;
    941     }
    942   }     
    943   //c      endif
    944  
    945 }
    946 
    947 int SphereGorski::ang2pix_ring(int nside, double theta, double phi) const {
    948     /*
    949     c=======================================================================
     882    ==================================================
    950883    c     gives the pixel number ipix (RING)
    951884    c     corresponding to angles theta and phi
    952     c=======================================================================
    953   */
     885    c==================================================
     886    */
    954887  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    955888  //  (16/12/98)
    … …  
    1017950  }
    1018951    return (ipix1 - 1);// ! in {0, npix-1}
    1019  
    1020 }
    1021 
    1022 int SphereGorski::ang2pix_nest(int nside, double theta, double phi) const {
     952}
     953
     954template<class T>
     955int SphereGorski<T>::ang2pix_nest(int nside, double theta, double phi) const
     956{
    1023957  /*
    1024     c=======================================================================
     958    ==================================================
    1025959    subroutine ang2pix_nest(nside, theta, phi, ipix)
    1026     c=======================================================================
     960    ==================================================
    1027961    c     gives the pixel number ipix (NESTED)
    1028962    c     corresponding to angles theta and phi
    … …  
    1033967    c     that the treatement of round-off will be consistent
    1034968    c     for every resolution
    1035     c=======================================================================
    1036   */
     969    ==================================================
     970    */
    1037971  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    1038972  //  (16/12/98)
    1039    
    1040       double    z, za, z0, tt, tp, tmp;
    1041       int face_num,jp,jm;
    1042       int ifp, ifm;
    1043       int  ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
    1044       double piover2(Pi/2.), twopi(2.*Pi);
    1045       int ns_max(8192);
    1046 
    1047       if( nside<1 || nside>ns_max ) {
    1048         cout << "nside out of range" << endl;
    1049         exit(0);
    1050       }
    1051       if( theta<0 || theta>Pi ) {
    1052         cout << "theta out of range" << endl;
    1053         exit(0);
    1054       }
    1055       z  = cos(theta);
    1056       za = fabs(z);
    1057       z0 = 2./3.;
    1058       if( phi>=twopi ) phi = phi - twopi;
    1059       if( phi<0. )    phi = phi + twopi;
    1060       tt = phi / piover2;// ! in [0,4[
    1061       if( za<=z0 ) { // then ! equatorial region
    1062        
    1063         //(the index of edge lines increase when the longitude=phi goes up)
    1064         jp = (int)floor(ns_max*(0.5 + tt - z*0.75));// !  ascending edge line index
    1065         jm = (int)floor(ns_max*(0.5 + tt + z*0.75));// ! descending edge line index
    1066 
    1067         //c        finds the face
    1068         ifp = jp / ns_max;//  ! in {0,4}
    1069         ifm = jm / ns_max;
    1070         if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then          ! faces 4 to 7
    1071         else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
    1072         else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
    1073 
    1074         ix = (int)fmod(jm, ns_max);
    1075         iy = ns_max - (int)fmod(jp, ns_max) - 1;
    1076       }
    1077       else { //! polar region, za > 2/3
    1078              
    1079         ntt = (int)floor(tt);
    1080         if( ntt>=4 ) ntt = 3;
    1081         tp = tt - ntt;
    1082         tmp = sqrt( 3.*(1. - za) );//  ! in ]0,1]
    1083 
    1084         //(the index of edge lines increase when distance from the closest pole goes up)
    1085         jp = (int)floor( ns_max * tp          * tmp );// ! line going toward the pole as phi increases
    1086         jm = (int)floor( ns_max * (1. - tp) * tmp );// ! that one goes away of the closest pole
    1087         jp = (int)min(ns_max-1, jp);// ! for points too close to the boundary
    1088         jm = (int)min(ns_max-1, jm);
    1089 
    1090          // finds the face and pixel's (x,y)
    1091          if( z>=0 ) {
    1092            face_num = ntt;//  ! in {0,3}
    1093             ix = ns_max - jm - 1;
    1094             iy = ns_max - jp - 1;
    1095          }
    1096          else {
    1097            face_num = ntt + 8;// ! in {8,11}
    1098             ix =  jp;
    1099             iy =  jm;
    1100          }
    1101       }
    1102 
    1103       ix_low = (int)fmod(ix,128);
    1104       ix_hi  =     ix/128;
    1105       iy_low = (int)fmod(iy,128);
    1106       iy_hi  =     iy/128;
    1107       ipf =  (x2pix_[ix_hi]+y2pix_[iy_hi]) * (128 * 128)+ (x2pix_[ix_low]+y2pix_[iy_low]);
    1108       ipf = ipf / pow(ns_max/nside,2);//  ! in {0, nside**2 - 1}
    1109       return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
    1110 }
    1111 
    1112 
    1113 void SphereGorski::pix2ang_ring(int nside, int ipix, double& theta,  double& phi) const {
     973 
     974  const PIXELS_XY& PXY= PIXELS_XY::instance();
     975
     976  double    z, za, z0, tt, tp, tmp;
     977  int face_num,jp,jm;
     978  int ifp, ifm;
     979  int  ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
     980  double piover2(Pi/2.), twopi(2.*Pi);
     981  int ns_max(8192);
     982 
     983  if( nside<1 || nside>ns_max ) {
     984    cout << "nside out of range" << endl;
     985    exit(0);
     986  }
     987  if( theta<0 || theta>Pi ) {
     988    cout << "theta out of range" << endl;
     989    exit(0);
     990  }
     991  z  = cos(theta);
     992  za = fabs(z);
     993  z0 = 2./3.;
     994  if( phi>=twopi ) phi = phi - twopi;
     995  if( phi<0. )    phi = phi + twopi;
     996  tt = phi / piover2;// ! in [0,4[
     997  if( za<=z0 ) { // then ! equatorial region
     998   
     999    //(the index of edge lines increase when the longitude=phi goes up)
     1000    jp = (int)floor(ns_max*(0.5 + tt - z*0.75));// !  ascending edge line index
     1001    jm = (int)floor(ns_max*(0.5 + tt + z*0.75));// ! descending edge line index
     1002   
     1003    //c        finds the face
     1004    ifp = jp / ns_max;//  ! in {0,4}
     1005    ifm = jm / ns_max;
     1006    if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then          ! faces 4 to 7
     1007    else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
     1008    else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
     1009   
     1010    ix = (int)fmod(jm, ns_max);
     1011    iy = ns_max - (int)fmod(jp, ns_max) - 1;
     1012  }
     1013  else { //! polar region, za > 2/3
     1014   
     1015    ntt = (int)floor(tt);
     1016    if( ntt>=4 ) ntt = 3;
     1017    tp = tt - ntt;
     1018    tmp = sqrt( 3.*(1. - za) );//  ! in ]0,1]
     1019   
     1020    //(the index of edge lines increase when distance from the closest pole goes up)
     1021    jp = (int)floor(ns_max*tp*tmp); // ! line going toward the pole as phi increases
     1022    jm = (int)floor(ns_max*(1.-tp)*tmp); // ! that one goes away of the closest pole
     1023    jp = (int)min(ns_max-1, jp);// ! for points too close to the boundary
     1024    jm = (int)min(ns_max-1, jm);
     1025   
     1026    // finds the face and pixel's (x,y)
     1027    if( z>=0 ) {
     1028      face_num = ntt;//  ! in {0,3}
     1029      ix = ns_max - jm - 1;
     1030      iy = ns_max - jp - 1;
     1031    }
     1032    else {
     1033      face_num = ntt + 8;// ! in {8,11}
     1034      ix =  jp;
     1035      iy =  jm;
     1036    }
     1037  }
     1038 
     1039  ix_low = (int)fmod(ix,128);
     1040  ix_hi  =     ix/128;
     1041  iy_low = (int)fmod(iy,128);
     1042  iy_hi  =     iy/128;
     1043  ipf= (PXY.x2pix_(ix_hi)+PXY.y2pix_(iy_hi))*(128*128)+(PXY.x2pix_(ix_low)+PXY.y2pix_(iy_low));
     1044  ipf = ipf / pow(ns_max/nside,2);//  ! in {0, nside**2 - 1}
     1045  return ( ipf + face_num*pow(nside,2));//    ! in {0, 12*nside**2 - 1}
     1046}
     1047
     1048template<class T>
     1049void SphereGorski<T>::pix2ang_ring(int nside,int ipix,double& theta,double& phi) const {
    11141050  /*
    1115     c=======================================================================
     1051    ===================================================
    11161052    c     gives theta and phi corresponding to pixel ipix (RING)
    11171053    c     for a parameter nside
    1118     c=======================================================================
    1119   */
     1054    ===================================================
     1055    */
    11201056  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    11211057  //  (16/12/98)
    … …  
    11781114}
    11791115
    1180 void SphereGorski::pix2ang_nest(int nside, int ipix, double& theta, double& phi) const {
     1116template<class T>
     1117void SphereGorski<T>::pix2ang_nest(int nside,int ipix,double& theta,double& phi) const {
    11811118  /*
    1182     c=======================================================================
     1119    ==================================================
    11831120    subroutine pix2ang_nest(nside, ipix, theta, phi)
    1184     c=======================================================================
     1121    ==================================================
    11851122    c     gives theta and phi corresponding to pixel ipix (NESTED)
    11861123    c     for a parameter nside
    1187     c=======================================================================
    1188   */
     1124    ==================================================
     1125    */
    11891126  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    11901127  //  (16/12/98)
    1191    
    1192       int npix, npface, face_num;
    1193       int  ipf, ip_low, ip_trunc, ip_med, ip_hi;
    1194       int     ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
    1195       double z, fn, fact1, fact2;
    1196       double piover2(Pi/2.);
    1197       int ns_max(8192);
    1198      
    1199      
    1200       // ! coordinate of the lowest corner of each face
    1201       int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
    1202       int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
    1203      
    1204       if( nside<1 || nside>ns_max ) {
    1205         cout << "nside out of range" << endl;
    1206         exit(0);
    1207       }
    1208       npix = 12 * nside*nside;
    1209       if( ipix<0 || ipix>npix-1 ) {
    1210         cout << "ipix out of range" << endl;
    1211         exit(0);
    1212       }
    1213 
    1214       fn = 1.*nside;
    1215       fact1 = 1./(3.*fn*fn);
    1216       fact2 = 2./(3.*fn);
    1217       nl4   = 4*nside;
    1218 
    1219       //c     finds the face, and the number in the face
    1220       npface = nside*nside;
    1221 
    1222       face_num = ipix/npface;//  ! face number in {0,11}
    1223       ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
    1224 
    1225       //c     finds the x,y on the face (starting from the lowest corner)
    1226       //c     from the pixel number
    1227       ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
    1228       ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
    1229       ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
    1230       ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
    1231 
    1232       ix = 1024*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
    1233       iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + pix2y_[ip_low];
    1234 
    1235       //c     transforms this in (horizontal, vertical) coordinates
    1236       jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
    1237       jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
    1238 
    1239       //c     computes the z coordinate on the sphere
    1240       //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
    1241       jr =  jrll[face_num]*nside - jrt - 1;
    1242       nr = nside;//                  ! equatorial region (the most frequent)
    1243       z  = (2*nside-jr)*fact2;
    1244       kshift = (int)fmod(jr - nside, 2);
    1245       if( jr<nside ) { //then     ! north pole region
    1246          nr = jr;
    1247          z = 1. - nr*nr*fact1;
    1248          kshift = 0;
    1249       }
    1250       else {
    1251         if( jr>3*nside ) {// then ! south pole region
    1252          nr = nl4 - jr;
    1253          z = - 1. + nr*nr*fact1;
    1254          kshift = 0;
     1128
     1129  const PIXELS_XY& PXY= PIXELS_XY::instance();
     1130   
     1131  int npix, npface, face_num;
     1132  int ipf, ip_low, ip_trunc, ip_med, ip_hi;
     1133  int ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
     1134  double z, fn, fact1, fact2;
     1135  double piover2(Pi/2.);
     1136  int ns_max(8192);
     1137         
     1138  // ! coordinate of the lowest corner of each face
     1139  int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
     1140  int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
     1141 
     1142  if( nside<1 || nside>ns_max ) {
     1143    cout << "nside out of range" << endl;
     1144    exit(0);
     1145  }
     1146  npix = 12 * nside*nside;
     1147  if( ipix<0 || ipix>npix-1 ) {
     1148    cout << "ipix out of range" << endl;
     1149    exit(0);
     1150  }
     1151 
     1152  fn = 1.*nside;
     1153  fact1 = 1./(3.*fn*fn);
     1154  fact2 = 2./(3.*fn);
     1155  nl4   = 4*nside;
     1156
     1157  //c     finds the face, and the number in the face
     1158  npface = nside*nside;
     1159 
     1160  face_num = ipix/npface;//  ! face number in {0,11}
     1161  ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
     1162 
     1163  //c     finds the x,y on the face (starting from the lowest corner)
     1164  //c     from the pixel number
     1165  ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
     1166  ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
     1167  ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
     1168  ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
     1169 
     1170  ix = 1024*PXY.pix2x_(ip_hi)+32*PXY.pix2x_(ip_med)+PXY.pix2x_(ip_low);
     1171  iy = 1024*PXY.pix2y_(ip_hi)+32*PXY.pix2y_(ip_med)+PXY.pix2y_(ip_low);
     1172 
     1173  //c     transforms this in (horizontal, vertical) coordinates
     1174  jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
     1175  jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
     1176 
     1177  //c     computes the z coordinate on the sphere
     1178  //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
     1179  jr =  jrll[face_num]*nside - jrt - 1;
     1180  nr = nside;//                  ! equatorial region (the most frequent)
     1181  z  = (2*nside-jr)*fact2;
     1182  kshift = (int)fmod(jr - nside, 2);
     1183  if( jr<nside ) { //then     ! north pole region
     1184    nr = jr;
     1185    z = 1. - nr*nr*fact1;
     1186    kshift = 0;
     1187  }
     1188  else {
     1189    if( jr>3*nside ) {// then ! south pole region
     1190      nr = nl4 - jr;
     1191      z = - 1. + nr*nr*fact1;
     1192      kshift = 0;
     1193    }
     1194  }
     1195  theta = acos(z);
     1196 
     1197  //c     computes the phi coordinate on the sphere, in [0,2Pi]
     1198  //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
     1199  jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
     1200  if( jp>nl4 ) jp = jp - nl4;
     1201  if( jp<1 )   jp = jp + nl4;
     1202  phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
     1203}
     1204
     1205// retourne le nom du fichier qui contient le spectre de puissance
     1206template<class T>
     1207void SphereGorski<T>::powfile(char filename[]) const
     1208{
     1209  bool status = false;
     1210  for (int k=0; k< 128; k++)
     1211    {
     1212      if( powFile_[k] != ' ' )
     1213        {
     1214          status = true;
     1215          break;
    12551216        }
    1256       }
    1257       theta = acos(z);
    1258      
    1259       //c     computes the phi coordinate on the sphere, in [0,2Pi]
    1260       //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
    1261       jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
    1262       if( jp>nl4 ) jp = jp - nl4;
    1263       if( jp<1 )   jp = jp + nl4;
    1264 
    1265       phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
    1266 
    1267 }
    1268 
    1269 
    1270 void SphereGorski::Pix2XY::pix2ang_nest(int nside, int ipix, double& theta, double& phi) const {
    1271   /*
    1272     c=======================================================================
    1273     subroutine pix2ang_nest(nside, ipix, theta, phi)
    1274     c=======================================================================
    1275     c     gives theta and phi corresponding to pixel ipix (NESTED)
    1276     c     for a parameter nside
    1277     c=======================================================================
    1278   */
    1279   // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    1280   //  (16/12/98)
    1281    
    1282       int npix, npface, face_num;
    1283       int  ipf, ip_low, ip_trunc, ip_med, ip_hi;
    1284       int     ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4;
    1285       double z, fn, fact1, fact2;
    1286       double piover2(Pi/2.);
    1287       int ns_max(8192);
    1288      
    1289      
    1290       // ! coordinate of the lowest corner of each face
    1291       int jrll[12]={2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};//! in unit of nside
    1292       int jpll[12]={1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};// ! in unit of nside/2
    1293      
    1294       if( nside<1 || nside>ns_max ) {
    1295         cout << "nside out of range" << endl;
    1296         exit(0);
    1297       }
    1298       npix = 12 * nside*nside;
    1299       if( ipix<0 || ipix>npix-1 ) {
    1300         cout << "ipix out of range" << endl;
    1301         exit(0);
    1302       }
    1303 
    1304       fn = 1.*nside;
    1305       fact1 = 1./(3.*fn*fn);
    1306       fact2 = 2./(3.*fn);
    1307       nl4   = 4*nside;
    1308 
    1309       //c     finds the face, and the number in the face
    1310       npface = nside*nside;
    1311 
    1312       face_num = ipix/npface;//  ! face number in {0,11}
    1313       ipf = (int)fmod(ipix,npface);//  ! pixel number in the face {0,npface-1}
    1314 
    1315       //c     finds the x,y on the face (starting from the lowest corner)
    1316       //c     from the pixel number
    1317       ip_low = (int)fmod(ipf,1024);//       ! content of the last 10 bits
    1318       ip_trunc =   ipf/1024 ;//       ! truncation of the last 10 bits
    1319       ip_med = (int)fmod(ip_trunc,1024);//  ! content of the next 10 bits
    1320       ip_hi  =     ip_trunc/1024   ;//! content of the high weight 10 bits
    1321 
    1322       ix = 1024*pix2x_[ip_hi] + 32*pix2x_[ip_med] + pix2x_[ip_low];
    1323       iy = 1024*pix2y_[ip_hi] + 32*pix2y_[ip_med] + pix2y_[ip_low];
    1324 
    1325       //c     transforms this in (horizontal, vertical) coordinates
    1326       jrt = ix + iy;//  ! 'vertical' in {0,2*(nside-1)}
    1327       jpt = ix - iy;//  ! 'horizontal' in {-nside+1,nside-1}
    1328 
    1329       //c     computes the z coordinate on the sphere
    1330       //      jr =  jrll[face_num+1]*nside - jrt - 1;//   ! ring number in {1,4*nside-1}
    1331       jr =  jrll[face_num]*nside - jrt - 1;
    1332       nr = nside;//                  ! equatorial region (the most frequent)
    1333       z  = (2*nside-jr)*fact2;
    1334       kshift = (int)fmod(jr - nside, 2);
    1335       if( jr<nside ) { //then     ! north pole region
    1336          nr = jr;
    1337          z = 1. - nr*nr*fact1;
    1338          kshift = 0;
    1339       }
    1340       else {
    1341         if( jr>3*nside ) {// then ! south pole region
    1342          nr = nl4 - jr;
    1343          z = - 1. + nr*nr*fact1;
    1344          kshift = 0;
    1345         }
    1346       }
    1347       theta = acos(z);
    1348      
    1349       //c     computes the phi coordinate on the sphere, in [0,2Pi]
    1350       //      jp = (jpll[face_num+1]*nr + jpt + 1 + kshift)/2;//  ! 'phi' number in the ring in {1,4*nr}
    1351       jp = (jpll[face_num]*nr + jpt + 1 + kshift)/2;
    1352       if( jp>nl4 ) jp = jp - nl4;
    1353       if( jp<1 )   jp = jp + nl4;
    1354 
    1355       phi = (jp - (kshift+1)*0.5) * (piover2 / nr);
    1356 
    1357 }
    1358 
    1359 
    1360 
     1217    }
     1218  if ( status )
     1219    {
     1220      strcpy(filename,powFile_);
     1221    }
     1222  else
     1223    {
     1224      strcpy(filename,"no file");
     1225    }
     1226}
     1227
     1228template<class T>
     1229void SphereGorski<T>::getParafast(int_4& nlmax,int_4& nmmax,int_4& iseed,float& fwhm,float& quadr,float& cut) const
     1230{
     1231  nlmax= nlmax_;
     1232  nmmax= nmmax_;
     1233  iseed= iseed_;
     1234  fwhm = fwhm_;
     1235  quadr= quadrupole_;
     1236  cut  = sym_cut_deg_;
     1237}
     1238
     1239template<class T>
     1240void SphereGorski<T>::setParafast(int_4 nlmax,int_4 nmmax,int_4 iseed,float fwhm,float quadr,float cut,char* filename)
     1241{
     1242  nlmax_= nlmax;
     1243  nmmax_= nmmax;
     1244  iseed_= iseed;
     1245  fwhm_ = fwhm;
     1246  quadrupole_ = quadr;
     1247  sym_cut_deg_= cut;
     1248  strcpy(powFile_,filename);
     1249}
     1250
     1251template <class T>
     1252void SphereGorski<T>::print(ostream& os) const
     1253{
     1254  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
     1255  //
     1256  os << " nSide_ = " << nSide_  << endl;
     1257  os << " nPix_   = " << nPix_   << endl;
     1258  os << " omeg_ = " << omeg_   << endl;
     1259
     1260  os << " contenu de pixels : ";
     1261  for(int i=0; i < nPix_; i++)
     1262    {
     1263      if(i%5 == 0) os << endl;
     1264      os <<  pixels_(i) <<", ";
     1265    }
     1266  os << endl;
     1267
     1268  os << endl;
     1269  const PIXELS_XY& PXY= PIXELS_XY::instance();
     1270
     1271  os << endl;  os << " contenu des tableaux conversions "<<endl;
     1272  for(int i=0; i < 5; i++)
     1273    {
     1274      os<<PXY.pix2x_(i)<<", "<<PXY.pix2y_(i)<<", "<<PXY.x2pix_(i)<<", "<<PXY.y2pix_(i)<<endl;
     1275    }
     1276  os << endl;
     1277
     1278  os << " les parametres des modules anafast et synfast " <<endl;
     1279  os << " nlmax, nmmax & iseed= " <<nlmax_<<", "<<nmmax_<<", "<<iseed_<<endl;
     1280  os << " fwhm, quadr & cut = "<<fwhm_<<", "<<quadrupole_<<", "<<sym_cut_deg_<<endl;
     1281  os << " powfile= " << powFile_<<endl;
     1282}
     1283
     1284//*******************************************************************
     1285// Class FIO_SphereGorski<T>
     1286//  Les objets delegues pour la gestion de persistance
     1287//*******************************************************************
     1288
     1289template <class T>
     1290FIO_SphereGorski<T>::FIO_SphereGorski()
     1291{
     1292  dobj= new SphereGorski<T>;
     1293  ownobj= true;
     1294}
     1295
     1296template <class T>
     1297FIO_SphereGorski<T>::FIO_SphereGorski(string const& filename)
     1298{
     1299  dobj= new SphereGorski<T>;
     1300  ownobj= true;
     1301  Read(filename);
     1302}
     1303
     1304template <class T>
     1305FIO_SphereGorski<T>::FIO_SphereGorski(const SphereGorski<T>& obj)
     1306{
     1307  dobj= new SphereGorski<T>(obj);
     1308  ownobj= true;
     1309}
     1310
     1311template <class T>
     1312FIO_SphereGorski<T>::FIO_SphereGorski(SphereGorski<T>* obj)
     1313{
     1314  dobj= obj;
     1315  ownobj= false;
     1316}
     1317
     1318template <class T>
     1319FIO_SphereGorski<T>::~FIO_SphereGorski()
     1320{
     1321  if (ownobj && dobj) delete dobj;
     1322}
     1323
     1324template <class T>
     1325AnyDataObj* FIO_SphereGorski<T>::DataObj()
     1326{
     1327  return(dobj);
     1328}
     1329
     1330template <class T>
     1331void FIO_SphereGorski<T>::ReadSelf(PInPersist& is)
     1332{
     1333  cout << " FIO_SphereGorski:: ReadSelf " << endl;
     1334
     1335  if(dobj == NULL)
     1336    {
     1337      dobj= new SphereGorski<T>;
     1338    }
     1339
     1340  // Pour savoir s'il y avait un DVList Info associe
     1341  char strg[256];
     1342  is.GetLine(strg, 255);
     1343  bool hadinfo= false;
     1344  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
     1345  if(hadinfo)
     1346    {    // Lecture eventuelle du DVList Info
     1347      is >> dobj->Info();
     1348    }
     1349
     1350  int_4 nSide;
     1351  is.GetI4(nSide);
     1352  dobj->setSizeIndex(nSide);
     1353
     1354  int_4 nPix;
     1355  is.GetI4(nPix);
     1356  dobj->setNbPixels(nPix);
     1357
     1358  r_8 Omega;
     1359  is.GetR8(Omega);
     1360  dobj->setPixSolAngle(Omega);
     1361
     1362  T* pixels= new T[nPix];
     1363  PIOSReadArray(is, pixels, nPix);
     1364  dobj->setDataBlock(pixels, nPix);
     1365  delete [] pixels;
     1366
     1367  int_4 nlmax,nmmax,iseed;
     1368  is.GetI4(nlmax);
     1369  is.GetI4(nmmax);
     1370  is.GetI4(iseed);
     1371 
     1372  float fwhm,quadr,cut;
     1373  is.GetR4(fwhm);
     1374  is.GetR4(quadr);
     1375  is.GetR4(cut);
     1376
     1377  char powfl[128];
     1378  is.GetLine(powfl, 127);
     1379
     1380  dobj->setParafast(nlmax,nmmax,iseed,fwhm,quadr,cut,powfl);
     1381}
     1382
     1383template <class T>
     1384void FIO_SphereGorski<T>::WriteSelf(POutPersist& os) const
     1385{
     1386  cout << " FIO_SphereGorski:: WriteSelf " << endl;
     1387
     1388  if(dobj == NULL)
     1389    {
     1390      cout << " WriteSelf:: dobj= null " << endl;
     1391      return;
     1392    }
     1393
     1394  char strg[256];
     1395  int_4 nSide= dobj->SizeIndex();
     1396  int_4 nPix = dobj->NbPixels();
     1397 
     1398  if(dobj->ptrInfo())
     1399    {
     1400      sprintf(strg,"SphereGorski: NSide=%6d  NPix=%9d HasInfo",nSide,nPix);
     1401      os.PutLine(strg);
     1402      os << dobj->Info();
     1403    }
     1404  else
     1405    {
     1406      sprintf(strg,"SphereGorski: NSide=%6d  NPix=%9d ",nSide,nPix);
     1407      os.PutLine(strg); 
     1408    }
     1409
     1410  os.PutI4(nSide);
     1411  os.PutI4(nPix);
     1412  os.PutR8(dobj->PixSolAngle(0));
     1413
     1414  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), nPix);
     1415
     1416  int_4 nlmax,nmmax,iseed;
     1417  float fwhm,quadr,cut;
     1418  dobj->getParafast(nlmax,nmmax,iseed,fwhm,quadr,cut);
     1419  os.PutI4(nlmax);
     1420  os.PutI4(nmmax);
     1421  os.PutI4(iseed);
     1422  os.PutR4(fwhm);
     1423  os.PutR4(quadr);
     1424  os.PutR4(cut);
     1425
     1426  char powfl[128];
     1427  dobj->powfile(powfl);
     1428  os.PutLine(powfl);
     1429}
     1430
     1431#ifdef __CXX_PRAGMA_TEMPLATES__
     1432#pragma define_template SphereGorski<double>
     1433#pragma define_template SphereGorski<float>
     1434#pragma define_template SphereGorski< complex<float> >
     1435#pragma define_template SphereGorski< complex<double> >
     1436#pragma define_template FIO_SphereGorski<double>
     1437#pragma define_template FIO_SphereGorski<float>
     1438#pragma define_template FIO_SphereGorski< complex<float> >
     1439#pragma define_template FIO_SphereGorski< complex<double> >
     1440#endif
     1441#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
     1442template class SphereGorski<double>;
     1443template class SphereGorski<float>;
     1444template class SphereGorski< complex<float> >;
     1445template class SphereGorski< complex<double> >;
     1446template class FIO_SphereGorski<double>;
     1447template class FIO_SphereGorski<float>;
     1448template class FIO_SphereGorski< complex<float> >;
     1449template class FIO_SphereGorski< complex<double> >;
     1450#endif
     1451
Note: See TracChangeset for help on using the changeset viewer.