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


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

modifs francois : passage en double etc. 18-OCT-99

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/SophyaLib/Samba/spheregorski.cc

    r470 r473  
    1 //
    21#include "spheregorski.h"
    32#include "strutil.h"
     
    1413extern "C"
    1514{
    16 void anafast_(int&, int&, int&,double&,float*,float*,float*,float*,float*,float*,float*);
    17 void synfast_(int&, int&, int&,int&, float&,float*,float*,float*,double*, double*,double*,double*,double*,float*);
     15void anafast_(int&,int&,int&,double&,float*,float*,float*,float*,float*,float*,float*);
     16void synfast_(int&,int&,int&,int&,float&,float*,float*,float*,double*,double*,double*,double*,double*,float*);
    1817}
    1918
    2019//*******************************************************************
    2120// Class PIXELS_XY
    22 //  Construction des tableaux necessaires a la traduction des indices RING en
    23 //  indices NESTED (ou l'inverse)
     21// Construction des tableaux necessaires a la traduction des indices RING en
     22// indices NESTED (ou l'inverse)
    2423//*******************************************************************
    2524
     
    5756    c     one breaks up the pixel number by even and odd bits
    5857    ==================================================
    59     */
     58  */
    6059  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    6160  //  (16/12/98)
     
    101100    c     ix + iy in {0, 128**2 -1}
    102101    =================================================
    103     */
     102  */
    104103  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    105104  //  (16/12/98)
     
    196195//++
    197196template<class T>
    198 SphereGorski<T>::SphereGorski(int_4 m)
     197SphereGorski<T>::SphereGorski(int m)
    199198
    200199//    Constructeur : m est la variable nside de l'algorithme de Gorski
     
    211210    }
    212211  // verifier que m est une puissance de deux
    213   int_4 x=m;
     212  int x= m;
    214213  while(x%2 == 0) x/=2;
    215214  if(x != 1)
     
    261260//++
    262261template<class T>
    263 void SphereGorski<T>::Resize(int_4 m)
     262void SphereGorski<T>::Resize(int m)
    264263
    265264//    m est la variable nside de l'algorithme de Gorski
     
    273272  }
    274273  // verifier que m est une puissance de deux
    275   int_4 x=m;
     274  int x= m;
    276275  while (x%2==0) x/=2;
    277276  if(x != 1)
     
    285284
    286285template<class T>
    287 void  SphereGorski<T>::Pixelize( int_4 m)
     286void  SphereGorski<T>::Pixelize( int m)
    288287
    289288//    prépare la pixelisation Gorski (m a la même signification
     
    294293{
    295294  // On memorise les arguments d'appel
    296   nSide_ = m; 
     295  nSide_= m; 
    297296
    298297  // Nombre total de pixels sur la sphere entiere
    299   nPix_=12*nSide_*nSide_;
     298  nPix_= 12*nSide_*nSide_;
    300299
    301300  // pour le moment les tableaux qui suivent seront ranges dans l'ordre
     
    310309
    311310  // solid angle per pixel   
    312   omeg_= 4*Pi/nPix_;
     311  omeg_= 4.0*Pi/nPix_;
    313312}
    314313
     
    335334//++
    336335template<class T>
    337 int_4 SphereGorski<T>::NbPixels() const
     336int SphereGorski<T>::NbPixels() const
    338337
    339338//    Retourne le nombre de pixels du découpage
     
    345344//++
    346345template<class T>
    347 int_4 SphereGorski<T>::NbThetaSlices() const
     346int SphereGorski<T>::NbThetaSlices() const
    348347
    349348//    Retourne le nombre de tranches en theta sur la sphere
    350349//--
    351350{
    352   return int_4(4*nSide_-1);
    353 }
    354 
    355 //++
    356 template<class T>
    357 void  SphereGorski<T>::GetThetaSlice(int_4 index, r_4& theta, TVector<float>& phi, TVector<T>& value) const
     351  return int(4*nSide_-1);
     352}
     353
     354//++
     355template<class T>
     356void SphereGorski<T>::GetThetaSlice(int index,double& theta,TVector<double>& phi,TVector<T>& value) const
    358357
    359358//    Retourne, pour la tranche en theta d'indice 'index' le theta
     
    372371    }
    373372
    374   int_4 iring= 0;
     373  int iring= 0;
    375374  int lring  = 0;
    376375  if(index < nSide_-1)
     
    393392  phi.ReSize(lring);
    394393  value.ReSize(lring);
    395   float TH=0.;
    396   float F =0.;
     394  double TH= 0.;
     395  double FI= 0.;
    397396  for(int kk = 0; kk < lring;kk++)
    398397    {
    399       PixThetaPhi(kk+iring,TH,F);
    400       phi(kk)= F;
     398      PixThetaPhi(kk+iring,TH,FI);
     399      phi(kk)= FI;
    401400      value(kk)= PixVal(kk+iring);
    402401    }
     
    407406//++
    408407template<class T>
    409 T& SphereGorski<T>::PixVal(int_4 k)
     408T& SphereGorski<T>::PixVal(int k)
    410409
    411410//    Retourne la valeur du contenu du pixel d'indice "RING" k
     
    424423//++
    425424template<class T>
    426 T const& SphereGorski<T>::PixVal(int_4 k) const
     425T const& SphereGorski<T>::PixVal(int k) const
    427426
    428427//    Retourne la valeur du contenu du pixel d'indice "RING" k
     
    440439//++
    441440template<class T>
    442 T& SphereGorski<T>::PixValNest(int_4 k)
     441T& SphereGorski<T>::PixValNest(int k)
    443442
    444443//    Retourne la valeur du contenu du pixel d'indice "NESTED" k
     
    456455
    457456template<class T>
    458 T const& SphereGorski<T>::PixValNest(int_4 k) const
     457T const& SphereGorski<T>::PixValNest(int k) const
    459458
    460459//    Retourne la valeur du contenu du pixel d'indice "NESTED" k
     
    467466      THROW(rangeCheckErr);
    468467  }
    469   int_4 pix= nest2ring(nSide_,k);
     468  int pix= nest2ring(nSide_,k);
    470469  return *(pixels_.Data()+pix);
    471470}
     
    474473//++
    475474template<class T>
    476 int_4 SphereGorski<T>::PixIndexSph(r_4 theta, r_4 phi) const
     475int SphereGorski<T>::PixIndexSph(double theta,double phi) const
    477476
    478477//    Retourne l'indice "RING" du pixel vers lequel pointe une direction
     
    480479//--
    481480{
    482   return ang2pix_ring(nSide_,double(theta),double(phi));
    483 }
    484 
    485 //++
    486 template<class T>
    487 int_4 SphereGorski<T>::PixIndexSphNest(r_4 theta, r_4 phi) const
     481  return ang2pix_ring(nSide_,theta,phi);
     482}
     483
     484//++
     485template<class T>
     486int SphereGorski<T>::PixIndexSphNest(double theta,double phi) const
    488487
    489488//    Retourne l'indice NESTED" du pixel vers lequel pointe une direction
     
    491490//--
    492491{
    493   return ang2pix_nest(nSide_,double(theta),double(phi));
     492  return ang2pix_nest(nSide_,theta,phi);
    494493}
    495494
     
    498497//++
    499498template<class T>
    500 void SphereGorski<T>::PixThetaPhi(int_4 k, r_4& teta, r_4& phi) const
    501 
    502 //    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice "RING" k
    503 //--
    504 {
    505   double t;
    506   double p;
    507   pix2ang_ring(nSide_,k, t, p);
    508   teta= (r_4)t;
    509   phi = (r_4)p;
    510 }
    511 
    512 //++
    513 template<class T>
    514 r_8 SphereGorski<T>::PixSolAngle(int_4 dummy) const
     499void SphereGorski<T>::PixThetaPhi(int k,double& theta,double& phi) const
     500
     501//   Retourne les coordonnées (theta,phi) du milieu du pixel d'indice "RING" k
     502//--
     503{
     504  pix2ang_ring(nSide_,k,theta,phi);
     505}
     506
     507//++
     508template<class T>
     509double SphereGorski<T>::PixSolAngle(int dummy) const
    515510//    Pixel Solid angle  (steradians)
    516511//    All the pixels have the same solid angle. The dummy argument is
     
    524519//++
    525520template<class T>
    526 void SphereGorski<T>::PixThetaPhiNest(int_4 k, float& teta, float& phi) const
    527 
    528 //    Retourne les coordonnées (teta,phi) du milieu du pixel d'indice
     521void SphereGorski<T>::PixThetaPhiNest(int k,double& theta,double& phi) const
     522
     523//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice
    529524//    NESTED k
    530525//--
    531526{   
    532   double t;
    533   double p;
    534   pix2ang_nest(nSide_, k, t, p);
    535   teta= (r_4)t;
    536   phi = (r_4)p;
    537 }
    538 
    539 //++
    540 template<class T>
    541 int_4 SphereGorski<T>::NestToRing(int_4 k) const
     527  pix2ang_nest(nSide_,k,theta,phi);
     528}
     529
     530//++
     531template<class T>
     532int SphereGorski<T>::NestToRing(int k) const
    542533
    543534//    conversion d'index NESTD en un index RING
     
    550541//++
    551542template<class T>
    552 int_4 SphereGorski<T>::RingToNest(int_4 k) const
     543int SphereGorski<T>::RingToNest(int k) const
    553544//
    554545//    conversion d'index RING en un index NESTED
     
    687678
    688679template<class T>
    689 int SphereGorski<T>::nest2ring(int nside, int ipnest) const {
     680int SphereGorski<T>::nest2ring(int nside, int ipnest) const
     681{
    690682  /*
    691683    ====================================================
     
    694686    c     conversion from NESTED to RING pixel number
    695687    ====================================================
    696     */
     688  */
    697689  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    698690  //  (16/12/98)
     
    776768    c     conversion from RING to NESTED pixel number
    777769    ==================================================
    778     */
     770  */
    779771  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    780772  //  (16/12/98)
     
    884876    c     corresponding to angles theta and phi
    885877    c==================================================
    886     */
     878  */
    887879  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    888880  //  (16/12/98)
     
    968960    c     for every resolution
    969961    ==================================================
    970     */
     962  */
    971963  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    972964  //  (16/12/98)
     
    1004996    ifp = jp / ns_max;//  ! in {0,4}
    1005997    ifm = jm / ns_max;
    1006     if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then          ! faces 4 to 7
     998    if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; //then  ! faces 4 to 7
    1007999    else if( ifp<ifm ) face_num = (int)fmod(ifp,4); // (half-)faces 0 to 3
    10081000    else face_num = (int)fmod(ifm,4) + 8;//! (half-)faces 8 to 11
     
    10531045    c     for a parameter nside
    10541046    ===================================================
    1055     */
     1047  */
    10561048  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    10571049  //  (16/12/98)
     
    11231115    c     for a parameter nside
    11241116    ==================================================
    1125     */
     1117  */
    11261118  // tranlated from FORTRAN (Gorski) to C, by B. Revenu, revised Guy Le Meur
    11271119  //  (16/12/98)
     
    12271219
    12281220template<class T>
    1229 void SphereGorski<T>::getParafast(int_4& nlmax,int_4& nmmax,int_4& iseed,float& fwhm,float& quadr,float& cut) const
     1221void SphereGorski<T>::getParafast(int& nlmax,int& nmmax,int& iseed,float& fwhm,float& quadr,float& cut) const
    12301222{
    12311223  nlmax= nlmax_;
     
    12381230
    12391231template<class T>
    1240 void SphereGorski<T>::setParafast(int_4 nlmax,int_4 nmmax,int_4 iseed,float fwhm,float quadr,float cut,char* filename)
     1232void SphereGorski<T>::setParafast(int nlmax,int nmmax,int iseed,float fwhm,float quadr,float cut,char* filename)
    12411233{
    12421234  nlmax_= nlmax;
     
    13481340    }
    13491341
    1350   int_4 nSide;
     1342  int nSide;
    13511343  is.GetI4(nSide);
    13521344  dobj->setSizeIndex(nSide);
    13531345
    1354   int_4 nPix;
     1346  int nPix;
    13551347  is.GetI4(nPix);
    13561348  dobj->setNbPixels(nPix);
    13571349
    1358   r_8 Omega;
     1350  double Omega;
    13591351  is.GetR8(Omega);
    13601352  dobj->setPixSolAngle(Omega);
     
    13651357  delete [] pixels;
    13661358
    1367   int_4 nlmax,nmmax,iseed;
     1359  int nlmax,nmmax,iseed;
    13681360  is.GetI4(nlmax);
    13691361  is.GetI4(nmmax);
     
    13931385
    13941386  char strg[256];
    1395   int_4 nSide= dobj->SizeIndex();
    1396   int_4 nPix = dobj->NbPixels();
     1387  int nSide= dobj->SizeIndex();
     1388  int nPix = dobj->NbPixels();
    13971389 
    13981390  if(dobj->ptrInfo())
     
    14141406  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), nPix);
    14151407
    1416   int_4 nlmax,nmmax,iseed;
     1408  int nlmax,nmmax,iseed;
    14171409  float fwhm,quadr,cut;
    14181410  dobj->getParafast(nlmax,nmmax,iseed,fwhm,quadr,cut);
Note: See TracChangeset for help on using the changeset viewer.