Index: trunk/SophyaProg/Tests/tspm.cc
===================================================================
--- trunk/SophyaProg/Tests/tspm.cc	(revision 575)
+++ trunk/SophyaProg/Tests/tspm.cc	(revision 606)
@@ -2,16 +2,33 @@
 #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,11 +49,6 @@
 
   // 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,7 +56,7 @@
   {
   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,15 +64,68 @@
   {
   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;
