#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>

#include <typeinfo>

#include <vector>
#include <string>

#include "piacmd.h"
#include "nobjmgr.h"
#include "pistdimgapp.h"
#include "servnobjm.h"
#include "tvector.h"
#include "pitvmaad.h"
#include "fftpserver.h"
#include "bruit.h"
#include "piscdrawwdg.h"
#include "ctimer.h"

#include "pidrawer.h"
#include "nomgadapter.h"

extern "C" {
void sopiamodule_init();
void sopiamodule_end();
}

void SophyaFFT(string& nom, string& nomout, string dopt);

class sopiamoduleExecutor : public CmdExecutor {
public:
                sopiamoduleExecutor();
  virtual       ~sopiamoduleExecutor();
  virtual int   Execute(string& keyw, vector<string>& args, string& toks);
};

//  Classe pour trace de grille en projection spherique Mollweide 
class MollweideGridDrawer : public PIDrawer {
public:
                    MollweideGridDrawer(int np, int nm, int nx, int ny);
                    MollweideGridDrawer(int np, int nm, double x_larg=0., double x_0=0., 
					double y_larg=0., double y_0=0.);
  virtual           ~MollweideGridDrawer();
  virtual void      Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax);
  virtual void      UpdateLimits();

protected:
  double xlarg, x0;
  double ylarg, y0;
  int nbparall, nbmerid;
  int nbpt;
  bool fgrevtheta;
};

/* --Methode-- */
MollweideGridDrawer::MollweideGridDrawer(int np, int nm, int nx, int ny)
{
  x0 = (double)nx/2.; 
  y0 = (double)ny/2.;
  xlarg = nx;  
  ylarg = ny;
  nbparall = np;
  nbmerid = nm;
  nbpt = ny/5;
  if ((nbpt % 20) == 0) nbpt++;  
  fgrevtheta = true;
}

/* --Methode-- */
MollweideGridDrawer::MollweideGridDrawer(int np, int nm, double x_larg, double x_0, 
					 double y_larg, double y_0)
{
  if (x_larg > 0.) {
    xlarg = x_larg; x0 = x_0;
  }
  else {
    xlarg = M_PI*2.;  x0 = M_PI;
  }
  if (y_larg > 0.) {
    ylarg = y_larg; y0 = y_0;
  }
  else {
    ylarg = M_PI;  y0 = 0.;
  }
  nbparall = np;
  nbmerid = nm;
  nbpt = 101;
  fgrevtheta = false;
}

/* --Methode-- */
MollweideGridDrawer::~MollweideGridDrawer()
{
}

