// Pour boucher les trous d'une sphere HEALPIX utilisant une autre sphere
//             cmv 13/6/01
#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 -0 nulval]"<<endl
    <<"           sphin.fits sphred.fits sphout.fits"<<endl
    <<" -w sphinw.fits : input sphere  weight (def=NULL)"<<endl
    <<" -0 nulval : value to identify empty pixel in sphin (def=0.)"<<endl
    <<"             (used 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[])
{
double nulval=0.;
char * fsphinw = NULL;
int c;
while((c = getopt(narg,arg,"h0:w:")) != -1) {
  switch (c) {
  case '0' :
    sscanf(optarg,"%lf",&nulval);
    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<<"Sphere Input : "<<fsphin<<endl
    <<"Weight Sphere Input : "<<fsphinw<<endl
    <<"Reducted Sphere Input : "<<fsphred<<endl
    <<"Sphere Output : "<<fsphout<<endl
    <<"nulval : "<<nulval<<endl;

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

// Creation de la sphere output
cout<<"Creating Sphere for output"<<endl;
SphereHEALPix<r_8> sphout(sphin,false);

// Filling hole for Output Sphere
cout<<"...Filling hole for Output Sphere"<<endl;
uint_4 n=0;
for(int_4 i=0;i<sphout.NbPixels();i++) {
  if(fsphinw) {if(sphinw(i)>0.) continue;}
     else     {if(sphout(i)!=nulval) continue;}
  double theta,phi;
  sphout.PixThetaPhi(i,theta,phi);
  int_4 ir = sphred.PixIndexSph(theta,phi);
  sphout(i) = sphred(ir);
  n++;
}
cout<<"      Number of filled pixels "<<n<<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<<sphout;
}

exit(0);
}