source: Sophya/trunk/SophyaPI/ProgPI/skymapmodule.cc@ 3416

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

#include <typeinfo>
#include <vector>
#include <string>

#include "piacmd.h"
#include "nobjmgr.h"
#include "pistdimgapp.h"
#include "servnobjm.h"
#include "pidrawer.h"
#include "nomgadapter.h"

#include "tvector.h"
#include "ntuple.h"
#include "slinparbuff.h"
#include "localmap.h"
#include "spherehealpix.h"
#include "spherethetaphi.h"
#include "sphericaltransformserver.h"

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

class skymapmoduleExecutor : public CmdExecutor {
public:
  enum {UnKnown=(uint_2) 0,
        TypeR8=(uint_2) 16, TypeR4=(uint_2) 32,
        Spherical=(uint_2) 1, HealPix=(uint_2) 2, ThetaPhi=(uint_2) 4,
        Spherical8=(uint_2) Spherical|TypeR8,
        Spherical4=(uint_2) Spherical|TypeR4,
        HealPix8  =(uint_2) Spherical|HealPix|TypeR8,
        HealPix4  =(uint_2) Spherical|HealPix|TypeR4,
        ThetaPhi8 =(uint_2) Spherical|ThetaPhi|TypeR8,
        ThetaPhi4 =(uint_2) Spherical|ThetaPhi|TypeR4};

                skymapmoduleExecutor();
  virtual       ~skymapmoduleExecutor();
  virtual int   Execute(string& keyw, vector<string>& args, string& toks);

  void Resol2SizeIndex(vector<string>& tokens);
  void SizeIndex2Resol(vector<string>& tokens);
  void TypeMap(vector<string>& tokens);
  void Map2Double(vector<string>& tokens);
  void Map2Float(vector<string>& tokens);
  void Map2Map(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);
  void Alm2Cl(vector<string>& tokens);
  void Cl2llCl(vector<string>& tokens);
  void ClMean(vector<string>& tokens);
  void ClMult(vector<string>& tokens);
  void ClOper(vector<string>& tokens);
  void ClRebin(vector<string>& tokens);
protected:
  void SetTypeMap(vector<string>& tokens);
  bool IsNSideGood(int_4 nside);
  void GetNSideGood(int_4 &nside);
  uint_2 typemap(AnyDataObj* obj);

  uint_2 DefTypMap;
};

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

string hgrp = "SkyMapCmd";
string gdesc = "4 Pi spherical, local maps manipulation commands \n";
gdesc += "and spherical harmonics analysis Ylm. \n This is mainly "; 
gdesc += "Based on SOPHYA SkyMap and Samba modules. ";
mpiac->AddHelpGroup(hgrp, gdesc);
string kw,usage;

kw = "settypemap";
usage = "Set le type de map(double) par default";
usage += "\n Usage: settypemap type";
usage += "\n   type= H for Healpix (default)";
usage += "\n         T for ThetaPhi";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "resol2szidx";
usage = "Compute SizeIndex value (=nside for HEALPix) for a";
usage += "\n given resolution, (resol in arcminutes)";
usage += "\n Usage: resol2szidx resol";
mpiac->RegisterCommand(kw, usage, this, hgrp);
kw = "szidx2resol";
usage = "Compute resolution for a given SizeIndex (=nside for HEALPix)";
usage += "\n Usage: szidx2resol szidx_m";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "typemap";
usage = "Imprime le type de map";
usage += "\n Usage: typemap map";
mpiac->RegisterCommand(kw, usage, this, hgrp);

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

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

kw = "map2map";
usage = "Convert a map(double) into an other map type";
usage += "\n Usage: map2map map(double) type [nside]";
usage += "\n   type= H for Healpix";
usage += "\n         T for ThetaPhi";
usage += "\n   if nside <=1 use nside from map";
mpiac->RegisterCommand(kw, usage, this, hgrp);

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

