// Pour boucher les trous d'une sphere HEALPIX en clusterisant dans une sphere
// plus petite             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<<"extrapsph [-r reduc] sphin.fits sphinw.fits"<<endl
    <<"                     sphout.fits [sphout_w.fits"
    <<" sphred.fits sphred_w.fits]"<<endl
    <<" -r reduc : reduction factor for nlat in clustering must be 2^n (def=2)"<<endl
    <<" sphin.fits : input sphere"<<endl
    <<" sphin_w.fits : input sphere filling weight"<<endl
    <<" sphout.fits  : output sphere"<<endl
    <<" sphout_w.fits  : output sphere filling weight"<<endl
    <<" sphred.fits  : output reducted sphere"<<endl
    <<" sphred_w.fits  : output reducted sphere filling weight"<<endl;
}

int main(int narg, char* arg[])
{
int red=2;
int c;
while((c = getopt(narg,arg,"hr:")) != -1) {
  switch (c) {
  case 'r' :
    sscanf(optarg,"%d",&red);
    break;
  case 'h' :
  default:
    usage(); exit(1);
  }
}

if(optind+2>=narg) {usage(); exit(1);}
char * fsphin   = arg[optind];
char * fsphinw  = arg[optind+1];
char * fsphout  = arg[optind+2];
char * fsphoutw = NULL;
char * fsphred = NULL;
char * fsphredw = NULL;
if(optind+3<narg) fsphoutw = arg[optind+3];
if(optind+4<narg) fsphred = arg[optind+4];
if(optind+5<narg) fsphredw = arg[optind+5];

cout<<"Sphere Input : "<<fsphin<<endl
    <<"Weight Sphere Input : "<<fsphinw<<endl
    <<"Sphere Output : "<<fsphout<<endl
    <<"Weight Sphere Output : "<<fsphoutw<<endl
    <<"Reducted Sphere Output : "<<fsphred<<endl
    <<"Reducted Weight Sphere Output : "<<fsphredw<<endl
    <<"Reduction : "<<red<<endl;

// Lecture des spheres
SphereHEALPix<r_8> sphin;
{FitsInFile sfits(fsphin); sfits >> sphin;}
int nlat = sphin.SizeIndex();
int nlatr = nlat/red;
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;
{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);
}

// Creation des spheres output
cout<<"Creating Spheres for output"<<endl;
SphereHEALPix<r_8> sphout(sphin,false);
SphereHEALPix<r_8> sphoutw(sphinw,false);

// Creation des spheres reduites
cout<<"Creating Reducted Spheres : nlatr = "<<nlatr<<endl;
SphereHEALPix<r_8> sphred(nlatr);
  sphred.SetPixels(0.);
SphereHEALPix<r_8> sphredw(nlatr);
  sphredw.SetPixels(0.);
cout<<"Creating Reducted Sphere :"<<endl
    <<"     Number of pixels : "<<sphred.NbPixels()<<endl
    <<"                 Nlat : "<<sphred.SizeIndex()<<endl;

// Filling reducted spheres
cout<<"...Filling Reducted Spheres"<<endl;
for(int_4 i=0;i<sphin.NbPixels();i++) {
  if(sphinw(i)<=0.) continue;
  double theta,phi;
  sphin.PixThetaPhi(i,theta,phi);
  int_4 ir = sphred.PixIndexSph(theta,phi);
  sphred(ir)  += sphin(i)*sphinw(i);
  sphredw(ir) += sphinw(i);
}
cout<<"...Computing Reducted Spheres"<<endl;
uint_4 n=0;
for(int_4 ir=0;ir<sphred.NbPixels();ir++)
  if(sphredw(ir) > 0.) {sphred(ir) /= sphredw(ir); n++;}
cout<<"      Number of filled pixels "<<n<<endl;

// Filling hole for Output Spheres
cout<<"...Filling hole for Output Spheres"<<endl;
n=0;
for(int_4 i=0;i<sphout.NbPixels();i++) {
  if(sphoutw(i)>0.) continue;
  double theta,phi;
  sphout.PixThetaPhi(i,theta,phi);
  int_4 ir = sphred.PixIndexSph(theta,phi);
  if(sphredw(ir)<=0.) continue;
  sphout(i)  = sphred(ir);
  sphoutw(i) = -sphredw(ir);
  n++;
}
cout<<"      Number of filled pixels "<<n<<endl;

// Ecriture des spheres
{
string dum = "rm -f "; dum += fsphout; system(dum.c_str());
FitsOutFile sfits(fsphout);
cout<<"Writing Output Sphere Fits file"<<endl;
sfits<<sphout;
}
if(fsphoutw) {
string dum = "rm -f "; dum += fsphoutw; system(dum.c_str());
FitsOutFile sfits(fsphoutw);
cout<<"Writing Output Sphere Weight Fits file"<<endl;
sfits<<sphoutw;
}
if(fsphred) {
string dum = "rm -f "; dum += fsphred; system(dum.c_str());
FitsOutFile sfits(fsphred);
cout<<"Writing Reducted Sphere Fits file"<<endl;
sfits<<sphred;
}
if(fsphredw) {
string dum = "rm -f "; dum += fsphredw; system(dum.c_str());
FitsOutFile sfits(fsphredw);
cout<<"Writing Reducted Sphere Weight Fits file"<<endl;
sfits<<sphredw;
}

exit(0);
}