Changeset 2313 in Sophya for trunk/SophyaLib/Samba

Timestamp:

Jan 22, 2003, 5:23:53 PM (23 years ago)

Author:

lemeur

Message:

transf. fourier ajustee pour pixelisations autres que healpix

Location:

trunk/SophyaLib/Samba

Files:

• sphericaltransformserver.cc (modified) (9 diffs)
• sphericaltransformserver.h (modified) (1 diff)

trunk/SophyaLib/Samba/sphericaltransformserver.cc

@@ -247 +247 @@
       //    sum_m  b(m,theta)*exp(i*m*phi)   -> f(phi,theta)
       // ---------------------------------------------------------------*/
-      TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
-      for (int i=0;i< nph;i++)
-        {
-          map(pixNumber(i))=Temp(i);
+
+
+      if (sphere_type.substr(0,4) == "RING")
+        {
+          TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
+          for (int i=0;i< nph;i++) map(pixNumber(i))=Temp(i);
+        }
+      else 
+        // pour des pixelisations quelconques (autres que HEALPix
+        //  nph n'est pas toujours pair
+        // ca fait des problemes pour les transformees de Fourier
+        // car le server de TF ajuste la longueur du vecteur reel 
+        // en sortie de TF, bref, la securite veut qu'on prenne une 
+        // TF complexe
+        {
+          TVector<complex<T> > Temp = fourierSynthesisFromB(b_m_theta,nph,phi0);
+          for (int i=0;i< nph;i++) map(pixNumber(i))=Temp(i).real();
         }
     }
@@ -377 +390 @@
   // somme_m' : de 0 a nph-1
   // On echange les sommations :
-  // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
+  // DT/T(j) = somme_m' (exp(i*m'*phi(j)) somme_k b(k*nph + m')*exp(i*(k*nph*phi(j))
   // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle 
   // vaut 1.
-  // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
+  // Il reste a calculer les transformees de Fourier de somme_k b(k*nph + m')
   // si phi0 n'est pas nul, il y a juste un decalage a faire.
   //**********************************************************************
@@ -407 +420 @@
 
   //     applies the shift in position <-> phase factor in Fourier space
+  // cout << " TF : nph= " << nph << " vec. entree " << data.Size() << endl;
   for (int mprime=0; mprime <= nph/2; mprime++)
     {
@@ -715 +729 @@
       
     }
+  bool healpix = true;
   if (sphere_type.substr(0,4) == "RING")
     {
@@ -724 +739 @@
     // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
     // c'est approximatif ; a raffiner.
+    healpix = false;
     if (sphere_type.substr(0,6) == "TETAFI")
       {
@@ -749 +765 @@
       throw SzMismatchError("SphericalTransformServer: arrays Alme and Almb have not the same size ?  ");
     }
-    mapFromWX(nlmax, nmmax, mapq, mapu, alme, almb);
+    mapFromWX(nlmax, nmmax, mapq, mapu, alme, almb, healpix);
     // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
 }
@@ -1059 +1075 @@
                                          SphericalMap<T>& mapu, 
                                          const Alm<T>& alme,
-                                         const Alm<T>& almb) const
+                                         const Alm<T>& almb, bool healpix) const
 
 synthesis of Stokes parameters following formulae : 
@@ -1086 +1102 @@
                                          SphericalMap<T>& mapu, 
                                          const Alm<T>& alme,
-                                         const Alm<T>& almb) const
-{
+                                         const Alm<T>& almb, bool healpix) const
+{
+  int i;
+
   Bm<complex<T> > b_m_theta_q(nmmax);
   Bm<complex<T> > b_m_theta_u(nmmax);
@@ -1137 +1155 @@
           b_m_theta_u(-m) = conj(b_m_theta_u(m));
         }
-
-      //    TVector<complex<T> > Tempq = fourierSynthesisFromB(b_m_theta_q,nph,phi0);  
-      //     TVector<complex<T> > Tempu = fourierSynthesisFromB(b_m_theta_u,nph,phi0); 
-      TVector<T> Tempq = RfourierSynthesisFromB(b_m_theta_q,nph,phi0);  
-      TVector<T> Tempu = RfourierSynthesisFromB(b_m_theta_u,nph,phi0); 
-      for (int i=0;i< nph;i++)
-        {
-          //      mapq(pixNumber(i))=Tempq(i).real();
-          //      mapu(pixNumber(i))=Tempu(i).real();
-          mapq(pixNumber(i))=Tempq(i);
-          mapu(pixNumber(i))=Tempu(i);
-          
+      if (healpix) 
+        {
+          TVector<T> Tempq = RfourierSynthesisFromB(b_m_theta_q,nph,phi0);  
+          TVector<T> Tempu = RfourierSynthesisFromB(b_m_theta_u,nph,phi0); 
+          for (i=0;i< nph;i++)
+            {
+              mapq(pixNumber(i))=Tempq(i);
+              mapu(pixNumber(i))=Tempu(i);
+            }
+        }
+      else
+        // pour des pixelisations quelconques (autres que HEALPix
+        //  nph n'est pas toujours pair
+        // ca fait des problemes pour les transformees de Fourier
+        // car le server de TF ajuste la longueur du vecteur reel 
+        // en sortie de TF, bref, la securite veut qu'on prenne une 
+        // TF complexe
+        {
+          TVector<complex<T> > Tempq = fourierSynthesisFromB(b_m_theta_q,nph,phi0);  
+          TVector<complex<T> > Tempu = fourierSynthesisFromB(b_m_theta_u,nph,phi0); 
+          for (i=0;i< nph;i++)
+            {
+              mapq(pixNumber(i))=Tempq(i).real();
+              mapu(pixNumber(i))=Tempu(i).real();
+            }
         }
     } 

trunk/SophyaLib/Samba/sphericaltransformserver.h

@@ -107 +107 @@
   void mapFromWX(int_4 nlmax, int_4 nmmax,
                SphericalMap<T>& mapq, SphericalMap<T>& mapu, 
-               const Alm<T>& alme, const Alm<T>& almb) const;
+               const Alm<T>& alme, const Alm<T>& almb, bool healpix) const;