#include "machdefs.h"    // Definitions specifiques SOPHYA

#include <math.h>
#include <iostream.h>

#include "nbmath.h"
#include "timing.h"

#include "array.h"
#include "skymap.h"
#include "samba.h"
#include "sambainit.h"
#include "fitsspherehealpix.h"
#include "fitstarray.h"

/*!
  \ingroup PrgMap
  \file prjsmap.cc
  \brief \b prjsmap: Molleweide and sinus projection of sky maps

  \verbatim 

  csh> prjsmap -h
  SOPHYA Version  1.1 Revision 0 (V_Fev2001) -- Mar  9 2001 15:45:31 cxx

  prjsmap : Molleweide and Sinus projection for sky maps
  Usage: prjsmap [-float/-r4] [-sinus] [-sx sizeX] [-sy sizeY] [-scale fact]
         [-off val][-ump val] [-invy] [-fitsin] [-fitsout] InFileName OutFileName
   -float (-r4): single precision sky map and projection (default = double)
   -sinus: Selects Sinus projection (default Molleweide projection)
   -sx sizeX: Projected map size in X
   -sy sizeY: Projected map size in Y
       default: sizeX = 2*sizeY ; sizeX*sizeY ~ map.NbPixels()/2
   -scale fact: scale factor applied to skymap pixels (def=1.)
   -off val: offset value for projected map (def=0.)
             prj(ix, jy) = offset + scale*skymap(theta, phi)
   -ump val: Unmapped pixel value in projected map (def=0.)
   -invy: reverses the projection direction along Y (Theta)
   -fitsout: Select the FITS format for the output map (default PPF format)
   -fitsin : Select the FITS format for the input vector (default PPF format)
   InFileName : Input file name (HEALPix map)
   OutFileName : Output file name (the C_l vector)

  \endverbatim
*/


/* Nouvelle-Fonction */
void Usage(bool fgerr)
{
  cout << endl;
  if (fgerr) {
    cout << " prjsmap : Argument Error ! prjsmap -h for usage " << endl; 
    exit(1);
  }
  else {
    cout << " prjsmap : Molleweide and Sinus projection for sky maps \n" 
	 << " Usage: prjsmap [-float/-r4] [-sinus] [-sx sizeX] [-sy sizeY] [-scale fact] \n" 
	 << "        [-off val][-ump val] [-invy] [-fitsin] [-fitsout] InFileName OutFileName\n" 
	 << "   -float (-r4): single precision sky map and projection (default = double) \n"
	 << "   -sinus: Selects Sinus projection (default Molleweide projection) \n"
	 << "   -sx sizeX: Projected map size in X \n"
	 << "   -sy sizeY: Projected map size in Y \n"
	 << "       default: sizeX = 2*sizeY ; sizeX*sizeY ~ map.NbPixels()/2 \n"
	 << "   -scale fact: scale factor applied to skymap pixels (def=1.) \n"
	 << "   -off val: offset value for projected map (def=0.) \n"
	 << "             prj(ix, jy) = offset + scale*skymap(theta, phi) \n"
	 << "   -ump val: Unmapped pixel value in projected map (def=0.) \n" 
	 << "   -invy: reverses the projection direction along Y (Theta) \n"
	 << "   -fitsout: Select the FITS format for the output map (default PPF format) \n"
	 << "   -fitsin : Select the FITS format for the input vector (default PPF format) \n"
	 << "   InFileName : Input file name (HEALPix map) \n" 
	 << "   OutFileName : Output file name (the C_l vector) \n" << endl; 
    exit(0);
  }  
}

