Changeset 606 in Sophya for trunk/SophyaProg/Tests/tspm.cc
- Timestamp:
- Nov 20, 1999, 9:47:42 PM (26 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaProg/Tests/tspm.cc
r575 r606 2 2 #include "sambainit.h" 3 3 #include "spheregorski.h" 4 #include "spherethetaphi.h" 4 5 5 6 #include "tod.h" 6 7 #include "timing.h" 7 8 9 template <class T> 10 void MeanSig(PixelMap<T> const & map, double& gmoy, double& gsig) 11 { 12 gmoy=0.; 13 gsig = 0.; 14 double valok; 15 for(int k=0; k<map.NbPixels(); k++) { 16 valok = map(k); 17 gmoy += valok; gsig += valok*valok; 18 } 19 gmoy /= (double)map.NbPixels(); 20 gsig = gsig/(double)map.NbPixels() - gmoy*gmoy; 21 if (gsig >= 0.) gsig = sqrt(gsig); 22 } 8 23 9 24 int main(int narg, char* arg[]) 10 25 { 11 26 double teta,phi; 27 double gmoy, gsig; 28 12 29 PeidaInit(); 13 30 InitTim(); // Initializing the CPU timer 14 31 if ((narg > 1) && (strcmp(arg[1],"-h") == 0) ) { 15 cout << " tspm [Gorski_M=32] : GorskiSpherical Map Test " << endl;32 cout << " tspm [Gorski_M=32] [M_TetaPhi=64] : Spherical Map Test " << endl; 16 33 exit(0); 17 34 } … … 32 49 33 50 // Computing mean and sigma on the sphere 34 double gmoy=0. , gsig = 0.; 35 double valok; 36 for(int k=0; k<sph.NbPixels(); k++) { 37 valok = sph(k); 38 gmoy += valok; gsig += valok*valok; 39 } 40 cout << "SphMap Mean= " << gmoy << " Sigma = " << gsig << endl; 51 MeanSig(sph, gmoy, gsig); 52 cout << "SphereGorski<double> Mean= " << gmoy << " Sigma = " << gsig << endl; 41 53 PrtTim("End of Mean-Sig "); 42 54 … … 44 56 { 45 57 POutPersist s("sphg.ppf"); 46 FIO_SphereGorski<double> fiog( sph) ;58 FIO_SphereGorski<double> fiog(&sph) ; 47 59 fiog.Write(s); 48 cout << "Sph Mapwritten to sphg.ppf " << endl;60 cout << "SphereGorski<double> written to sphg.ppf " << endl; 49 61 } 50 62 … … 52 64 { 53 65 FIO_SphereGorski<double> fiog("sphg.ppf"); 54 double gmoy=0. , gsig = 0.;55 double valok;56 66 SphereGorski<double> sph2 = fiog; 57 67 PrtTim("End of Write/Read "); 58 68 69 cout << " Spherical map from file sphg.ppf NPixels= " << sph2.NbPixels() << endl; 70 59 71 int ndiff = 0; 60 72 for(int k=0; k<sph2.NbPixels(); k++) { 61 valok = sph2(k); 62 gmoy += valok; gsig += valok*valok; 63 if ((sph2(k)-sph(k)) > 1.e-49) ndiff++; 73 // if ((sph2(k)-sph(k)) > 1.e-49) { 74 if ( sph2(k) != sph(k) ) { 75 ndiff++; 76 if (ndiff < 20) cout << "!!!Diff: K= " << k << " SPH2= " << sph2(k) << " SPH= " << sph(k) << endl; 77 } 64 78 } 79 MeanSig(sph, gmoy, gsig); 80 cout << "SphMapFromFile Mean= " << gmoy << " Sigma = " << gsig << endl; 81 cout << " NDiff = " << ndiff << " (should be zero = 0) " << endl; 82 PrtTim("End of Mean-Sig "); 83 } 84 85 86 int mt=64; 87 if (narg > 2) mt = atoi(arg[2]); 88 cout << "\n ===== ThetaPhi Spherical Map Test MT= " << mt << endl; 89 90 SphereThetaPhi<float> spht(m); 91 92 cout << "Filling spherical map NPixels= " << spht.NbPixels() << endl; 93 for (int j=0;j<spht.NbPixels();j++) 94 { 95 spht.PixThetaPhi(j,teta,phi); 96 spht(j)= 0.5* cos(5.*teta)*sin(10.*phi); 97 } 98 PrtTim("End of Fill "); 99 100 // Computing mean and sigma on the sphere 101 MeanSig(spht, gmoy, gsig); 102 cout << "SphereThetaPhi<float> Mean= " << gmoy << " Sigma = " << gsig << endl; 103 PrtTim("End of Mean-Sig "); 104 105 // Writing to a PPF file 106 { 107 POutPersist s("spht.ppf"); 108 FIO_SphereThetaPhi<float> fiog(spht) ; 109 fiog.Write(s); 110 cout << "SphereThetaPhi<float> written to spht.ppf " << endl; 111 } 112 113 // Reading from the file 114 { 115 FIO_SphereThetaPhi<float> fiog("spht.ppf"); 116 SphereThetaPhi<float> sph2 = fiog; 117 PrtTim("End of Write/Read "); 118 119 cout << " Spherical map from file sph.ppf NPixels= " << sph2.NbPixels() << endl; 120 121 int ndiff = 0; 122 for(int k=0; k<sph2.NbPixels(); k++) { 123 // if ((sph2(k)-sph(k)) > 1.e-49) { 124 if ( sph2(k) != spht(k) ) { 125 ndiff++; 126 if (ndiff < 20) cout << "!!!Diff: K= " << k << " SPH2= " << sph2(k) << " SPH= " << sph(k) << endl; 127 } 128 } 129 MeanSig(sph, gmoy, gsig); 65 130 cout << "SphMapFromFile Mean= " << gmoy << " Sigma = " << gsig << endl; 66 131 cout << " NDiff = " << ndiff << " (should be zero = 0) " << endl;
Note:
See TracChangeset
for help on using the changeset viewer.