source: Sophya/trunk/ArchTOIPipe/TestPipes/map_comp.cc@ 3366

//#define TOISEQBUFFERED

#include <unistd.h>
#include "toi.h"
#include "fitstoirdr.h"
#include "fitstoiwtr.h"
#include "fitsspherehealpix.h" 
#include "sambainit.h"
#include <stdexcept>

void usage(void);
void usage(void) {
 cout<<"map_comp [-h]  fitsphere1 fitshere2 fitsphere2w fitsout"<<endl;
 return;
}

////////////////////////////////////////////////////////////////
int main(int narg, char** arg) {


//-- Decodage arguments

if(optind+3>=narg) {usage(); exit(2);}
string const fitsphere1  = arg[optind];
string const fitsphere2  = arg[optind+1];
string const fitsphere2w = arg[optind+2];
string const fitsphout   = arg[optind+3];

cout<<">>>> map_comp:"<<endl
    <<"Fits OutFile "<<fitsphout<<endl
    <<"Fits Healpix Spheres "<<fitsphere1<< " " <<fitsphere2<< " " << fitsphere2w <<endl;

SophyaInit();

//--------------------------------------------------------------------
try {
//--------------------------------------------------------------------


 // Lecture de la sphere Healpix
 SphereHEALPix<r_8> sph1,sph2,sph2w,sphdiff;        
 FitsInFile sfits1(fitsphere1);
 sfits1 >> sph1;
 cout<<"SphereHEALPix: Type de map : "<<sph1.TypeOfMap()<<endl
     <<"               Nombre de pixels : "<<sph1.NbPixels()<<endl;
 FitsInFile sfits2(fitsphere2);
 sfits2 >> sph2;
 cout<<"SphereHEALPix: Type de map : "<<sph2.TypeOfMap()<<endl
     <<"               Nombre de pixels : "<<sph2.NbPixels()<<endl;

 if((sph1.TypeOfMap() != sph2.TypeOfMap()) || (sph1.NbPixels() !=sph2.NbPixels())) {
   cout << " Maps differ " << endl;
   exit(-1);
 }
 FitsInFile sfits2w(fitsphere2w);
 sfits2w >> sph2w;
 cout<<"SphereHEALPix: Type de map : "<<sph2w.TypeOfMap()<<endl
     <<"               Nombre de pixels : "<<sph2w.NbPixels()<<endl;
 if((sph2w.TypeOfMap() != sph2.TypeOfMap()) || (sph2w.NbPixels() !=sph2.NbPixels())) {
   cout << " Maps differ" << endl;
   exit(-1);
 }
 int nlat = sph1.SizeIndex();
 int npix = sph1.NbPixels();
 // Creation de la sphere Healpix
 SphereHEALPix<r_8>* sph = new SphereHEALPix<r_8>(nlat);
 cout<<"SphereHEALPix: Type de map : "<<sph->TypeOfMap()<<endl
     <<"               Nombre de pixels : "<<sph->NbPixels()<<endl
     <<"               Nlat : "<<sph->SizeIndex()<<endl;


 TArray<r_8> diff(npix);
 // for (int j=0;j<npix;j++) {cout << j << " " << sph1.PixVal(j)<<endl;}
 int nok=0;

 for (int i=0;i<npix;i++) if ((sph2w.PixVal(i) > 0) && (sph1.PixVal(i) > -999.)) {
   diff[i] = sph1.PixVal(i) - sph2.PixVal(i);
   sph->PixVal(i)=diff[i];
   nok ++;
 } else sph->PixVal(i)=-99999999.;

 r_8 min,max;
 diff.MinMax(min,max);
 cout << " MinMax = " << min << " " << max << endl;
 cout << " Mean =  " << diff.Sum()/npix << endl;
 cout << " Var = " << diff.SumX2()/npix << endl;
 cout << "Values computed on " << double(nok)/npix*100. << "% of the map" << endl;

 // Ecriture de la sphere Healpix sur fits
 {
 FitsOutFile sfits(fitsphout,FitsFile::clear);
 cout<<"tsttoi2map: Creating sphere fits file "<<fitsphout<<endl;
 sfits << *sph;
 }


//--------------------------------------------------------------------
} catch (PThrowable & exc) {
 cout<<"\ntstmap2toi: Catched Exception \n"<<(string)typeid(exc).name() 
     <<" - Msg= "<<exc.Msg()<<endl;
} catch (const std::exception & sex) {
 cout<<"\ntstmap2toi: Catched std::exception \n" 
     <<(string)typeid(sex).name()<<endl; 
} catch (...) {
 cout<<"\ntstmap2toi: some other exception was caught ! "<<endl;
}
//--------------------------------------------------------------------

exit(0);
}