Changeset 606 in Sophya for trunk/SophyaProg/Tests

Timestamp:

Nov 20, 1999, 9:47:42 PM (26 years ago)

Author:

ansari

Message:

Amelioration test spheres - Reza 20/11/99

Location:

trunk/SophyaProg/Tests

Files:

• tspm.cc (modified) (4 diffs)
• tspm2.cc (modified) (5 diffs)

trunk/SophyaProg/Tests/tspm.cc

@@ -2 +2 @@
 #include "sambainit.h"
 #include "spheregorski.h"
+#include "spherethetaphi.h"
 
 #include "tod.h"
 #include "timing.h"
 
+template <class T>
+void  MeanSig(PixelMap<T> const & map, double& gmoy, double& gsig)
+{
+  gmoy=0.;
+  gsig = 0.;
+  double valok;
+  for(int k=0; k<map.NbPixels(); k++) {
+    valok = map(k);
+    gmoy += valok;  gsig += valok*valok;
+  }
+  gmoy /= (double)map.NbPixels();
+  gsig = gsig/(double)map.NbPixels() - gmoy*gmoy;
+  if (gsig >= 0.) gsig = sqrt(gsig);
+}
 
 int main(int narg, char* arg[]) 
 {
   double teta,phi;
+  double gmoy, gsig;
+
   PeidaInit();
   InitTim();   // Initializing the CPU timer 
   if ((narg > 1) && (strcmp(arg[1],"-h") == 0) )  {
-    cout << " tspm [Gorski_M=32] : Gorski Spherical Map Test " << endl;
+    cout << " tspm [Gorski_M=32] [M_TetaPhi=64] : Spherical Map Test " << endl;
     exit(0);
     }
@@ -32 +49 @@
 
   // Computing mean and sigma on the sphere 
-  double gmoy=0. , gsig = 0.; 
-  double valok;
-  for(int k=0; k<sph.NbPixels(); k++) {
-    valok = sph(k);
-    gmoy += valok;  gsig += valok*valok; 
-  }
-  cout << "SphMap Mean= " << gmoy << "  Sigma = " << gsig << endl;
+  MeanSig(sph, gmoy, gsig);
+  cout << "SphereGorski<double> Mean= " << gmoy << "  Sigma = " << gsig << endl;
   PrtTim("End of Mean-Sig ");
 
@@ -44 +56 @@
   {
   POutPersist s("sphg.ppf"); 
-  FIO_SphereGorski<double> fiog(sph) ;
+  FIO_SphereGorski<double> fiog(&sph) ;
   fiog.Write(s);
-  cout << "SphMap written to sphg.ppf "  << endl;
+  cout << "SphereGorski<double>  written to sphg.ppf "  << endl;
   }
 
@@ -52 +64 @@
   {
   FIO_SphereGorski<double> fiog("sphg.ppf");
-  double gmoy=0. , gsig = 0.; 
-  double valok;
   SphereGorski<double>  sph2 = fiog;
   PrtTim("End of Write/Read ");
 
+  cout << " Spherical map from file sphg.ppf NPixels= " << sph2.NbPixels() << endl;
+
   int ndiff = 0;
   for(int k=0; k<sph2.NbPixels(); k++) {
-    valok = sph2(k);
-    gmoy += valok;  gsig += valok*valok; 
-    if ((sph2(k)-sph(k)) > 1.e-49)  ndiff++;
+//    if ((sph2(k)-sph(k)) > 1.e-49)  { 
+    if ( sph2(k) != sph(k) )   { 
+      ndiff++;
+      if (ndiff < 20)  cout << "!!!Diff: K= " << k << " SPH2= " << sph2(k) << " SPH= " << sph(k) << endl;
+    }
   }
+  MeanSig(sph, gmoy, gsig);
+  cout << "SphMapFromFile Mean= " << gmoy << "  Sigma = " << gsig << endl;
+  cout << " NDiff = " << ndiff << " (should be zero = 0) " << endl;
+  PrtTim("End of Mean-Sig ");
+  }
+
+
+  int  mt=64;
+  if (narg > 2)  mt = atoi(arg[2]);
+  cout <<  "\n ===== ThetaPhi Spherical Map Test MT= " << mt << endl;
+
+  SphereThetaPhi<float> spht(m);
+
+  cout << "Filling spherical map NPixels= " << spht.NbPixels() << endl;
+  for (int j=0;j<spht.NbPixels();j++) 
+    {
+      spht.PixThetaPhi(j,teta,phi);
+      spht(j)= 0.5* cos(5.*teta)*sin(10.*phi);
+    }
+  PrtTim("End of Fill ");
+
+  // Computing mean and sigma on the sphere 
+  MeanSig(spht, gmoy, gsig);
+  cout << "SphereThetaPhi<float> Mean= " << gmoy << "  Sigma = " << gsig << endl;
+  PrtTim("End of Mean-Sig ");
+
+//  Writing to a PPF file
+  {
+  POutPersist s("spht.ppf"); 
+  FIO_SphereThetaPhi<float> fiog(spht) ;
+  fiog.Write(s);
+  cout << "SphereThetaPhi<float> written to spht.ppf "  << endl;
+  }
+
+  // Reading from the file
+  {
+  FIO_SphereThetaPhi<float> fiog("spht.ppf");
+  SphereThetaPhi<float>  sph2 = fiog;
+  PrtTim("End of Write/Read ");
+
+  cout << " Spherical map from file sph.ppf NPixels= " << sph2.NbPixels() << endl;
+
+  int ndiff = 0;
+  for(int k=0; k<sph2.NbPixels(); k++) {
+//    if ((sph2(k)-sph(k)) > 1.e-49)  { 
+    if ( sph2(k) != spht(k) )   { 
+      ndiff++;
+      if (ndiff < 20)  cout << "!!!Diff: K= " << k << " SPH2= " << sph2(k) << " SPH= " << sph(k) << endl;
+    }
+  }
+  MeanSig(sph, gmoy, gsig);
   cout << "SphMapFromFile Mean= " << gmoy << "  Sigma = " << gsig << endl;
   cout << " NDiff = " << ndiff << " (should be zero = 0) " << endl;

trunk/SophyaProg/Tests/tspm2.cc

@@ -58 +58 @@
   printf("Xsize= %d   Ysize= %d  NPix= %d\n",img.XSize(),img.YSize(),img.XSize()*img.YSize() );
   n = 0;
+  img.Zero();
   for(j=0; j<img.YSize(); j++) {
     yd = (r_4)(j+0.5)/(r_4)img.YSize()-0.5;
@@ -112 +113 @@
    nomobj = "sphg1";
    fiog.Write(s, nomobj);
-   cout << "SphMap written to POutPersist with name "  << nomobj << endl;
-   }
-
+   cout << "SphMap SphereGorski<double> written to POutPersist with name "  << nomobj << endl;
+   }
+
+  FitsIoServer fios;
   // On projete dans un fichier FITS
-  FitsIoServer fios;
   fios.sinus_picture_projection(sph, "gsin1.fits");
   fios.save(sph, "sph1.fits");
@@ -146 +147 @@
   PrtTim("End of Mean-Sig ");
   }
