#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 "pidrawer.h"
#include "nomgadapter.h"

#include "skymap.h"
#include "sphericaltransformserver.h"

static bool IsNSideGood(int_4 nside) {
 if(nside<=1) return false;
 while(nside>1) {if(nside%2!=0) return false; else nside/=2;}
 return true;
}

extern "C" {
  void skymapmodule_init();
  void skymapmodule_end();
}

class skymapmoduleExecutor : public CmdExecutor {
public:
                skymapmoduleExecutor();
  virtual       ~skymapmoduleExecutor();
  virtual int   Execute(string& keyw, vector<string>& args, string& toks);
  void Map2Double(vector<string>& tokens);
  void MapMult(vector<string>& tokens);
  void MapCover(vector<string>& tokens);
  void Map2Cl(vector<string>& tokens);
  void Cl2Map(vector<string>& tokens);
  void Map2Alm(vector<string>& tokens);
  void Alm2Map(vector<string>& tokens);
  void CrMapMask(vector<string>& tokens);
  void CrMaskFrMap(vector<string>& tokens);
  void MaskMap(vector<string>& tokens);
  void MapProj(vector<string>& tokens);
  void Map2Local(vector<string>& tokens);
  void MapOper(vector<string>& tokens);
  void MapStat(vector<string>& tokens);
};

/* --Methode-- */
skymapmoduleExecutor::skymapmoduleExecutor()
{
PIACmd * mpiac;
NamedObjMgr omg;
mpiac = omg.GetImgApp()->CmdInterpreter();

string hgrp = "SkyMapCmd";

string kw = "map2double";
string usage = "Convert a float map to a double map";
usage += "\n Usage: map2double map";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "mapmult";
usage = "Multiply a map by a constant";
usage += "\n Usage: mapmult map cste";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "mapcover";
usage = "Print the covering of a map";
usage += "\n Usage: mapcover map v1,v2 [varname]";
usage += "\n    v1,v2: good pixels have v1<=val<=v2 (def=0.99,1.01)";
usage += "\n          (use default (0.99,1.01): \"!\")";
usage += "\n    $varname: percent of sky covering";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "map2cl";
usage = "SkyMap to Cl";
usage += "\n Usage: map2cl map(double) clvec [lmax] [niter]";
usage += "\n   lmax: if <=0 or \"!\" use default (3*nside-1)";
usage += "\n   niter: number of iterations (<=0 or def=3)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "cl2map";
usage = "Generate SkyMap from Cl";
usage += "\n Usage: cl2map clvec map nside [fwhm]";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "map2alm";
usage = "SkyMap to Alm";
usage += "\n Usage: map2alm map(double) almmat [lmax] [niter]";
usage += "\n   lmax: if <=0 or \"!\" use default (3*nside-1)";
usage += "\n   niter: number of iterations (<=0 or def=3)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "alm2map";
usage = "Generate SkyMap from Alm";
usage += "\n Usage: alm2map almmat map nside";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "crmapmask";
usage = "Create a map mask (nside) between latitudes lat1,lat2 longitudes lon1,lon2";
usage += "\n Usage: crmapmask mapmsk nside lat1,lat2 lon1,lon2 [valmsk,valnomsk]";
usage += "\n    mask pixels between [lat1,lat2] ([-90,90] in degrees)";
usage += "\n                and    [lon1,lon2] ([0,360[ in degrees)";
usage += "\n                (for no mask \"!\")";
usage += "\n    valmsk=value to be written into masked pixels (def=0)";
usage += "\n    valnomsk=value to be written into unmasked pixels (def=1)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "crmaskfrmap";
usage = "Create a map mask (nside) relative to an other map pixels values";
usage += "\n Usage: crmaskfrmap mapmsk nside map(double) v1,v2 [valmsk,valnomsk]";
usage += "\n    mask if v1<mapmsk(i)<v2";
usage += "\n    valmsk=value to be written into masked pixels (def=0)";
usage += "\n    valnomsk=value to be written into unmasked pixels (def=1)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "maskmap";
usage = "Mask a map with a mask map";
usage += "\n Usage: maskmap map msk";
usage += "\n    operation is map(i) *= msk(theta,phi)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "maproj";
usage = "Project a map into an other one";
usage += "\n Usage: maproj map mapproj [nside]";
usage += "\n   Create mapproj(nside) and project map into mapproj";
usage += "\n   (use \"maproj map mapproj\" to make a copy)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "map2local";
usage = "Project a map into a local map";
usage += "\n Usage: map2local map localmap nx,ny angleX,angleY phi0,theta0 [x0,y0] [angle]";
usage += "\n   nx,ny: number of pixels in x(col),y(row) direction";
usage += "\n          X-axis==phi, Y-axis==theta";
usage += "\n   angleX,angleY: total angle aperture in x,y direction (degrees)";
usage += "\n   phi0,theta0: define origin in phi,theta (degrees)";
usage += "\n   x0,y0: map phi0,theta0 to pixel x0,y0 (\"!\" or def=middle of the localmap)";
usage += "\n   angle: angle (degrees) of rotation between x-axis of  map-coordinate system";
usage += "\n           and the tangent to parallel on the sphere (def: 0.)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "mapop";
usage = "operation on a map";
usage += "\n Usage: mapop map op1 map1 [op2 map2] [op2 map3] [...]";
usage += "\n   Do \"map op1= map1\", \"map op2=map2\", ...";
usage += "\n      where op = +,-,*,/";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "mapstat";
usage = "Statistics from a map";
usage += "\n Usage: mapstat map [msk] [meanvarname] [sigmavarname]";
usage += "\n   msk = mask sphere used as weight";
usage += "\n     if msk(i)<=0 do not use that pixel";
usage += "\n     if msk(i)>0  use that pixel with weight msk(i)";
usage += "\n   msk = \"!\" means no mask sphere";
usage += "\n   meanvarname  = variable name to store mean";
usage += "\n   sigmavarname = variable name to store sigma";
usage += "\n   (\"!\" means do not store)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

}

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

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