/* --Methode-- */
void MollweideGridDrawer::Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax)
{
  PIGrCoord * xxl = new PIGrCoord[nbpt];
  PIGrCoord * xxr = new PIGrCoord[nbpt];
  PIGrCoord * yy = new PIGrCoord[nbpt];

  char buff[64];

  // Trace des paralleles
  g->DrawLine(xmin, y0, xmax, y0);
  g->DrawString(x0+xlarg/100, y0+ylarg/100, "0");

  //  PIGrCoord asc, PIGrCoord desc;
  int i,k;

  int nkp = (nbparall-1)/2;
  if (nkp > 0) {
    double dy = 0.5*ylarg/(double)(nkp+1);
    double dtd = 90./(double)(nkp+1);
    for(k=1; k<=nkp; k++) {
      double fx = sin(acos(2.*(double)k*dy/ylarg))*0.5;
      double yc = k*dy+y0;
      g->DrawLine(x0-fx*xlarg, yc, x0+fx*xlarg, yc); 
      double ctd = (fgrevtheta) ? -dtd*k : dtd*k;
      sprintf(buff, "%5.1f", ctd);
      g->DrawString(x0, yc, buff);
      g->DrawString(x0-0.3*xlarg, yc, buff);
      g->DrawString(x0+0.25*xlarg, yc, buff);
      yc = -k*dy+y0;
      g->DrawLine(x0-fx*xlarg, yc, x0+fx*xlarg, yc); 
      sprintf(buff, "%5.1f", -ctd);
      g->DrawString(x0, yc, buff);
      g->DrawString(x0-0.3*xlarg, yc, buff);
      g->DrawString(x0+0.25*xlarg, yc, buff);
    }
  }

  // Trace des meridiennes
  g->DrawLine(x0, ymin, x0, ymax);

  int nkm = (nbmerid-1)/2;
  if (nkm < 1)  nkm = 1;
  double dphi = M_PI/nkm;
  double dpd = 180./(double)(nkm);
  double dy = ylarg/(nbpt-1);
  for(k=0; k<nkm; k++) {
    double phi0 = k*dphi;
    double dphimil = fabs(phi0-M_PI);
    if (dphimil < 1.e-6)  continue;
    for(i=0; i<nbpt; i++) {
      double yc = (double)i*dy+y0-0.5*ylarg;
      double fx = (2.*fabs(yc-y0)/ylarg <= 1.) ? 
	sin(acos(2.*fabs(yc-y0)/ylarg)) * dphimil/M_PI * 0.5 : 0.;
      yy[i] = yc;
      xxl[i] = x0-fx*xlarg;
      xxr[i] = x0+fx*xlarg;
    }

    g->DrawPolygon(xxl, yy, nbpt, false);
    g->DrawPolygon(xxr, yy, nbpt, false);

    if (k == 0)  continue;
    for(i=-1; i<=1; i++) {
      double yc = 0.20*i*ylarg+y0;
      double fx = sin(acos(2.*fabs(yc-y0)/ylarg)) * dphimil/M_PI * 0.5;
      double xc = x0-fx*xlarg;
      sprintf(buff, "%4.1f", dpd*k);
      g->DrawString(xc, yc, buff);
      xc = x0+fx*xlarg;
      sprintf(buff, "%4.1f", 360.-dpd*k);
      g->DrawString(xc, yc, buff);
    }
  }

  delete [] xxl;
  delete [] xxr;
  delete [] yy;
}

/* --Methode-- */
void MollweideGridDrawer::UpdateLimits()
{
  SetLimits(x0-xlarg/2., x0+xlarg/2., y0-ylarg/2., y0+ylarg/2.);
}

/* --Methode-- */
sopiamoduleExecutor::sopiamoduleExecutor()
{

PIACmd * mpiac;
NamedObjMgr omg;
mpiac = omg.GetImgApp()->CmdInterpreter();

string hgrp = "SophyaCmd";
string kw = "powerspec";
string usage = "FFT on a vector -> Plots power spectrum ";
usage += "\n Usage: fftp vecName vecFFT [graphic_att] ";
mpiac->RegisterCommand(kw, usage, this, hgrp);
kw = "mollgridsph";
usage = "Creates a spherical coordinate grid in Molleweide projection ";
usage += "\n Usage: mollgridsph NameSphericalMap [Nb_Parallel Nb_Meridien graphic_att] ";
mpiac->RegisterCommand(kw, usage, this, hgrp);
kw = "mollgrid";
usage = "Creates a spherical coordinate grid in Molleweide projection ";
usage += "\n Usage: mollgrid [Nb_Parallel Nb_Meridien graphic_att] ";
mpiac->RegisterCommand(kw, usage, this, hgrp);

}

/* --Methode-- */
sopiamoduleExecutor::~sopiamoduleExecutor()
{
}