+
   {
-  SphereGorski<double> sph(m);
+  SphereGorski<float> sph(m);
 
   cout << "Filling spherical map2 NPixels= " << sph.NbPixels() << endl;
@@ -164 +166 @@
 
    {
-   FIO_SphereGorski<double> fiog(&sph) ;
+   FIO_SphereGorski<float> fiog(&sph) ;
    nomobj = "sphg2";
    fiog.Write(s, nomobj);
-   cout << "SphMap written to POutPersist with name "  << nomobj << endl;
+   cout << "SphMap SphereGorski<float> written to POutPersist with name "  << nomobj << endl;
    }
 
@@ -173 +175 @@
   // On projete dans un fichier FITS
   FitsIoServer fios;
+   {
+   cout << "Test of Write/Read SphereGorski<float> to FITS (sphg_r4.fits) " << endl;
+   fios.save(sph, "sphg_r4.fits");
+ 
+   SphereGorski<float> sphr(4);
+   fios.load(sphr, "sphg_r4.fits", 2);
+   cout << " Read from file - SphereGorski<float> NPixels= " << sphr.NbPixels() << endl;
+   int ndiff = 0;
+   for(int k=0; k<sphr.NbPixels(); k++) {
+    if ( sphr(k) != sph(k) )   { 
+      ndiff++;
+      if (ndiff < 20)  cout << "!!!Diff: K= " << k << " SPHR= " << sphr(k) << " SPH= " << sph(k) << endl;
+      }
+    }
+   cout << " ReadFrom FITS NDiff = " << ndiff << " (should be zero = 0) " << endl;
+   
+   }
+
+
   fios.sinus_picture_projection(sph, "gsin2.fits");
   fios.save(sph, "sph2.fits");
 
-  TMatrix<double> mtx(3*m, 6*m);
+  TMatrix<float> mtx(3*m, 6*m);
   Project_Mol(sph, mtx);
   fios.save(mtx, "mtx2.fits");