if(kw == "map2double") {
  if(tokens.size()<1) {
    cout<<"Usage: map2double map"<<endl;
    return(0);
  }
  Map2Double(tokens);
} else if(kw == "mapmult") {
  if(tokens.size()<2) {
    cout<<"Usage: mapmult map cste"<<endl;
    return(0);
  }
  MapMult(tokens);
} else if(kw == "mapcover") {
  if(tokens.size()<1) {
    cout<<"Usage: mapcover map v1,v2 [varname]"<<endl;
    return(0);
  }
  MapCover(tokens);
} else if(kw == "map2cl") {
  if(tokens.size()<2) {
    cout<<"Usage: map2cl map clvec [lmax] [niter]"<<endl;
    return(0);
  }
  Map2Cl(tokens);
} else if(kw == "cl2map") {
  if(tokens.size()<2) {
    cout<<"Usage: cl2map clvec map nside [fwhm]"<<endl;
    return(0);
  }
  Cl2Map(tokens);
} else if(kw == "map2alm") {
  if(tokens.size()<2) {
    cout<<"Usage: map2alm map(double) almmat [lmax] [niter]"<<endl;
    return(0);
  }
  Map2Alm(tokens);
} else if(kw == "alm2map") {
  if(tokens.size()<3) {
    cout<<"Usage: alm2map almmat map nside"<<endl;
    return(0);
  }
  Alm2Map(tokens);
} else if(kw == "crmapmask") {
  if(tokens.size()<2) {
    cout<<"Usage: crmapmask mapmsk nside lat1,lat2 lon1,lon2 [valmsk,valnomsk]"<<endl;
    return(0);
  }
  CrMapMask(tokens);
} else if(kw == "crmaskfrmap") {
  if(tokens.size()<3) {
    cout<<"Usage: crmaskfrmap mapmsk nside map(double) v1,v2 [valmsk,valnomsk]"<<endl;
    return(0);
  }
  CrMaskFrMap(tokens);
} else if(kw == "maskmap") {
  if(tokens.size()<2) {
    cout<<"Usage: maskmap map msk"<<endl;
    return(0);
  }
  MaskMap(tokens);
} else if(kw == "maproj") {
  if(tokens.size()<2) {
    cout<<"Usage: maproj map mapproj [nside]"<<endl;
    return(0);
  }
  MapProj(tokens);
} else if(kw == "map2local") {
  if(tokens.size()<5) {
    cout<<"Usage: map2local map localmap nx,ny angleX,angleY phi0,theta0 [x0,y0] [angle]"<<endl;
    return(0);
  }
  Map2Local(tokens);
} else if(kw == "mapop") {
  if(tokens.size()<3) {
    cout<<"Usage: mapop map op map1 [op map2] [op map3] [...]"<<endl;
    return(0);
  }
  MapOper(tokens);
} else if(kw == "mapstat") {
  if(tokens.size()<1) {
    cout<<"Usage: Usage: mapstat map [msk] [meanvarname] [sigmavarname]"<<endl;
    return(0);
  }
  MapStat(tokens);
}