kw = "mapcover";
usage = "Print the covering of a map(double)";
usage += "\n Usage: mapcover map(double) 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 map(double) 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 map(double) from Cl";
usage += "\n Usage: cl2map clvec map(double) nside [fwhm]";
usage += "\n    (see command settypemap)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "map2alm";
usage = "SkyMap map(double) 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 map(double) from Alm";
usage += "\n Usage: alm2map almmat map(double) nside";
usage += "\n    (see command settypemap)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "crmapmask";
usage = "Create a map mask(double) (nside) between latitudes lat1,lat2 longitudes lon1,lon2";
usage += "\n Usage: crmapmask mapmsk(double) 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)";
usage += "\n          (see command settypemap)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "crmaskfrmap";
usage = "Create a map mask(double) (nside) relative to an other map(double) pixels values";
usage += "\n Usage: crmaskfrmap mapmsk(double) 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(double) with a mask map(double)";
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(double) into an other one";
usage += "\n Usage: maproj map(double) mapproj(double) [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(double) into a local map(double)";
usage += "\n Usage: map2local map(double) localmap(double) 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(double)";
usage += "\n Usage: mapop map(double) op1 map1(double) [op2 map2(double)] [op2 map3(double)] [...]";
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(double)";
usage += "\n Usage: mapstat map(double) [msk(double)] [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);

kw = "alm2cl";
usage = "Project alm to cl";
usage += "\n Usage: alm2cl alm cl";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "cl2llcl";
usage = "Project cl to l*(l+1)*cl";
usage += "\n Usage: cl2llcl cl llcl";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "clmean";
usage = "Compute mean for a cl vector";
usage += "\n Usage: clmean cl imin,imax [meanvarname]";
usage += "\n   imin,imax  = compute between these indices";
usage += "\n   meanvarname  = variable name to store mean";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "clmult";
usage = "Multiply a cl by a constant";
usage += "\n Usage: clmult cl val";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "clop";
usage = "operation on a cl";
usage += "\n Usage: clop cl(double) op1 cl1(double) [op2 cl2(double)] [op2 cl3(double)] [...]";
usage += "\n   Do \"cl op1= cl1\", \"cl op2=cl2\", ...";
usage += "\n      where op = +,-,*,/";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "clrebin";
usage = "rebin a cl into an ntuple";
usage += "\n Usage: clrebin cl ntuple nbin,bin0";
usage += "\n   nbin: rebin by nbin";
usage += "\n   bin0: begin rebinning at bin bin0 (def=0)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

DefTypMap = HealPix;
}

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

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

if(kw == "settypemap") {
  SetTypeMap(tokens);
} else if(kw == "resol2szidx") {
  if(tokens.size()<1) {
    cout<<"Usage: resol2szidx resol"<<endl;
    return(0);
  }
  Resol2SizeIndex(tokens);
} else if(kw == "szidx2resol") {
  if(tokens.size()<1) {
    cout<<"Usage: szidx2resol szidx_m"<<endl;
    return(0);
  }
  SizeIndex2Resol(tokens);
} else if(kw == "typemap") {
  if(tokens.size()<1) {
    cout<<"Usage: typemap map"<<endl;
    return(0);
  }
  TypeMap(tokens);
} else if(kw == "map2double") {
  if(tokens.size()<1) {
    cout<<"Usage: map2double map"<<endl;
    return(0);
  }
  Map2Double(tokens);
} else if(kw == "map2float") {
  if(tokens.size()<1) {
    cout<<"Usage: map2float map"<<endl;
    return(0);
  }
  Map2Float(tokens);
} else if(kw == "map2map") {
  if(tokens.size()<2) {
    cout<<"Usage: map2map map type [nside]"<<endl;
    return(0);
  }
  Map2Map(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);
} else if(kw == "alm2cl") {
  if(tokens.size()<2) {
    cout<<"Usage: alm2cl alm cl"<<endl;
    return(0);
  }
  Alm2Cl(tokens);
} else if(kw == "cl2llcl") {
  if(tokens.size()<2) {
    cout<<"Usage: cl2llcl cl llcl"<<endl;
    return(0);
  }
  Cl2llCl(tokens);
} else if(kw == "clmean") {
  if(tokens.size()<1) {
    cout<<"Usage: clmean cl imin,imax [meanvarname]"<<endl;
    return(0);
  }
  ClMean(tokens);
} else if(kw == "clmult") {
  if(tokens.size()<1) {
    cout<<"Usage: clmult cl val"<<endl;
    return(0);
  }
  ClMult(tokens);
} else if(kw == "clop") {
  if(tokens.size()<3) {
    cout<<"Usage: clop cl op1 cl1 [op2 cl2] [op2 cl3] [...]"<<endl;
    return(0);
  }
  ClOper(tokens);
} else if(kw == "clrebin") {
  if(tokens.size()<2) {
    cout<<"Usage: clrebin cl ntuple nbin,bin0"<<endl;
    return(0);
  }
  ClRebin(tokens);
}

