Index: trunk/SophyaPI/ProgPI/skymapmodule.cc
===================================================================
--- trunk/SophyaPI/ProgPI/skymapmodule.cc	(revision 2161)
+++ trunk/SophyaPI/ProgPI/skymapmodule.cc	(revision 2162)
@@ -19,4 +19,10 @@
 #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" {
@@ -34,10 +40,14 @@
   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);
 };
 
@@ -63,6 +73,7 @@
 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";
+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);
@@ -71,4 +82,11 @@
 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);
 
@@ -76,20 +94,29 @@
 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";
-usage += "\n Usage: crmapmask mapmsk nside lat1 lat2 [valmsk=0] [valnomsk=1]";
+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    valmsk=value to be written into masked pixels";
-usage += "\n    valnomsk=value to be written into unmasked pixels";
+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=0] [valnomsk=1]";
+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";
-usage += "\n    valnomsk=value to be written into unmasked pixels";
+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);
 
@@ -107,4 +134,16 @@
 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";
@@ -114,4 +153,16 @@
 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);
+
 }
 
@@ -139,5 +190,5 @@
 } else if(kw == "mapcover") {
   if(tokens.size()<1) {
-    cout<<"Usage: mapcover map v1 v2 [varname]"<<endl;
+    cout<<"Usage: mapcover map v1,v2 [varname]"<<endl;
     return(0);
   }
@@ -149,4 +200,10 @@
   }
   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) {
@@ -155,7 +212,13 @@
   }
   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 [valmsk=0] [valnomsk=1]"<<endl;
+    cout<<"Usage: crmapmask mapmsk nside lat1,lat2 lon1,lon2 [valmsk,valnomsk]"<<endl;
     return(0);
   }
@@ -163,5 +226,5 @@
 } else if(kw == "crmaskfrmap") {
   if(tokens.size()<3) {
-    cout<<"Usage: crmaskfrmap mapmsk nside map(double) v1 v2 [valmsk=0] [valnomsk=1]"<<endl;
+    cout<<"Usage: crmaskfrmap mapmsk nside map(double) v1,v2 [valmsk,valnomsk]"<<endl;
     return(0);
   }
@@ -179,4 +242,10 @@
   }
   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) {
@@ -185,4 +254,10 @@
   }
   MapOper(tokens);
+} else if(kw == "mapstat") {
+  if(tokens.size()<1) {
+    cout<<"Usage: Usage: mapstat map [msk] [meanvarname] [sigmavarname]"<<endl;
+    return(0);
+  }
+  MapStat(tokens);
 }
 
@@ -246,6 +321,6 @@
 
  double v1=0.99, v2=1.01;
- if(tokens.size()>1) v1=atof(tokens[1].c_str());
- if(tokens.size()>2) v2=atof(tokens[2].c_str());
+ 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;
@@ -255,9 +330,9 @@
  cout<<"NGood="<<ngood<<" NTot="<<map->NbPixels()
      <<" -> "<<100.*per<<" %";
