// Pour boucher les trous d'une sphere HEALPIX utilisant une autre sphere
//             cmv 13/6/01
// extrap2sph -w sph143k05_e.fits sph143k05.fits sphred.fits sphout_2.fits
// extrap2sph -m -1.e-30 -M +1.e-30 sph143k05.fits sphred.fits sphout_22.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" 

void usage();
void usage()
{
cout<<"extrap2sph [-w sphinw.fits -m valmin -M valmax]"<<endl
    <<"           sphin.fits sphred.fits sphout.fits"<<endl
    <<" -w sphinw.fits : input sphere  weight (def=NULL)"<<endl
    <<" -m minval : value to identify EMPTY pixel (def=-1.e30.)"<<endl
    <<" -M maxval : value to identify EMPTY pixel (def=+1.e30)"<<endl
    <<"    condition for EMPTY pixel is (minval<val<maxval)"<<endl
    <<"    (bounds are excluded, used only if sphinw.fits==NULL)"<<endl
    <<" sphin.fits   : input sphere"<<endl
    <<" sphred.fits  : input reducted sphere"<<endl
    <<" sphout.fits  : output sphere"<<endl;
}

int main(int narg, char* arg[])
{
bool tstmin=false, tstmax=false;
double valmin=-1.e30, valmax=+1.e30;
char * fsphinw = NULL;
int c;
while((c = getopt(narg,arg,"hm:M:w:")) != -1) {
  switch (c) {
  case 'm' :
    sscanf(optarg,"%lf",&valmin);
    tstmin=true;
    break;
  case 'M' :
    sscanf(optarg,"%lf",&valmax);
    tstmax=true;
    break;
  case 'w' :
    fsphinw = optarg;
    break;
  case 'h' :
  default:
    usage(); exit(1);
  }
}

if(optind+2>=narg) {usage(); exit(1);}
char * fsphin   = arg[optind];
char * fsphred = arg[optind+1];
char * fsphout  = arg[optind+2];

cout<<"Input Sphere          : "<<fsphin<<endl
    <<"Input Weight Sphere   : "<<fsphinw<<endl
    <<"Input Reducted Sphere : "<<fsphred<<endl
    <<"Output Sphere         : "<<fsphout<<endl;
if(!fsphinw)
  cout<<"- test min("<<tstmin<<") : "<<valmin<<endl
      <<"  test max("<<tstmax<<") : "<<valmax<<endl
      <<"Condition for EMPTY pixel is :\n    ("<<valmin
      <<" < sphere value < "<<valmax<<") bounds are excluded!"<<endl;

// Lecture des spheres
SphereHEALPix<r_8> sphin;
{FitsInFile sfits(fsphin); sfits >> sphin;}
int nlat = sphin.SizeIndex();
cout<<"Opening Input Sphere :"<<endl
    <<"          Type of map : "<<sphin.TypeOfMap()<<endl
    <<"     Number of pixels : "<<sphin.NbPixels()<<endl
    <<"                 Nlat : "<<sphin.SizeIndex()<<endl;

SphereHEALPix<r_8> sphinw;
if(fsphinw) {
  FitsInFile sfits(fsphinw);
  sfits >> sphinw;
  cout<<"Opening Input Weight Sphere :"<<endl
      <<"          Type of map : "<<sphinw.TypeOfMap()<<endl
      <<"     Number of pixels : "<<sphinw.NbPixels()<<endl
      <<"                 Nlat : "<<sphinw.SizeIndex()<<endl;
  if(sphinw.SizeIndex()!=nlat) {
    cout<<"Incompatible sphin sphinw"<<endl;
    exit(2);
  }
}

SphereHEALPix<r_8> sphred;
{FitsInFile sfits(fsphred); sfits >> sphred;}
cout<<"Opening Input Reducted Sphere :"<<endl
    <<"          Type of map : "<<sphred.TypeOfMap()<<endl
    <<"     Number of pixels : "<<sphred.NbPixels()<<endl
    <<"                 Nlat : "<<sphred.SizeIndex()<<endl;

// Filling hole for Output Sphere
cout<<"...Filling hole for Output Sphere"<<endl;
uint_4 n=0, n0=0;
for(int_4 i=0;i<sphin.NbPixels();i++) {
  if(fsphinw) {
    if(sphinw(i)>0.) continue;
  } else {
    // Not Filled pixels ]valmin,valmax[
    // Filled pixels ]-Infinity,valmin] or [valmax,+Infinity[
    if(tstmin && sphin(i)<=valmin) continue;
    if(tstmax && sphin(i)>=valmax) continue;
  }
  double theta,phi;
  sphin.PixThetaPhi(i,theta,phi);
  int_4 ir = sphred.PixIndexSph(theta,phi);
  sphin(i) = sphred(ir);
  n++; if(sphred(ir)!=0.) n0++;
}
cout<<"      Number of pixels filled "
    <<n<<" (not null values "<<n0<<")"<<endl;

// Ecriture de la sphere
{
string dum = "rm -f "; dum += fsphout; system(dum.c_str());
FitsOutFile sfits(fsphout);
cout<<"Writing Output Sphere Fits file"<<endl;
sfits<<sphin;
}

exit(0);
}