return(0);
}

/* --Methode-- */
void skymapmoduleExecutor::Resol2SizeIndex(vector<string>& tokens)
{
 double resamin = atof(tokens[0].c_str());
 double resrad = Angle(resamin,Angle::ArcMin).ToRadian();
 int_4 mh = SphereHEALPix<r_4>::ResolToSizeIndex(resrad);
 int_4 mt = SphereThetaPhi<r_4>::ResolToSizeIndex(resrad);
 cout << "Resol2SizeIndex: Resol= " << resamin << " arcmin ->" 
      << resrad << " radian  NPix ~= " << 4.*M_PI/(resrad*resrad) 
      << " \n  ===> nside=m_HEALPix= " << mh << "  m_ThetaPhi= " << mt << endl;
 return;
}

void skymapmoduleExecutor::SizeIndex2Resol(vector<string>& tokens)
{
 int_4 m = atoi(tokens[0].c_str());
 double resh = SphereHEALPix<r_4>::SizeIndexToResol(m);
 double rest = SphereThetaPhi<r_4>::SizeIndexToResol(m);
 cout << "SizeIndex2Resol: SizeIndex (m)= " << m
      << "\n HEALPix:  " << Angle(resh).ToArcMin() << " arcmin "
      << " NPix= " << 12*m*m 
      << "\n ThetaPhi: " << Angle(rest).ToArcMin() << " arcmin "
      << " NPix ~= " <<  (int_4)(4.*M_PI/(rest*rest)) << endl;
 return;
}

/* --Methode-- */
void skymapmoduleExecutor::TypeMap(vector<string>& tokens)
{
 NamedObjMgr omg;
 AnyDataObj* obj=omg.GetObj(tokens[0]);
 uint_2 t = typemap(obj);
 if(t==UnKnown)  {cout<<"TypeMap: UnKnown"<<endl; return;}

 int nside = 0;
 if(t&TypeR8) {
   SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
   nside = map->SizeIndex();
 } else if(t&TypeR4) {
   SphericalMap<r_4>* map = dynamic_cast<SphericalMap<r_4>*>(obj);
   nside = map->SizeIndex();
 }

 string dum = "";
 if(t==HealPix8)        dum = "SphereHEALPix<r_8>";
 else if(t==HealPix4)   dum = "SphereHEALPix<r_4>";
 else if(t==ThetaPhi8)  dum = "SphereThetaPhi<r_8>";
 else if(t==ThetaPhi4)  dum = "SphereThetaPhi<r_4>";
 else if(t==Spherical8) dum = "SphericalMap<r_8>";
 else if(t==Spherical4) dum = "SphericalMap<r_4>";

 cout<<"TypeMap: "<<dum<<" nside="<<nside<<endl;
}