- if(tokens.size()>3) {
-   if(omg.HasVar(tokens[3])) omg.DeleteVar(tokens[3]);
+ if(tokens.size()>2) {
+   if(omg.HasVar(tokens[2])) omg.DeleteVar(tokens[2]);
    char str[128]; sprintf(str,"%g",per);
-   omg.SetVar(tokens[3],(string)str);
-   cout<<" stored into $"<<tokens[3];
+   omg.SetVar(tokens[2],(string)str);
+   cout<<" stored into $"<<tokens[2];
  }
  cout<<endl;
@@ -270,5 +345,6 @@
 
  int_4 lmax=0, niter=3;
- if(tokens.size()>2) lmax = atoi(tokens[2].c_str());
+ 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;
@@ -301,4 +377,35 @@
 
 /* --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)
 {
@@ -306,5 +413,6 @@
 
  int_4 lmax=0, niter=3;
- if(tokens.size()>2) lmax = atoi(tokens[2].c_str());
+ 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;
@@ -341,5 +449,78 @@
 
 /* --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;
@@ -347,64 +528,25 @@
  int_4 nside = atoi(tokens[1].c_str());
  if(nside<=1) return;
- int_4 n=nside; bool ok=true;
- while(n>1) {if(n%2!=0) {ok=false; break;} else n/=2;}
- if(!ok) {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}
-
- SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>(nside);
-
- double lat1=-1.e30, lat2=1.e30;
- if(tokens.size()>2) lat1=atof(tokens[2].c_str());
- if(tokens.size()>3) lat2=atof(tokens[3].c_str());
+ 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) valmsk=atof(tokens[4].c_str());
- if(tokens.size()>5) unvalmsk=atof(tokens[5].c_str());
-
- cout<<"CrMapMask "<<tokens[0]<<" nside="<<nside
-     <<" for latitude in ["<<lat1<<","<<lat2<<"]"
-     <<" valmsk="<<valmsk<<" unvalmsk="<<unvalmsk<<endl;
-
- lat1 = (90.-lat1)*M_PI/180.; lat2 = (90.-lat2)*M_PI/180.;
- for(int_4 i=0;i<map->NbPixels();i++) {
-   double theta,phi;
-   map->PixThetaPhi(i,theta,phi);
-   if(theta<lat1 && theta>lat2) (*map)(i)=valmsk;
-     else                       (*map)(i)=unvalmsk;
- }
-
- 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;
- int_4 n=nside; bool ok=true;
- while(n>1) {if(n%2!=0) {ok=false; break;} else n/=2;}
- if(!ok) {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}
-
- AnyDataObj* obj=omg.GetObj(tokens[2]);
- 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[1]<<" is not a SphereHEALPix<r_8>"<<endl;
-   return;
- }
-
- SphereHEALPix<r_8>* msk = new SphereHEALPix<r_8>(nside);
-
- double v1=-1.e50, v2=1.e50;
- if(tokens.size()>3) v1 = atof(tokens[3].c_str());
- if(tokens.size()>4) v2 = atof(tokens[4].c_str());
-
- double valmsk=0., unvalmsk=1.;
- if(tokens.size()>5) valmsk=atof(tokens[5].c_str());
- if(tokens.size()>6) unvalmsk=atof(tokens[6].c_str());
+ if(tokens.size()>4) sscanf(tokens[4].c_str(),"%lf,%lf",&valmsk,&unvalmsk);
 
  cout<<"CrMaskFrMap "<<tokens[0]<<" nside"<<nside<<" with "<<tokens[2]<<endl
@@ -413,6 +555,5 @@
 
  for(int_4 i=0;i<msk->NbPixels();i++) {
-   double theta,phi;
-   msk->PixThetaPhi(i,theta,phi);
+   double theta,phi; msk->PixThetaPhi(i,theta,phi);
    int_4 ip = map->PixIndexSph(theta,phi);
    if(ip<0 || ip>=map->NbPixels()) continue;
@@ -455,6 +596,5 @@
 
  for(int_4 i=0;i<map->NbPixels();i++) {
-   double theta,phi;
-   map->PixThetaPhi(i,theta,phi);
+   double theta,phi; map->PixThetaPhi(i,theta,phi);
    int_4 ip = msk->PixIndexSph(theta,phi);
    if(ip<0 || ip>=msk->NbPixels()) continue;
@@ -483,8 +623,5 @@
  int_4 nside = nsidemap;
  if(tokens.size()>2) nside = atoi(tokens[2].c_str());
- if(nside<=1) return;
- int_4 n=nside; bool ok=true;
- while(n>1) {if(n%2!=0) {ok=false; break;} else n/=2;}
- if(!ok) {cout<<"MapProj: Bad nside "<<nside<<endl; return;}
+ if(!IsNSideGood(nside)) {cout<<"MapProj: Bad nside "<<nside<<endl; return;}
 
  SphereHEALPix<r_8>* mapproj = new SphereHEALPix<r_8>(nside);
@@ -495,6 +632,5 @@
   } else if(nsidemap<nside) {  // mapproj>map
    for(int_4 i=0;i<mapproj->NbPixels();i++) {
-     double theta,phi;
-     mapproj->PixThetaPhi(i,theta,phi);
+     double theta,phi; mapproj->PixThetaPhi(i,theta,phi);
      int_4 ip = map->PixIndexSph(theta,phi);
      if(ip<0 || ip>=map->NbPixels()) continue;
@@ -505,6 +641,5 @@
    nproj = 0; *mapproj = 0.;
    for(int_4 i=0;i<map->NbPixels();i++) {
-     double theta,phi;
-     map->PixThetaPhi(i,theta,phi);
+     double theta,phi; map->PixThetaPhi(i,theta,phi);
      int_4 ip = mapproj->PixIndexSph(theta,phi);
      if(ip<0 || ip>=mapproj->NbPixels()) continue;
@@ -517,4 +652,79 @@
  
  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]); 
 }
 
@@ -549,6 +759,5 @@
    cout<<" "<<op<<"= "<<tokens[iop+1];
    for(int_4 i=0;i<map->NbPixels();i++) {
-     double theta,phi;
-     map->PixThetaPhi(i,theta,phi);
+     double theta,phi; map->PixThetaPhi(i,theta,phi);
      int_4 ip = map1->PixIndexSph(theta,phi);
      if(ip<0 || ip>=map1->NbPixels()) continue;
@@ -564,4 +773,82 @@
 }
 
+/* --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;