/* Nouvelle-Fonction */
template <class T> 
class _PrjMap {
public :
/* Methode */
static void PM_MeanSigma(PixelMap<T> const & map, double& gmoy, double& gsig, T& min, T& max)
{
  min = max = map(0);
  gmoy=0.;
  gsig = 0.; 
  double valok;
  for(int k=0; k<map.NbPixels(); k++) {
    valok = map(k);
    gmoy += valok;  gsig += valok*valok; 
    if (map(k)<min)  min = map(k);
    else if (map(k)>max)  max = map(k);      
  }
  gmoy /= (double)map.NbPixels();
  gsig = gsig/(double)map.NbPixels() - gmoy*gmoy;
  if (gsig >= 0.) gsig = sqrt(gsig);
  return;
}

/* Methode */
static TArray<T> Project_Molleweide(SphericalMap<T> const & map, int sx=0, int sy=0, T ump=0, 
				    T offset=0, T scale=1, bool invy=false)
{
  if ( (sx <= 0) && (sy <= 0) ) {
    sy = sqrt((double)(map.NbPixels()/2));
    sx = sy*2;    
  }
  else {
    if (sx <= 0) sx = 2*sy;
    else sy = sx/2;
  }
  TArray<T> prj(sx, sy);
  prj = ump;   // On met tout a ump
  
  cout << " Molleweide projection (sizeX= " << sx << " sizeY= " << sy << " )" 
       << " Scale=" << scale << " Offset=" << offset << endl;
  r_8 xa, yd, teta,phi, facteur;
  int_4 ix, jy, k;
  for(jy=0; jy<sy; jy++) {
    yd = (r_8)(jy+0.5)/(r_8)sy-0.5;
    if (invy) teta = (0.5-yd)*M_PI;
    else teta = (yd+0.5)*M_PI;
    facteur=2.*M_PI/sin(acos((double)yd*2));
    for(ix=0; ix<sx; ix++)  {
      xa = (r_8)(ix+0.5)/(r_8)sx-0.5;
      phi = xa*facteur+M_PI;
      if ( (phi <= 2*M_PI) && (phi >= 0.) ) {
        k = map.PixIndexSph(teta, phi);
        prj(ix,jy,0) = offset + scale*map(k);
      }
    }
  }

  return prj;   
}

/* Methode */
static TArray<T> Project_Sinus(SphericalMap<T> const & map, int sx=0, int sy=0, T ump=0, 
			       T offset=0, T scale=1, bool invy=false)
{
  if ( (sx <= 0) && (sy <= 0) ) {
    sy = sqrt((double)(map.NbPixels()/2));
    sx = sy*2;    
  }
  else {
    if (sx <= 0) sx = 2*sy;
    else sy = sx/2;
  }
  TArray<T> prj(sx, sy);
  prj = ump;   // On met tout a ump

  cout << " Sinus projection (sizeX= " << sx << " sizeY= " << sy << " )"
       << " Scale=" << scale << " Offset=" << offset << endl;
  r_8 xa, yd, teta,phi, facteur;
  int_4 ix, jy, k;
  for(jy=0; jy<sy; jy++) {
    yd = (r_8)(jy+0.5)/(r_8)sy-0.5;
    if (invy) teta = (0.5-yd)*M_PI;
    else teta = (yd+0.5)*M_PI;
    facteur=1./sin(teta);
    for(ix=0; ix<sx; ix++)  {
      xa = (r_8)(ix+0.5)/(r_8)sx-0.5;
      phi = 2*M_PI*xa*facteur+M_PI;
      if ( (phi <= 2*M_PI) && (phi >= 0.) ) {
        k = map.PixIndexSph(teta, phi);
        prj(ix,jy,0) = offset + scale*map(k);
      }
    }
  }
  return prj;   
}

/* Methode */
static void Project(string & infile, string & outfile, int sx, int sy, T ump, T offset, 
                    T scale, bool fgsin, bool invy, bool fgfitsin, bool fgfitsout) 
{

  double deg2rad =  M_PI/180.;
  double minute2rad =  M_PI/(180.*60.);

  SphericalMap<T> * sphp = NULL;
  SphereHEALPix<T> sphh;

  PPersist * pp = NULL;

  if (fgfitsin) {
    cout << "--- Reading Input FITS file " << infile << endl; 
    FitsInFile fii(infile);
    fii >> sphh;
    sphp = &sphh;
  }
  else {
    cout << "--- Reading Input PPF file " << infile << endl; 
    PInPersist ppi(infile);
    pp = ppi.ReadObject();
    if (pp) 
      sphp = dynamic_cast< SphericalMap<T> * > (pp->DataObj());
    else sphp = NULL;
    if (sphp == NULL) {
      cout << " Object in PPF file is not a SphericalMap<T> or wrong data type ! " << endl;
      throw TypeMismatchExc("_PrjMap::Project(...) Error reading input PPF file !");
    }
  }

  SphericalMap<T> & sph = *sphp;
  T min, max;
  double mean, sigma;
  PM_MeanSigma(sph, mean, sigma, min, max);
  
  cout << " Input map type= " << sph.TypeOfMap() << " PixParm (=NSide for HEALPix)= " 
       << sph.SizeIndex() << endl;
  cout << " Input map : NbPixels= " << sph.NbPixels() << " Resolution= "
       << sqrt(sph.PixSolAngle(0))/minute2rad << " Arcminutes " << endl;
  cout << " map.Min= " << min << " map.Max= " << max 
       << " map.Mean= " << mean << " map.Sigma= " << sigma << endl;
  
  cout << "--- Sky map projection " << endl;
  TArray<T> prj;

  if (fgsin) prj = Project_Sinus(sph, sx, sy, ump, offset, scale, invy); 
  else prj = Project_Molleweide(sph, sx, sy, ump, offset, scale, invy); 

  cout << "--- Statistics on the projected map :" ; 
  prj.Show();
  prj.MinMax(min, max);
  MeanSigma(prj, mean, sigma);  
  cout << " prj.Min= " << min << " prj.Max= " << max 
       << " prj.Mean= " << mean << " prj.Sigma= " << sigma << endl;
  
  if (fgfitsout) {
    cout << "--- Writing projection to Output FITS file " << outfile << endl; 
    FitsOutFile fio(outfile);
    fio << prj;
  }
  else {
    cout << "--- Writing projection to Output PPF file " << outfile << endl; 
    POutPersist ppo(outfile);
    ppo << prj;
  }

  if (pp) delete pp;
}

};