/* --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;
 }
 SphericalMap<r_4>* map = dynamic_cast<SphericalMap<r_4>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_4>"<<endl;
   return;
 }
 uint_2 t = typemap(obj);
 int_4 nside = map->SizeIndex();

 SphericalMap<r_8>* map8 = NULL;
 if(t&HealPix)       map8 = new SphereHEALPix<r_8>(nside);
 else if(t&ThetaPhi) map8 = new SphereThetaPhi<r_8>(nside);
 else return;

 cout<<"Map2Double: converting to double "<<tokens[0]<<" nside="<<nside<<endl;
 for(int_4 i=0;i<map8->NbPixels();i++) (*map8)(i) = (r_8) (*map)(i);

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

/* --Methode-- */
void skymapmoduleExecutor::Map2Float(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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
   return;
 }
 uint_2 t = typemap(obj);
 int_4 nside = map->SizeIndex();

 SphericalMap<r_4>* map4 = NULL;
 if(t&HealPix)       map4 = new SphereHEALPix<r_4>(nside);
 else if(t&ThetaPhi) map4 = new SphereThetaPhi<r_4>(nside);
 else return;

 cout<<"Map2Float: converting to float "<<tokens[0]<<" nside="<<nside<<endl;
 for(int_4 i=0;i<map4->NbPixels();i++) (*map4)(i) = (r_4) (*map)(i);

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

