Ignore:
Timestamp:
Nov 13, 2001, 4:29:59 PM (24 years ago)
Author:
lemeur
Message:

adaptation aus spheres theta-phi

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/Samba/sphericaltransformserver.cc

    r1715 r1756  
    154154{
    155155  /*=======================================================================
    156     computes a map form its alm for the HEALPIX pixelisation
     156    computes a map from its alm for the HEALPIX pixelisation
    157157    map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
    158158    = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
    … …  
    167167  int_4 nmmax=nlmax;
    168168  int_4 nsmax=0;
     169  // le Resize est suppose mettre a zero
    169170  map.Resize(pixelSizeIndex);
    170171  char* sphere_type=map.TypeOfMap();
     172      int premiereTranche = 0;
     173      int derniereTranche = map.NbThetaSlices()-1;
    171174  if (strncmp(sphere_type,"RING",4) == 0)
    172     {
    173       nsmax=map.SizeIndex();
    174     }
     175  {
     176    nsmax=map.SizeIndex();
     177  }
    175178  else
     179  {
    176180    // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
    177181    // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
    … …  
    179183    // c'est approximatif ; a raffiner.
    180184    if (strncmp(sphere_type,"TETAFI",6) == 0)
    181       {
    182         nsmax=(int_4)sqrt(map.NbPixels()/12.);
    183       }
     185    {
     186      nsmax=(int_4)sqrt(map.NbPixels()/12.);
     187      premiereTranche++;
     188      derniereTranche--;
     189    }
    184190    else
    185       {
    186         cout << " unknown type of sphere : " << sphere_type << endl;
    187         throw IOExc(" unknown type of sphere: " + (string)sphere_type );
    188       }
     191    {
     192      cout << " unknown type of sphere : " << sphere_type << endl;
     193      throw IOExc(" unknown type of sphere: " + (string)sphere_type );
     194    }
    189195  //  cout << "GenerateFromAlm: the sphere is of type : " << sphere_type << endl;
    190196  //  cout << "GenerateFromAlm: size index (nside) of the sphere= " << nsmax << endl;
    191197  //  cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
    192   if (nlmax>3*nsmax-1)
    193     {
     198  //  if (nlmax>3*nsmax-1)
     199  //  {
    194200      //     cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
    195       if (strncmp(sphere_type,"TETAFI",6) == 0)
    196         {
    197           cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
    198           cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
    199         }
    200     }
     201  //   if (strncmp(sphere_type,"TETAFI",6) == 0)
     202  //    {
     203  //      cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
     204  //      cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
     205  //    }
     206  //}
     207  }
    201208  Bm<complex<T> > b_m_theta(nmmax);
    202209
    … …  
    205212
    206213  // pour chaque tranche en theta
    207   for (int_4 ith = 0; ith < map.NbThetaSlices();ith++)
     214  for (int_4 ith = premiereTranche; ith <= derniereTranche;ith++)
    208215    {
    209216      int_4 nph;
    … …  
    231238            }
    232239        }
    233 
    234240    //        obtains the negative m of b(m,theta) (= complex conjugate)
    235241
    … …  
    377383  // si phi0 n'est pas nul, il y a juste un decalage a faire.
    378384  //**********************************************************************
    379 
    380385  TVector< complex<T> > bw(nph);
    381386  TVector< complex<T> > dataout(nph);
    … …  
    410415
    411416  TVector<T> sortie;
    412 
    413417  fftIntfPtr_-> FFTBackward(data, sortie);
    414418 
    … …  
    419423 /*! \fn  Alm<T> SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
    420424
    421 \return the Alm coefficients from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0)
     425\return the Alm coefficients from analysis of a temperature map.
    422426
    423427    \param<nlmax> : maximum value of the l index
    … …  
    425429     \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
    426430
    427 \param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis)
    428   */
     431 */
     432template<class T>
     433void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut) const
     434{
     435  DecomposeToAlm(const_cast< SphericalMap<T>& >(map), alm, nlmax, cos_theta_cut, 0);
     436}
     437//*******************************************
     438
     439 /*! \fn  Alm<T> SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
     440
     441\return the Alm coefficients from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0)
     442
     443    \param<nlmax> : maximum value of the l index
     444
     445     \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
     446
     447\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis). If iterationOrder is not null, the method works with SphereHEALPix but NOT WITH SphereThetaPhi maps !  */
    429448template<class T>
    430449void SphericalTransformServer<T>::DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
    … …  
    606625      throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
    607626    }
    608   if (nph%2 != 0 )
    609     {
    610       throw PException("SphericalTransformServer<T>::CFromFourierAnalysis : longueur de datain impair ?");
    611     }
     627  // if (nph%2 != 0 )
     628  //  {
     629      //  throw PException("SphericalTransformServer<T>::CFromFourierAnalysis : longueur de datain impair ?");
     630  //  }
    612631  TVector<complex<T> > transformedData;
    613632
    614633  // la taille du vecteur complexe retourne est nph/2+1 (si la taille
    615634  // du vecteur reel entre est nph)
     635  //   cout << " longueur de datain  = " << nph << endl;
    616636  fftIntfPtr_-> FFTForward(datain, transformedData);
    617 
     637  //  cout <<  " taille de la transformee " << transformedData.Size() << endl;
    618638  //  TVector< complex<T> > dataout(nmmax+1);
    619639  dataout.ReSize(nmmax+1);
    … …  
    732752    // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
    733753}
    734 
    735 
    736754 /*! \fn void SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& mapq,
    737755                                              const SphericalMap<T>& mapu,
    … …  
    749767
    750768 \c cos_theta_cut : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
     769
     770
     771 */
     772template<class T>
     773void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& mapq,
     774                                              const SphericalMap<T>& mapu,
     775                                              Alm<T>& alme,
     776                                              Alm<T>& almb,
     777                                              int_4 nlmax,
     778                                              r_8 cos_theta_cut) const
     779{
     780  DecomposeToAlm(const_cast< SphericalMap<T>& >(mapq), const_cast< SphericalMap<T>& >(mapu), alme, almb, nlmax, cos_theta_cut);
     781}
     782
     783 /*! \fn void SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& mapq,
     784                                              const SphericalMap<T>& mapu,
     785                                              Alm<T>& alme,
     786                                              Alm<T>& almb,
     787                                              int_4 nlmax,
     788                                              r_8 cos_theta_cut,
     789                                              int iterationOrder) const
     790
     791analysis of a polarization map into Alm coefficients.
     792
     793 The spheres \c mapq and \c mapu contain respectively the Stokes parameters.
     794
     795 \c a2lme and \c a2lmb will receive respectively electric and magnetic Alm's
     796    nlmax : maximum value of the l index
     797
     798 \c cos_theta_cut : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
     799
     800\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis). If iterationOrder is not null, the method works with SphereHEALPix but NOT WITH SphereThetaPhi maps !
     801
     802THE INPUT MAPS CAN BE MODIFIED (only if iterationOrder >0)
     803
    751804 */
    752805template<class T>
    … …  
    12431296
    12441297
    1245 /*! \fn TVector<T>  SOPHYA::SphericalTransformServer::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
     1298/*! \fn TVector<T>  SOPHYA::SphericalTransformServer::DecomposeToCl(SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
    12461299
    12471300\return power spectrum from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0)
    … …  
    12511304     \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
    12521305
    1253 \param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis)
     1306\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. If iterationOrder is not null, the method works with SphereHEALPix but NOT WITH SphereThetaPhi maps !
    12541307
    12551308  */
    … …  
    12591312  Alm<T> alm;
    12601313  DecomposeToAlm( sph, alm, nlmax, cos_theta_cut, iterationOrder);
    1261     //    cout << " SphericalTransformServer : impression alm " << endl;
    1262     //     alm.Print();
     1314  // power spectrum
     1315     return  alm.powerSpectrum();
     1316}
     1317
     1318
     1319/*! \fn TVector<T>  SOPHYA::SphericalTransformServer::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
     1320
     1321\return power spectrum from analysis of a temperature map.
     1322
     1323     \param<nlmax> : maximum value of the l index
     1324
     1325     \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
     1326
     1327
     1328  */
     1329
     1330
     1331template <class T>
     1332TVector<T>  SphericalTransformServer<T>::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
     1333{
     1334  Alm<T> alm;
     1335  DecomposeToAlm( sph, alm, nlmax, cos_theta_cut);
    12631336  // power spectrum
    12641337     return  alm.powerSpectrum();
Note: See TracChangeset for help on using the changeset viewer.