/* Main-Program */
int main(int narg, char *arg[])
{
  if ((narg > 1) && (strcmp(arg[1],"-h") == 0) ) Usage(false);
  
  double ump = 0.;
  double off = 0.;
  double scale = 1.;
  int sx = 0;
  int sy = 0;
  string infile;
  string outfile;
  bool fgfitsin = false;
  bool fgfitsout = false;
  bool fgr4 = false;
  bool fgsinus = false;
  bool fginvy = false;

  cout << " prjsmap : Decoding command line options ... " << endl;

  int ko=1;
  for (int k=1; k<narg; k++)   {
    if (strcmp(arg[k], "-sx") == 0)  {
      if (k == narg-1) Usage(true);  // -sx est suivi d'un argument 
      sx = atoi(arg[k+1]);  k++;       // k++ pour sauter au suivant
    }
    if (strcmp(arg[k], "-sy") == 0)  {
      if (k == narg-1) Usage(true);  // -sy est suivi d'un argument 
      sy = atoi(arg[k+1]);  k++;       // k++ pour sauter au suivant
    }
    else if (strcmp(arg[k], "-scale") == 0)  {
      if (k == narg-1) Usage(true);  
      scale = atof(arg[k+1]);  k++;   
    }
    else if (strcmp(arg[k], "-off") == 0)  {
      if (k == narg-1) Usage(true);  
      off = atof(arg[k+1]);  k++;   
    }
    else if (strcmp(arg[k], "-ump") == 0)  {
      if (k == narg-1) Usage(true);  
      ump = atof(arg[k+1]);  k++;   
    }
    else if (strcmp(arg[k], "-sinus") == 0) {
      fgsinus = true;  
    }
    else if (strcmp(arg[k], "-invy") == 0) {
      fginvy = true;  
    }
    else if (strcmp(arg[k], "-fitsin") == 0) {
      fgfitsin = true;  
    }
    else if (strcmp(arg[k], "-fitsout") == 0) {
      fgfitsout = true;  
    }
    else if ((strcmp(arg[k], "-float") == 0) || (strcmp(arg[k], "-r4") == 0) ) {
      fgr4 = true;  
    }
   
    else { ko = k; break; }  // Debut des noms
  }

  if ((narg-ko) < 2)  Usage(true); 
  infile = arg[ko];
  outfile = arg[ko+1];
  
  // Bloc try de gestion des exception 
  try {
    InitTim();
    SophyaInit();
    if (fgr4) {
      cout << " Projecting SphericalMap<r_4>  -->  Array<r_4> (float)" << endl;
      _PrjMap<r_4>::Project(infile, outfile, sx, sy, ump, off, scale, 
			    fgsinus, fginvy, fgfitsin, fgfitsout);
    }
    else {
      cout << " Projecting SphericalMap<r_8>  -->  Array<r_8> (double)" << endl;
      _PrjMap<r_8>::Project(infile, outfile, sx, sy, ump, off, scale, 
			    fgsinus, fginvy, fgfitsin, fgfitsout);
    }
  }
  catch (PThrowable & exc) {   // Exceptions de SOPHYA
    cerr << " prjsmap: Catched Exception " << (string)typeid(exc).name()
         << " - Msg= " << exc.Msg() << endl;
  }
  catch (...) {    // Autres Exceptions
    cerr << " prjsmap: some other exception was caught ! " << endl;
  }
  
  PrtTim("End of prjsmap ");
  return(0);
}