/* --Methode-- */
void skymapmoduleExecutor::Map2Map(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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
   return;
 }
 uint_2 t = typemap(obj);
 int_4 nside = map->SizeIndex();

 uint_2 tomap = UnKnown;
 if(toupper(tokens[1][0])=='H')      tomap = HealPix;
 else if(toupper(tokens[1][0])=='T') tomap = ThetaPhi;
 else {cout<<"unkown map type "<<(char)toupper(tokens[1][0])<<" (H or T)!"<<endl; return;}
 if(t&tomap) {cout<<"map of same type, no conversion done"<<endl; return;}

 int_4 tonside = nside;
 if(tokens.size()>2) tonside = atoi(tokens[2].c_str());
 else {
   if(tomap&HealPix) {tonside=int(nside/1.5); GetNSideGood(tonside);}
   else tonside=int(map->NbThetaSlices()/2.5);
 }
 if(tomap&HealPix && !IsNSideGood(tonside))
   {cout<<"Bad nside for Healpix "<<tonside<<endl; return;}

 SphericalMap<r_8>* map8 = NULL;
 if(tomap&HealPix)       map8 = new SphereHEALPix<r_8>(tonside);
 else if(tomap&ThetaPhi) map8 = new SphereThetaPhi<r_8>(tonside);
 else return;

 cout<<"Map2Map: convert "<<tokens[0]<<" nside="<<nside
     <<" to type "<<tokens[1]<<" tonside="<<tonside
     <<"   ntheta="<<map->NbThetaSlices()<<" -> "<<map8->NbThetaSlices()<<endl;

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

 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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<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[512]; 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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
   return;
 }
 uint_2 t = typemap(obj);

 SphericalTransformServer<r_8> ylmserver;

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

 SphericalMap<r_8>* tmpmap = NULL;
 if(t&HealPix)       tmpmap = new SphereHEALPix<r_8>(*((SphereHEALPix<r_8>*)map),false);
 else if(t&ThetaPhi) tmpmap = new SphereThetaPhi<r_8>(*((SphereThetaPhi<r_8>*)map),false);
 else return;

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

 TVector<r_8> cltmp = ylmserver.DecomposeToCl(*tmpmap,lmax,0.,niter);
 cout<<"("<<cltmp.Size()<<")"<<endl;
 delete tmpmap;

 omg.AddObj(cltmp,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());
 if(DefTypMap&HealPix && !IsNSideGood(nside))
   {cout<<"Cl2Map: Bad nside for HealPix "<<nside<<endl; return;}

 double fwhm = 0.;
 if(tokens.size()>3) fwhm = atof(tokens[3].c_str());

 SphericalMap<r_8>* map = NULL;
 if(DefTypMap&HealPix)       map = new SphereHEALPix<r_8>;
 else if(DefTypMap&ThetaPhi) map = new SphereThetaPhi<r_8>;
 else return;

 SphericalTransformServer<r_8> ylmserver;
 cout<<"Cl2Map: Generate map "<<tokens[1]<<" with nside="<<nside
     <<" from cl "<<tokens[0]<<"("<<cl->Size()<<")"
     <<" 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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
   return;
 }
 uint_2 t = typemap(obj);

 SphericalTransformServer<r_8> ylmserver;

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

 SphericalMap<r_8>* tmpmap = NULL;
 if(t&HealPix)       tmpmap = new SphereHEALPix<r_8>(*((SphereHEALPix<r_8>*)map),false);
 else if(t&ThetaPhi) tmpmap = new SphereThetaPhi<r_8>(*((SphereThetaPhi<r_8>*)map),false);
 else return;

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

 Alm<r_8> almtmp;
 ylmserver.DecomposeToAlm(*tmpmap,almtmp,lmax,0.,niter);
 cout<<"("<<almtmp.rowNumber()<<")"<<endl;
 delete tmpmap;

 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());
 if(DefTypMap&HealPix && !IsNSideGood(nside))
   {cout<<"Alm2Map: Bad nside for HealPix "<<nside<<endl; return;}

 SphericalMap<r_8>* map = NULL;
 if(DefTypMap&HealPix)       map = new SphereHEALPix<r_8>;
 else if(DefTypMap&ThetaPhi) map = new SphereThetaPhi<r_8>;
 else return;

 SphericalTransformServer<r_8> ylmserver;
 cout<<"Alm2Map: Generate map "<<tokens[1]<<" with nside="<<nside
     <<" from alm "<<tokens[0]<<"("<<alm.rowNumber()<<")"<<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(DefTypMap&HealPix && !IsNSideGood(nside))
   {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}

 SphericalMap<r_8>* map = NULL;
 if(DefTypMap&HealPix)       map = new SphereHEALPix<r_8>(nside);
 else if(DefTypMap&ThetaPhi) map = new SphereThetaPhi<r_8>(nside);
 else return;

 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;

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

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

 SphericalMap<r_8>* msk = NULL;
 if(t&HealPix)       msk = new SphereHEALPix<r_8>(nside);
 else if(t&ThetaPhi) msk = new SphereThetaPhi<r_8>(nside);
 else return;

 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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
   return;
 }

 AnyDataObj* objm=omg.GetObj(tokens[1]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
   return;
 }
 SphericalMap<r_8>* msk = dynamic_cast<SphericalMap<r_8>*>(objm);
 if(msk==NULL) {
   cout<<"Error: "<<tokens[1]<<" is not a SphericalMap<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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
   return;
 }
 int_4 nsidemap = map->SizeIndex();
 uint_2 t = typemap(obj);

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

 SphericalMap<r_8>* mapproj = NULL;
 if(t&HealPix)       mapproj = new SphereHEALPix<r_8>(nside);
 else if(t&ThetaPhi) mapproj = new SphereThetaPhi<r_8>(nside);
 else return;

 cout<<"MapProj: project "<<tokens[0]<<" n="<<nsidemap
     <<" to "<<tokens[1]<<" n="<<nside<<endl;

 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
   SphericalMap<int_4>*  nproj= NULL;
   if(t&HealPix)       nproj= new SphereHEALPix<int_4>(nside);
   else if(t&ThetaPhi) nproj= new SphereThetaPhi<int_4>(nside);
   int_4 i;
   for(i=0;i<nproj->NbPixels();i++) {(*nproj)(i)=0;(*mapproj)(i)=0.;}
   for(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(i=0;i<mapproj->NbPixels();i++)
      (*mapproj)(i) /= (r_8) (*nproj)(i);
   delete nproj;
 }
 
 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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
   return;
 }

 cout<<"MapOper: "<<tokens[0];
 for(unsigned short iop=1;iop<tokens.size()-1;iop+=2) {
   string op = tokens[iop];
   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;}
   SphericalMap<r_8>* map1 = dynamic_cast<SphericalMap<r_8>*>(obj1);
   if(map1==NULL)
     {cout<<"Error: "<<tokens[iop+1]<<" is not a SphericalMap<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;
 }
 SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
   return;
 }

 SphericalMap<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<SphericalMap<r_8>*>(objm);
     if(msk==NULL) {
       cout<<"Error: "<<tokens[1]<<" is not a SphericalMap<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<<" sigma^2="<<sum2*sum2<<endl;

 if(tokens.size()>2) if(tokens[2][0]!='!') {
   if(omg.HasVar(tokens[2])) omg.DeleteVar(tokens[2]);
   char str[512]; 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[512]; sprintf(str,"%.8f",sum2);
   omg.SetVar(tokens[3],(string)str);
   cout<<" sigma stored into $"<<tokens[3];
 }
 if(tokens.size()>2) cout<<endl;

}