/* --Methode-- */
int sopiamoduleExecutor::Execute(string& kw, vector<string>& tokens, string& toks)
{

if (kw == "powerspec") {
  if (tokens.size() < 2) {
    cout << "Usage: powerspec nameVec vecFFT [graphic_att]" << endl;
    return(0);
    }
  if (tokens.size() < 3)  tokens.push_back((string)"n");
  SophyaFFT(tokens[0], tokens[1], tokens[2]);
  }
 else if ( (kw == "mollgridsph") || (kw == "mollgrid") ) {
   NamedObjMgr omg;
   MollweideGridDrawer * mollgrid = NULL;
   string dopt="";
   if (kw == "mollgridsph") {
     if (tokens.size() < 1) {
       cout << "Usage: mollgridsph NameSphericalMap [Nb_Parallel Nb_Meridien graphic_att]"
	    << endl;
       return(0);
     }
     int np = 5;
     int nm = 5;
     dopt = "same,thinline";
     if (tokens.size() > 1)  np = atoi(tokens[1].c_str());
     if (tokens.size() > 2)  nm = atoi(tokens[2].c_str());
     if (tokens.size() > 3)  dopt = tokens[3];
     NObjMgrAdapter* obja = omg.GetObjAdapter(tokens[0]);
     if (obja == NULL) {
       cout << "mollgridsph Error , No object with name " << tokens[0] << endl; 
       return(0);
     }
     P2DArrayAdapter* arr = obja->Get2DArray(dopt);
     if (arr == NULL) {
       cout << "mollgridsph Error , object has non 2DArrayAdapter - Not a spherical map "
	    << endl;
       return(0);       
     }
     int nx = arr->XSize();
     int ny = arr->YSize();
     delete arr;
     mollgrid = new MollweideGridDrawer(np, nm, nx, ny);
     cout << " mollgridsph: Creating MollweideGridDrawer() for object" << tokens[0] << endl;
   }
   else if (kw == "mollgrid") {
     int np = 5;
     int nm = 5;
     if (tokens.size() > 0)  np = atoi(tokens[0].c_str());
     if (tokens.size() > 1)  nm = atoi(tokens[1].c_str());
     if (tokens.size() > 2)  dopt = tokens[2];
     mollgrid = new MollweideGridDrawer(np, nm);
     cout << " mollgrid: Creating MollweideGridDrawer(" << np << "," << nm << ")" << endl;
   }

   if (mollgrid == NULL) {
     cout << "mollgridsph Error/BUG !!! mollgrid = NULL "  << endl; 
     return(0);
   }

   int wrsid = 0;
   //bool fgsr = true;
   
   string name = "mollgrid";
   wrsid = omg.GetImgApp()->DispScDrawer(mollgrid, name, dopt); 
 }

return(0);

}

static sopiamoduleExecutor * piaerex = NULL;
/* Nouvelle-Fonction */
void sopiamodule_init()
{
// Fonction d'initialisation du module
// Appele par le gestionnaire de modules de piapp (PIACmd::LoadModule())
if (piaerex) delete piaerex;
piaerex = new sopiamoduleExecutor;
}

/* Nouvelle-Fonction */
void sopiamodule_end()
{
// Desactivation du module
if (piaerex) delete piaerex;
piaerex = NULL;
}


/* --Methode-- */
void SophyaFFT(string& nom, string& nomout, string dopt)
{
Timer tm("powerspec");
NamedObjMgr omg;
AnyDataObj* obj=omg.GetObj(nom);
if (obj == NULL) {
  cout << "SophyaFFT() Error , Pas d'objet de nom " << nom << endl;
  return;
}

Vector* vin = dynamic_cast<Vector *>(obj);
if (vin == NULL) {
  cout << "SophyaFFT() Error , Objet n'est pas un vecteur " <<  endl;
  return;
  }
TVector< complex<double> > * vout = new TVector< complex<double> >  ;
FFTPackServer ffts;
cout << "SophyaFFT() - Computing FFT of vector size " << vin->NElts() << endl;
ffts.FFTForward(*vin, *vout);
// cout << " SophyaFFT - FFT done " << endl;
 // tm.Split(" FFT done ");

omg.AddObj(vout, nomout);
omg.DisplayObj(nomout, dopt);
return;

}