return(0);
}

/* --Methode-- */
void skymapmoduleExecutor::Map2Double(vector<string>& tokens)
{
 NamedObjMgr omg;
 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_4>* map = dynamic_cast<SphereHEALPix<r_4>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_4>"<<endl;
   return;
 }
 int_4 nside = map->SizeIndex();
 SphereHEALPix<r_8>* map8 = new SphereHEALPix<r_8>(nside);
 cout<<"Map2Double: converting to double "<<tokens[0]<<endl;
 for(int_4 i=0;i<map8->NbPixels();i++) (*map8)(i) = (r_8) (*map)(i);
 omg.AddObj(map8,tokens[0]);
}

/* --Methode-- */
void skymapmoduleExecutor::MapMult(vector<string>& tokens)
{
 NamedObjMgr omg;
 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }
 double v = atof(tokens[1].c_str());
 cout<<"MapMult: Multiply "<<tokens[0]<<" by "<<v<<endl;
 for(int_4 i=0;i<map->NbPixels();i++) (*map)(i) *= v;
}

/* --Methode-- */
void skymapmoduleExecutor::MapCover(vector<string>& tokens)
{
 NamedObjMgr omg;
 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 double v1=0.99, v2=1.01;
 if(tokens.size()>1) if(tokens[1][0]!='!')
   sscanf(tokens[1].c_str(),"%lf,%lf",&v1,&v2);
 cout<<"MapCover: "<<tokens[0]<<" good px=["<<v1<<","<<v2<<"]"<<endl;
 int_4 ngood=0;
 for(int_4 i=0;i<map->NbPixels();i++)
   if(v1<=(*map)(i) && (*map)(i)<=v2) ngood++;
 double per = (double)ngood/(double)map->NbPixels();
 cout<<"NGood="<<ngood<<" NTot="<<map->NbPixels()
     <<" -> "<<100.*per<<" %";
 if(tokens.size()>2) {
   if(omg.HasVar(tokens[2])) omg.DeleteVar(tokens[2]);
   char str[128]; sprintf(str,"%g",per);
   omg.SetVar(tokens[2],(string)str);
   cout<<" stored into $"<<tokens[2];
 }
 cout<<endl;
}

/* --Methode-- */
void skymapmoduleExecutor::Map2Cl(vector<string>& tokens)
{
 NamedObjMgr omg;

 int_4 lmax=0, niter=3;
 if(tokens.size()>2) if(tokens[2][0]!='!')
   lmax = atoi(tokens[2].c_str());
 if(tokens.size()>3) niter = atoi(tokens[3].c_str());
 if(niter<=0) niter = 3;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 SphericalTransformServer<r_8> ylmserver;

 int_4 nside = map->SizeIndex();
 if(lmax<=0) lmax = 3*nside-1;

 cout<<"Map2Cl: "<<tokens[0]<<" -> "<<tokens[1]
     <<" lmax="<<lmax<<" niter="<<niter<<endl;

 SphereHEALPix<r_8> tmpmap(*map,false);
 TVector<r_8> cltmp = ylmserver.DecomposeToCl(tmpmap,lmax,0.,niter);

 TVector<r_8>* cl = new TVector<r_8>(cltmp,false);
 omg.AddObj(cl,tokens[1]);
}

/* --Methode-- */
void skymapmoduleExecutor::Cl2Map(vector<string>& tokens)
{
 NamedObjMgr omg;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 TVector<r_8>* cl = dynamic_cast<TVector<r_8>*>(obj);
 if(cl==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a TVector<r_8>"<<endl;
   return;
 }

 int nside = 128;
 if(tokens.size()>2) nside = atoi(tokens[2].c_str());
 double fwhm = 0.;
 if(tokens.size()>3) fwhm = atof(tokens[3].c_str());

 SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>;

 SphericalTransformServer<r_8> ylmserver;
 cout<<"Cl2Map: Generate map "<<tokens[1]<<" with nside="<<nside
     <<" from cl "<<tokens[0]<<" with fwhm="<<fwhm<<endl;
 ylmserver.GenerateFromCl(*map,nside,*cl,fwhm);

 omg.AddObj(map,tokens[1]);
}