/* --Methode-- */
void skymapmoduleExecutor::Alm2Cl(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);

 TVector<r_8> cl(almmat->NRows());
 cl = 0.;

 cout<<"Alm2Cl: project alm "<<tokens[0]<<"("<<alm.rowNumber()
     <<") to cl "<<tokens[1]<<"("<<cl.Size()<<")"<<endl;

 for(int_4 l=0;l<alm.rowNumber();l++) {
   cl(l) = alm(l,0).real()*alm(l,0).real()+alm(l,0).imag()*alm(l,0).imag();
   if(l>0) for(int_4 m=1;m<=l;m++)
     cl(l) += 2.*(alm(l,m).real()*alm(l,m).real()+alm(l,m).imag()*alm(l,m).imag());
   cl(l) /= 2.*l+1;
 }

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

/* --Methode-- */
void skymapmoduleExecutor::Cl2llCl(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;
 }

 TVector<r_8> llcl(*cl,false);

 cout<<"Cl2llCl: project "<<tokens[0]<<"("<<cl->Size()
     <<") to l*(l+1)*cl "<<tokens[1]<<"("<<llcl.Size()<<")"<<endl;

 for(int_4 l=0;l<llcl.Size();l++) llcl(l) *= (double)l*(l+1.);

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

/* --Methode-- */
void skymapmoduleExecutor::ClMean(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 lmin=0, lmax=-1;
 if(tokens.size()>1) if(tokens[1][0]!='!')
   sscanf(tokens[1].c_str(),"%d,%d",&lmin,&lmax);
 if(lmin<0) lmin=0; if(lmax>=cl->Size()) lmax=cl->Size()-1;
 if(lmin>lmax) {lmin=0; lmax=cl->Size()-1;}

 cout<<"ClMean: compute mean "<<tokens[0]<<"("<<cl->Size()
     <<") between ["<<lmin<<","<<lmax<<"]";

 double mean = 0.;
 for(int_4 l=lmin;l<=lmax;l++) mean += (*cl)(l);
 mean /= lmax-lmin+1;

 cout<<" mean="<<mean;

 if(tokens.size()>2) {
   if(omg.HasVar(tokens[2])) omg.DeleteVar(tokens[2]);
   char str[512]; sprintf(str,"%.8f",mean);
   omg.SetVar(tokens[2],(string)str);
   cout<<" stored into $"<<tokens[2];
 }
 cout<<endl;

}

/* --Methode-- */
void skymapmoduleExecutor::ClMult(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;
 }

 double val = 1.;
 if(tokens.size()>1) val=atof(tokens[1].c_str());

 cout<<"ClMult: multiply "<<tokens[0]<<"("<<cl->Size()<<") by "<<val<<endl;

 *cl *= val;

}

