Changeset 1756 in Sophya for trunk/SophyaLib/Samba


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

adaptation aus spheres theta-phi

Location:
trunk/SophyaLib/Samba
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/Samba/alm.h

    r1683 r1756  
    1919    public :
    2020 
    21 Alm() : TriangularMatrix<complex<T> >()  {};
     21Alm() : TriangularMatrix<complex<T> >()  {;};
    2222    /* instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s */
    2323Alm(int nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;}

  • trunk/SophyaLib/Samba/lambdaBuilder.cc

    r1683 r1756  
    4040}
    4141
    42 TriangularMatrix<r_8>* LambdaLMBuilder::a_recurrence_ = NULL;
    43 TriangularMatrix<r_8>* LambdaLMBuilder::lam_fact_     = NULL;
     42TriangularMatrix<r_8> LambdaLMBuilder::a_recurrence_ = TriangularMatrix<r_8>();
     43TriangularMatrix<r_8> LambdaLMBuilder::lam_fact_     = TriangularMatrix<r_8>();
    4444TVector<r_8>* LambdaLMBuilder::normal_l_     = NULL;
    4545
    … …  
    8181void LambdaLMBuilder::array_init(int lmax, int mmax)
    8282   {
    83      if (a_recurrence_ == NULL)
    84        {
    85          a_recurrence_ = new TriangularMatrix<r_8>;
     83     if (a_recurrence_.Size() == 0)
     84       {
     85         //      a_recurrence_ = new TriangularMatrix<r_8>;
    8686         updateArrayRecurrence(lmax);
    8787       }
    8888     else
    89        if ( lmax > (*a_recurrence_).rowNumber()-1     )
     89       if ( lmax > (a_recurrence_).rowNumber()-1     )
    9090         {
    9191           cout << " WARNING : The classes LambdaXXBuilder will be more efficient if instanciated with parameter lmax = maximum value of l index which will be needed in the whole application (arrays not recomputed) " << endl;
    92            cout << "lmax= " << lmax << " previous instanciation with lmax= " <<  (*a_recurrence_).rowNumber() << endl;
     92           cout << "lmax= " << lmax << " previous instanciation with lmax= " <<  (a_recurrence_).rowNumber() << endl;
    9393         updateArrayRecurrence(lmax);
    9494         }
    … …  
    110110         r_8 lam_1=1. /bignorm2 ;
    111111         //      r_8 a_rec = LWK->a_recurr(m,m);
    112          r_8 a_rec = (*a_recurrence_)(m,m);
     112         r_8 a_rec = a_recurrence_(m,m);
    113113         r_8 b_rec = 0.;
    114114         for (int l=m+1; l<=lmax_; l++)
    … …  
    118118           b_rec=1./a_rec;
    119119           //      a_rec= LWK->a_recurr(l,m);
    120            a_rec= (*a_recurrence_)(l,m);
     120           a_rec= a_recurrence_(l,m);
    121121           lam_0 = lam_1;
    122122           lam_1 = lam_2;
    … …  
    135135void  LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)
    136136   {
    137      (*a_recurrence_).ReSizeRow(lmax+1);
     137     a_recurrence_.ReSizeRow(lmax+1);
    138138     for (int m=0; m<=lmax;m++)
    139139     {
    140        (*a_recurrence_)(m,m) = sqrt( 2.*m +3.);
     140       a_recurrence_(m,m) = sqrt( 2.*m +3.);
    141141       for (int l=m+1; l<=lmax; l++)
    142142         {
    143143           r_8 fl2 = (l+1.)*(l+1.);
    144            (*a_recurrence_)(l,m)=sqrt( (4.*fl2-1.)/(fl2-m*m) );
     144           a_recurrence_(l,m)=sqrt( (4.*fl2-1.)/(fl2-m*m) );
    145145         }
    146146     }
    … …  
    153153void  LambdaLMBuilder::updateArrayLamNorm()
    154154     {
    155        (*lam_fact_).ReSizeRow(lmax_+1);
     155       lam_fact_.ReSizeRow(lmax_+1);
    156156       for(int m = 0;m<= lmax_; m++)
    157157         {
    158158           for (int l=m; l<=lmax_; l++)
    159159             {
    160                (*lam_fact_)(l,m) =2.*(r_8)sqrt( (2.*l+1)*(l+m)*(l-m)/(2.*l-1)  );
     160               lam_fact_(l,m) =2.*(r_8)sqrt( (2.*l+1)*(l+m)*(l-m)/(2.*l-1)  );
    161161             }
    162162         }
    … …  
    214214void LambdaWXBuilder::array_init()
    215215   {
    216      if (lam_fact_ == NULL || normal_l_ == NULL)
    217        {
    218          lam_fact_ = new  TriangularMatrix<r_8>;
     216     if (lam_fact_.Size() < 1 || normal_l_ == NULL)
     217       {
     218         //      lam_fact_ = new  TriangularMatrix<r_8>;
    219219         normal_l_ = new TVector<r_8>;
    220220         updateArrayLamNorm();
    221221       }
    222222     else
    223        if ( lmax_ > (*lam_fact_).rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
     223       if ( lmax_ > lam_fact_.rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
    224224         {
    225225           updateArrayLamNorm();
    … …  
    247247          r_8 lam_lm1m = LambdaLMBuilder::lamlm(l-1,m);
    248248          r_8 lam_lm = LambdaLMBuilder::lamlm(l,m);
    249           r_8 lam_fact_l_m = (*lam_fact_)(l,m);
     249          r_8 lam_fact_l_m = lam_fact_(l,m);
    250250          r_8  a_w =  2. * (l - m*m) * one_on_s2 + l*(l-1.);
    251251          r_8  b_w =  c_on_s2 * lam_fact_l_m;
    … …  
    284284void LambdaPMBuilder::array_init()
    285285   {
    286      if (lam_fact_ == NULL || normal_l_ == NULL)
    287        {
    288          lam_fact_ = new  TriangularMatrix<r_8>;
     286     if (lam_fact_.Size() < 1 || normal_l_ == NULL)
     287       {
     288         //      lam_fact_ = new  TriangularMatrix<r_8>;
    289289         normal_l_ = new TVector<r_8>;
    290290         updateArrayLamNorm();
    291291       }
    292292     else
    293        if ( lmax_ > (*lam_fact_).rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
     293       if ( lmax_ > lam_fact_.rowNumber()-1  || lmax_ >  (*normal_l_).NElts()-1 )
    294294         {
    295295           updateArrayLamNorm();
    … …  
    318318             r_8 lam_lm1m = LambdaLMBuilder::lamlm(l-1,m);
    319319             r_8 lam_lm = LambdaLMBuilder::lamlm(l,m);
    320              r_8 lam_fact_l_m = (*lam_fact_)(l,m);
     320             r_8 lam_fact_l_m = lam_fact_(l,m);
    321321             r_8  a_w =  2. * (l - m*m) * one_on_s2 + l*(l-1.);
    322322             r_8 f_w =  lam_fact_l_m/(sth_*sth_);

  • trunk/SophyaLib/Samba/lambdaBuilder.h

    r1683 r1756  
    6060
    6161
    62  static TriangularMatrix<r_8>* a_recurrence_;
     62 // static TriangularMatrix<r_8>* a_recurrence_;
     63 static TriangularMatrix<r_8> a_recurrence_;
    6364 TriangularMatrix<r_8> lambda_;
    6465
    … …  
    6768 void updateArrayLamNorm();
    6869
    69  static  TriangularMatrix<r_8>* lam_fact_;
     70 // static  TriangularMatrix<r_8>* lam_fact_;
     71 static  TriangularMatrix<r_8> lam_fact_;
    7072 static  TVector<r_8>*  normal_l_;
    7173 int_4 lmax_;

  • 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();

  • trunk/SophyaLib/Samba/sphericaltransformserver.h

    r1683 r1756  
    4646
    4747
    48 void DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder = 0) const;
     48void DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const;
     49void DecomposeToAlm(const SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut) const;
    4950
    5051void DecomposeToAlm(SphericalMap<T>& mapq,
    … …  
    5455                    int_4 nlmax,
    5556                    r_8 cos_theta_cut,
    56                     int iterationOrder = 0) const;
     57                    int iterationOrder) const;
     58void DecomposeToAlm(const SphericalMap<T>& mapq,
     59                    const SphericalMap<T>& mapu,
     60                    Alm<T>& alme,
     61                    Alm<T>& almb,
     62                    int_4 nlmax,
     63                    r_8 cos_theta_cut) const;
    5764
    5865TVector<T>  DecomposeToCl(SphericalMap<T>& sph, 
    59                            int_4 nlmax, r_8 cos_theta_cut, int iterationOrder = 0) const;
     66                           int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const;
     67
     68TVector<T>  DecomposeToCl(const SphericalMap<T>& sph, 
     69                           int_4 nlmax, r_8 cos_theta_cut) const;
    6070
    6171
Note: See TracChangeset for help on using the changeset viewer.