/* --Methode-- */
void skymapmoduleExecutor::Map2Alm(vector<string>& tokens)
{
 NamedObjMgr omg;

 int_4 lmax=0, niter=3;
 if(tokens.size()>2) if(tokens[2][0]!='!')
   lmax = atoi(tokens[2].c_str());
 if(tokens.size()>3) niter = atoi(tokens[3].c_str());
 if(niter<=0) niter = 3;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 SphericalTransformServer<r_8> ylmserver;

 int_4 nside = map->SizeIndex();
 if(lmax<=0) lmax = 3*nside-1;

 cout<<"Map2Alm: "<<tokens[0]<<" -> "<<tokens[1]
     <<" lmax="<<lmax<<" niter="<<niter<<endl;

 SphereHEALPix<r_8> tmpmap(*map,false);
 Alm<r_8> almtmp;
 ylmserver.DecomposeToAlm(tmpmap,almtmp,lmax,0.,niter);

 int n = almtmp.rowNumber();
 TMatrix< complex<r_8> >* alm = new TMatrix< complex<r_8> >(n,n);
 *alm = (complex<r_8>)0.;
 for(int i=0;i<n;i++) for(int j=0;j<=i;j++) (*alm)(i,j)=almtmp(i,j);
 omg.AddObj(alm,tokens[1]);
}

/* --Methode-- */
void skymapmoduleExecutor::Alm2Map(vector<string>& tokens)
{
 NamedObjMgr omg;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 TMatrix< complex<r_8> >* almmat = dynamic_cast<TMatrix< complex<r_8> >*>(obj);
 if(almmat==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a TMatrix< complex<r_8> >"<<endl;
   return;
 }

 int lmax = almmat->NRows()-1;
 Alm<r_8> alm(lmax);
 alm = (complex<r_8>)0.;
 for(int i=0;i<lmax+1;i++) for(int j=0;j<=i;j++) alm(i,j)=(*almmat)(i,j);

 int nside = 128;
 if(tokens.size()>2) nside = atoi(tokens[2].c_str());
 SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>;

 SphericalTransformServer<r_8> ylmserver;
 cout<<"Alm2Map: Generate map "<<tokens[1]<<" with nside="<<nside
     <<" from alm "<<tokens[0]<<endl;
 ylmserver.GenerateFromAlm(*map,nside,alm);

 omg.AddObj(map,tokens[1]);
}

/* --Methode-- */
void skymapmoduleExecutor::CrMapMask(vector<string>& tokens)
{
 NamedObjMgr omg;

 int_4 nside = atoi(tokens[1].c_str());
 if(!IsNSideGood(nside)) {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}

 SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>(nside);

 bool lat=false, lon=false;
 double lat1=-99., lat2=99., lon1=-999., lon2=999.;
 if(tokens.size()>2) if(tokens[2][0]!='!')
   {lat=true; sscanf(tokens[2].c_str(),"%lf,%lf",&lat1,&lat2);}
 if(tokens.size()>3) if(tokens[3][0]!='!')
   {lon=true; sscanf(tokens[3].c_str(),"%lf,%lf",&lon1,&lon2);}

 double valmsk=0., unvalmsk=1.;
 if(tokens.size()>4) sscanf(tokens[4].c_str(),"%lf,%lf",&valmsk,&unvalmsk);

 cout<<"CrMapMask "<<tokens[0]<<" nside="<<nside;
 if(lat) cout<<" for latitude in ["<<lat1<<","<<lat2<<"]";
 if(lon) cout<<" for longitude in ["<<lon1<<","<<lon2<<"]";
 cout<<" valmsk="<<valmsk<<" unvalmsk="<<unvalmsk<<endl;

 //Attention: conversion en co-latitude ==> inversion: [lat2,lat1]
 lat1 = (90.-lat1)*M_PI/180.; lat2 = (90.-lat2)*M_PI/180.;
 lon1 *= M_PI/180.; lon2 *= M_PI/180.;
 for(int_4 i=0;i<map->NbPixels();i++) {
   (*map)(i)=unvalmsk;
   if(!lat && !lon) continue;
   double theta,phi; map->PixThetaPhi(i,theta,phi);
   if(lat && (theta<lat2 || lat1<theta)) continue;
   if(lon && (  phi<lon1 || lon2<phi  )) continue;
   (*map)(i)=valmsk;
 }

 omg.AddObj(map,tokens[0]);
}