/* --Methode-- */
void skymapmoduleExecutor::ClOper(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;
 }

 cout<<"ClOper: "<<tokens[0]<<"("<<cl->Size()<<")";
 for(unsigned short iop=1;iop<tokens.size()-1;iop+=2) {
   string op = tokens[iop];
   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;}
   TVector<r_8>* cl1 = dynamic_cast<TVector<r_8>*>(obj1);
   if(cl1==NULL)
     {cout<<"Error: "<<tokens[iop+1]<<" is not a TVector<r_8>"<<endl;
      continue;}
   cout<<" "<<op<<"= "<<tokens[iop+1]<<"("<<cl1->Size()<<")";
   int_4 n=cl->Size(); if(n>cl1->Size()) n=cl1->Size();
   for(int_4 i=0;i<n;i++) {
     if(op=="+")      (*cl)(i) += (*cl1)(i);
     else if(op=="-") (*cl)(i) -= (*cl1)(i);
     else if(op=="*") (*cl)(i) *= (*cl1)(i);
     else if(op=="/")
       {if((*cl1)(i)!=0.) (*cl)(i)/=(*cl1)(i); else (*cl)(i)=0.;}
   }
 }
 cout<<endl;
}
/* --Methode-- */
void skymapmoduleExecutor::ClRebin(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_4 nn = cl->Size();
 if(nn<=0) return;

 int_4 nrebin=1, ibin0=0;
 if(tokens.size()>2) sscanf(tokens[2].c_str(),"%d,%d",&nrebin,&ibin0);
 if(ibin0<0 || ibin0>=nn) ibin0=0;

 cout<<"ClRebin: rebin "<<tokens[0]<<"("<<nn<<") by nbin="<<nrebin
     <<" begin at "<<ibin0<<" -> "<<tokens[1]<<endl;
 
 const int ndim = 8;
 char *nament[ndim] =
   {"l","n","clmean","cllin","clpar","sclmean","scllin","sclpar"};
 NTuple clbin(ndim,nament);
 r_4 xnt[ndim];

 SLinParBuff slp(nrebin+2);

 for(int_4 i=ibin0;i<nn;i+=nrebin) {
   r_8 ll=0.;
   slp.Reset();
   int_4 j;
   for(j=0;j<ndim;j++) xnt[j]=0.;
   for(j=i;j<nn;j++) {
     ll += (r_8) j;
     slp.Push((r_8)j,(*cl)(j));
     if((int_4)slp.NPoints()>=nrebin) break;
   }
   if(slp.NPoints()<=0) continue;
   ll /=  (r_8) slp.NPoints();
   xnt[0] = ll;
   xnt[1] = slp.NPoints();
   r_8 s,a0,a1,a2;
   s = slp.Compute(a0);
   if(s>0.) {xnt[2]=a0; xnt[5]=s;}
   s = slp.Compute(a0,a1);
   if(s>0.) {xnt[3]=a0+a1*ll; xnt[6]=s;}
   s = slp.Compute(a0,a1,a2);
   if(s>0.) {xnt[4]=a0+a1*ll+a2*ll*ll; xnt[7]=s;}
   clbin.Fill(xnt);
 }

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

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

/* --Methode-- */
void skymapmoduleExecutor::SetTypeMap(vector<string>& tokens)
{
 if(tokens.size()<1) {
   cout<<"SetTypeMap: Usage: settypemap type"<<endl;
 } else {
   if(toupper(tokens[0][0])=='H')      DefTypMap = HealPix;
   else if(toupper(tokens[0][0])=='T') DefTypMap = ThetaPhi;
   else cout<<"unkown map type, should be (H or T)!"<<endl;
 }
 string dum;
 if(DefTypMap&HealPix)       dum="HealPix";
 else if(DefTypMap&ThetaPhi) dum="ThetaPhi";
 cout<<"SetTypeMap: type map is "<<dum<<" ("<<DefTypMap<<")"<<endl;
}

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

/* --Methode-- */
void skymapmoduleExecutor::GetNSideGood(int_4 &nside)
// get the nearest good nside for HealPix
{
 double n=1.e50; int_4 ns=nside;
 for(int_4 i=1;i<66000;i*=2) {
   double v = fabs((double)(nside-i));
   if(v>n) continue;
   n=v; ns=i;
 }
 nside = ns;
}

/* --Methode-- */
uint_2 skymapmoduleExecutor::typemap(AnyDataObj* obj)
{
  if(obj==NULL) return UnKnown;
  if(dynamic_cast<SphereHEALPix<r_4> *>(obj)) return HealPix4;
  if(dynamic_cast<SphereHEALPix<r_8> *>(obj)) return HealPix8;
  if(dynamic_cast<SphereThetaPhi<r_4>*>(obj)) return ThetaPhi4;
  if(dynamic_cast<SphereThetaPhi<r_8>*>(obj)) return ThetaPhi8;
  if(dynamic_cast<SphericalMap<r_4>  *>(obj)) return Spherical4;
  if(dynamic_cast<SphericalMap<r_8>  *>(obj)) return Spherical8;
  return UnKnown;
}

////////////////////////////////////////////////////////////
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;
}