// Creation d'une sphere de masque a partir d'une sphere de valeurs
//         cmv 13/6/01
// cremskfrsph -m 0.1 -v 1.,0. sph143k05_e.fits sphmskw.fits
// cremskfrsph -n -m -1.e-30. -M 1.e-30. -v 1.,0. sph143k05.fits sphmsk.fits
#include "machdefs.h"
#include <unistd.h>
#include <stdexcept>
#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include "skymapinit.h"
#include "skymap.h"
#include "fitsspherehealpix.h" 

/*!
  \ingroup PrgMap
  \file cremskfrsph.cc
  \brief \b cremskfrsph: Create a masked sphere from a sphere of datas
  \verbatim
csh> cremskfrsph -h
cremskfrsph [-n -m min -M max -v valmsk,inimsk] sphval.fits sphmask.fits
 -m min : min value
 -M max : max value
    condition for filling masque is: min<=sphere_value<=max
          (fill mask with valmsk)     (bounds are included)
 -n : negate the condition
      (switch to: sphere_value<min || max<sphere_value)
 -b bmin,bmax : on coupe en latitude entre [bmin,bmax] degres (def=no+cut)
 -v valmsk,inimsk:
    valmsk: value to be put into mask if ok (def=1.)
    inimsk: initialisation value for the mask (def=0.)
  \endverbatim
*/

void usage();
void usage()
{
cout<<"cremskfrsph [-n -m min -M max -v valmsk,inimsk] sphval.fits sphmask.fits"<<endl
    <<" -m min : min value"<<endl
    <<" -M max : max value"<<endl
    <<"    condition for filling masque is: min<=sphere_value<=max"<<endl
    <<"          (fill mask with valmsk)     (bounds are included)"<<endl
    <<" -n : negate the condition"<<endl
    <<"      (switch to: sphere_value<min || max<sphere_value)"<<endl
    <<" -b bmin,bmax : on coupe en latitude entre [bmin,bmax] degres (def=no+cut)"<<endl
    <<" -v valmsk,inimsk:"<<endl
    <<"    valmsk: value to be put into mask if ok (def=1.)"<<endl
    <<"    inimsk: initialisation value for the mask (def=0.)"<<endl;
}


int main(int narg, char* arg[])
{
double vmin=-1.e30,vmax=1.e30,vmask=1.,vmaskini=0., bmin=1.,bmax=-1.;
bool tstmin=false,tstmax=false,negate=false;
string dum;
int c;
while((c = getopt(narg,arg,"hnm:M:v:b:")) != -1) {
  switch (c) {
  case 'n' :
    negate = true;
    break;
  case 'm' :
    sscanf(optarg,"%lf",&vmin);
    tstmin = true;
    break;
  case 'M' :
    sscanf(optarg,"%lf",&vmax);
    tstmax = true;
    break;
  case 'v' :
    sscanf(optarg,"%lf,%lf",&vmask,&vmaskini);
    break;
  case 'b' :
    sscanf(optarg,"%lf,%lf",&bmin,&bmax);
    break;
  case 'h' :
  default:
    usage(); exit(1);
  }
}

if(optind+1>=narg) {usage(); exit(1);}
char * sphval = arg[optind];
char * sphmsk = arg[optind+1];

if(negate) dum = ".NOT."; else dum = "";
cout<<"Sphere values : "<<sphval<<endl
    <<"Sphere mask   : "<<sphmsk<<endl
    <<"  ...min("<<tstmin<<") "<<vmin<<endl
    <<"     max("<<tstmax<<") "<<vmax<<endl
    <<"  ...negate "<<negate<<endl
    <<"  ...mask set value "<<vmask<<endl
    <<"     mask init value "<<vmaskini<<endl
    <<"Cut latitude between "<<bmin<<" and "<<bmax<<" deg"<<endl
    <<"Condition on data to set the mask: "
    <<dum<<"( "<<vmin<<" <= V <= "<<vmax<<" )"<<endl;

// Lecture de la sphere Healpix des valeurs
SphereHEALPix<r_8> sph;
FitsInFile sfits(sphval);
sfits >> sph;
cout<<"Opening Sphere HEALPix for testing values :"<<endl
    <<"          Type of map : "<<sph.TypeOfMap()<<endl
    <<"     Number of pixels : "<<sph.NbPixels()<<endl
    <<"                 Nlat : "<<sph.SizeIndex()<<endl;

// Ouverture de la sphere Healpix pour le masque
cout<<"Creating Mask Sphere"<<endl;
SphereHEALPix<r_8> sphm(sph,false);

// Filling Mask Sphere with cut on values
cout<<"Filling Mask"<<endl;
double thetamin = (90.-bmax) * M_PI/180.;
double thetamax = (90.-bmin) * M_PI/180.;
uint_4 nmask = 0;
for(int_4 i=0;i<sph.NbPixels();i++) {
  sphm(i) = vmaskini;

  bool skp = (tstmin || tstmax) ? negate : false;
  if((tstmin && sph(i)<vmin) || (tstmax && sph(i)>vmax)) skp = !negate;
  if(skp) continue;   // Do nothing

  if(bmin<bmax) {
    double theta,phi;
    sphm.PixThetaPhi(i,theta,phi);
    if(theta>=thetamin && theta<=thetamax) continue;
  }

  sphm(i) = vmask;
  nmask++;
}
cout<<"    .... Number of values set in mask   : "<<nmask<<endl;
cout<<"    .... Fraction of values set in mask : "
    <<(r_8)nmask/sph.NbPixels()*100.<<" %"<<endl;

// Ecriture de la sphere Healpix sur fits
{
dum = "rm -f "; dum += sphmsk; system(dum.c_str());
FitsOutFile swfits(sphmsk,FitsFile::clear);
cout<<"Writing Mask Sphere Fits file"<<endl;
swfits<<sphm;
}

return 0;
}