/* --Methode-- */
void skymapmoduleExecutor::CrMaskFrMap(vector<string>& tokens)
{
 NamedObjMgr omg;

 int_4 nside = atoi(tokens[1].c_str());
 if(nside<=1) return;
 if(!IsNSideGood(nside)) {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}

 AnyDataObj* obj=omg.GetObj(tokens[2]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[2]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[2]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 SphereHEALPix<r_8>* msk = new SphereHEALPix<r_8>(nside);

 double v1=-1.e100, v2=1.e100;
 if(tokens.size()>3) if(tokens[3][0]!='!')
   sscanf(tokens[3].c_str(),"%lf,%lf",&v1,&v2);

 double valmsk=0., unvalmsk=1.;
 if(tokens.size()>4) sscanf(tokens[4].c_str(),"%lf,%lf",&valmsk,&unvalmsk);

 cout<<"CrMaskFrMap "<<tokens[0]<<" nside"<<nside<<" with "<<tokens[2]<<endl
     <<"   for value in ["<<v1<<","<<v2<<"]"
     <<" valmsk="<<valmsk<<" unvalmsk="<<unvalmsk<<endl;

 for(int_4 i=0;i<msk->NbPixels();i++) {
   double theta,phi; msk->PixThetaPhi(i,theta,phi);
   int_4 ip = map->PixIndexSph(theta,phi);
   if(ip<0 || ip>=map->NbPixels()) continue;
   double v = (*map)(ip);
   if(v>v1 && v<v2) (*msk)(i)=valmsk;
     else           (*msk)(i)=unvalmsk;
 }

 omg.AddObj(msk,tokens[0]);
}

/* --Methode-- */
void skymapmoduleExecutor::MaskMap(vector<string>& tokens)
{
 NamedObjMgr omg;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 AnyDataObj* objm=omg.GetObj(tokens[1]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
   return;
 }
 SphereHEALPix<r_8>* msk = dynamic_cast<SphereHEALPix<r_8>*>(objm);
 if(msk==NULL) {
   cout<<"Error: "<<tokens[1]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 cout<<"MskMap: "<<tokens[0]<<" with mask "<<tokens[1]<<endl;

 for(int_4 i=0;i<map->NbPixels();i++) {
   double theta,phi; map->PixThetaPhi(i,theta,phi);
   int_4 ip = msk->PixIndexSph(theta,phi);
   if(ip<0 || ip>=msk->NbPixels()) continue;
   (*map)(i) *= (*msk)(ip);
 }

}

/* --Methode-- */
void skymapmoduleExecutor::MapProj(vector<string>& tokens)
{
 NamedObjMgr omg;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }
 int_4 nsidemap = map->SizeIndex();

 int_4 nside = nsidemap;
 if(tokens.size()>2) nside = atoi(tokens[2].c_str());
 if(!IsNSideGood(nside)) {cout<<"MapProj: Bad nside "<<nside<<endl; return;}

 SphereHEALPix<r_8>* mapproj = new SphereHEALPix<r_8>(nside);

 if(nsidemap==nside) {  // tailles egales
   for(int_4 i=0;i<mapproj->NbPixels();i++)
      (*mapproj)(i) = (*map)(i);
  } else if(nsidemap<nside) {  // mapproj>map
   for(int_4 i=0;i<mapproj->NbPixels();i++) {
     double theta,phi; mapproj->PixThetaPhi(i,theta,phi);
     int_4 ip = map->PixIndexSph(theta,phi);
     if(ip<0 || ip>=map->NbPixels()) continue;
     (*mapproj)(i) = (*map)(ip);
   }
 } else {  // mapproj<map
   SphereHEALPix<int_4> nproj(nside);
   nproj = 0; *mapproj = 0.;
   for(int_4 i=0;i<map->NbPixels();i++) {
     double theta,phi; map->PixThetaPhi(i,theta,phi);
     int_4 ip = mapproj->PixIndexSph(theta,phi);
     if(ip<0 || ip>=mapproj->NbPixels()) continue;
     (*mapproj)(ip) += (*map)(i);
     nproj(ip)++;
   }
   for(int_4 i=0;i<mapproj->NbPixels();i++)
      (*mapproj)(i) /= (r_8) nproj(i);
 }
 
 omg.AddObj(mapproj,tokens[1]);
}

/* --Methode-- */
void skymapmoduleExecutor::Map2Local(vector<string>& tokens)
{
 NamedObjMgr omg;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 int_4 nx=10,ny=10;
 if(tokens.size()>2) sscanf(tokens[2].c_str(),"%d,%d",&nx,&ny);
 if(nx<=0)
   {cout<<"Error: nx="<<nx<<" value >0"<<endl; return;}
 if(ny<=0)
   {cout<<"Error: ny="<<ny<<" value >0"<<endl; return;}

 double angleX=1., angleY=1.;
 if(tokens.size()>3) sscanf(tokens[3].c_str(),"%lf,%lf",&angleX,&angleY);
 if(angleX<=0. || angleX>180.)
   {cout<<"Error: bad angleX="<<angleX<<" value #]0,180]"<<endl; return;}
 if(angleY<=0. || angleY>180.)
   {cout<<"Error: bad angleY="<<angleY<<" value #]0,180]"<<endl; return;}

 double phi0=0., theta0=0.;
 if(tokens.size()>4) sscanf(tokens[4].c_str(),"%lf,%lf",&phi0,&theta0);
 if(phi0<0. || phi0>=360.)
   {cout<<"Error: bad phi0="<<phi0<<" value #[0,360["<<endl; return;}
 if(theta0<-90. || theta0>90.)
   {cout<<"Error: bad theta0="<<theta0<<" value #[-90,90]"<<endl; return;}

 int_4 x0=nx/2, y0=ny/2;
 if(tokens.size()>5) if(tokens[5][0]!='!')
   sscanf(tokens[5].c_str(),"%d,%d",&x0,&y0);

 double angle=0.;
 if(tokens.size()>6) sscanf(tokens[6].c_str(),"%lf",&angle);

 cout<<"Map2Local: proj "<<tokens[0]<<" to local "<<tokens[1]<<endl
     <<" nx="<<nx<<" ny="<<ny<<" angleX="<<angleX<<" angleY="<<angleY
     <<" phi0="<<phi0<<" theta0="<<theta0<<" x0="<<x0<<" y0="<<y0
     <<" angle="<<angle<<endl;
 if(angle<-90. || angle>90.)
   {cout<<"Error: bad angle="<<angle<<" value #[-90,90]"<<endl; return;}

 theta0 = (90.-theta0);
 LocalMap<r_8>* lmap = new LocalMap<r_8>(nx,ny,angleX,angleY,theta0,phi0,x0,y0,angle);
 *lmap = 0.;  //lmap->print(cout);

 int_4 nbad=0;
 double tmax=-9999., pmax=-9999., tmin=9999., pmin=9999.;
 for(int_4 i=0;i<lmap->NbPixels();i++) {
   double theta,phi;
   lmap->PixThetaPhi(i,theta,phi);
   int_4 ip = map->PixIndexSph(theta,phi);
   if(ip<0 || ip>=map->NbPixels()) {nbad++; continue;}
   if(theta<tmin) tmin=theta; if(theta>tmax) tmax=theta; 
   if(phi<pmin) pmin=phi; if(phi>pmax) pmax=phi; 
   (*lmap)(i) = (*map)(ip);
 }
 if(nbad) cout<<"Warning: "<<nbad<<" bad pixels numbers"<<endl;
 cout<<" phi=["<<pmin/M_PI*180.<<","<<pmax/M_PI*180.<<"]="
               <<(pmax-pmin)/M_PI*180.
     <<" theta=["<<tmin/M_PI*180.<<","<<tmax/M_PI*180.<<"]="
               <<(tmax-tmin)/M_PI*180.<<endl;

 omg.AddObj(lmap,tokens[1]); 
}

/* --Methode-- */
void skymapmoduleExecutor::MapOper(vector<string>& tokens)
{
 NamedObjMgr omg;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 cout<<"MapOper: "<<tokens[0];
 for(unsigned short iop=1;iop<tokens.size()-1;iop+=2) {
   char op = tokens[iop][0];
   if(op!='+' && op!='-' && op!='*' && op!='/') continue;
   AnyDataObj* obj1=omg.GetObj(tokens[iop+1]);
   if(obj1==NULL)
     {cout<<"Error: Pas d'objet de nom "<<tokens[iop+1]<<endl;
      continue;}
   SphereHEALPix<r_8>* map1 = dynamic_cast<SphereHEALPix<r_8>*>(obj1);
   if(map1==NULL)
     {cout<<"Error: "<<tokens[iop+1]<<" is not a SphereHEALPix<r_8>"<<endl;
      continue;}
   cout<<" "<<op<<"= "<<tokens[iop+1];
   for(int_4 i=0;i<map->NbPixels();i++) {
     double theta,phi; map->PixThetaPhi(i,theta,phi);
     int_4 ip = map1->PixIndexSph(theta,phi);
     if(ip<0 || ip>=map1->NbPixels()) continue;
     if(op=='+')      (*map)(i) += (*map1)(ip);
     else if(op=='-') (*map)(i) -= (*map1)(ip);
     else if(op=='*') (*map)(i) *= (*map1)(ip);
     else if(op=='/')
       {if((*map1)(ip)!=0.) (*map)(i) /= (*map1)(ip);
        else (*map)(i) = 0.;}
   }
 }
 cout<<endl;
}

/* --Methode-- */
void skymapmoduleExecutor::MapStat(vector<string>& tokens)
{
 NamedObjMgr omg;

 AnyDataObj* obj=omg.GetObj(tokens[0]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
   return;
 }
 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 SphereHEALPix<r_8>* msk = NULL;
 if(tokens.size()>1) {
   if(tokens[1][0]!='!') {
     AnyDataObj* objm=omg.GetObj(tokens[1]);
     if(objm==NULL) {
       cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
       return;
     }
     msk = dynamic_cast<SphereHEALPix<r_8>*>(objm);
     if(msk==NULL) {
       cout<<"Error: "<<tokens[1]<<" is not a SphereHEALPix<r_8>"<<endl;
       return;
     }
   }
 }

 cout<<"MapStat for map "<<tokens[0];
 if(msk) cout<<" using mask "<<tokens[1];
 cout<<endl;

 double sum=0., sum2=0., sumw=0.;
 int npix = map->NbPixels(), npixgood=0;
 for(int_4 i=0;i<npix;i++) {
   double w = 1.;
   if(msk) {
     double theta,phi; map->PixThetaPhi(i,theta,phi);
     int_4 ip = msk->PixIndexSph(theta,phi);
     if(ip<0 || ip>=msk->NbPixels()) w=0.;
       else                          w=(*msk)(ip);
   }
   if(w<=0.) continue;
   npixgood++; sumw += w;
   sum += (*map)(i)*w;
   sum2 += (*map)(i)*(*map)(i)*w*w;
 }
 if(sumw<=0. || npixgood==0) {
   sum = sum2 = sumw=0.;
 } else {
   sum /= sumw;
   sum2 = sum2/sumw - sum*sum;
   if(sum2<=0.) sum2=0.; else sum2=sqrt(sum2);
 }
 cout<<"Number of good px="<<npixgood<<" / "<<npix
     <<" ("<<100.*npixgood/npix<<" %)"<<endl;
 cout<<" mean="<<sum<<" sigma="<<sum2<<endl;

 if(tokens.size()>2) if(tokens[2][0]!='!') {
   if(omg.HasVar(tokens[2])) omg.DeleteVar(tokens[2]);
   char str[128]; sprintf(str,"%.8f",sum);
   omg.SetVar(tokens[2],(string)str);
   cout<<" mean stored into $"<<tokens[2];
 }
 if(tokens.size()>3) if(tokens[3][0]!='!') {
   if(omg.HasVar(tokens[3])) omg.DeleteVar(tokens[3]);
   char str[128]; sprintf(str,"%.8f",sum2);
   omg.SetVar(tokens[3],(string)str);
   cout<<" sigma stored into $"<<tokens[3];
 }
 if(tokens.size()>2) cout<<endl;

}

////////////////////////////////////////////////////////////
static skymapmoduleExecutor * piaskmex = NULL;

void skymapmodule_init()
{
if (piaskmex) delete piaskmex;
piaskmex = new skymapmoduleExecutor;
}

void skymapmodule_end()
{
// Desactivation du module
if (piaskmex) delete piaskmex;
piaskmex = NULL;
}

////////////////////////////////////////////////////////////