Index: trunk/SophyaPI/ProgPI/skymapmodule.cc
===================================================================
--- trunk/SophyaPI/ProgPI/skymapmodule.cc	(revision 2162)
+++ trunk/SophyaPI/ProgPI/skymapmodule.cc	(revision 2175)
@@ -20,10 +20,4 @@
 #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();
@@ -33,8 +27,21 @@
 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 TypeMap(vector<string>& tokens);
   void Map2Double(vector<string>& tokens);
+  void Map2Map(vector<string>& tokens);
   void MapMult(vector<string>& tokens);
   void MapCover(vector<string>& tokens);
@@ -50,4 +57,11 @@
   void MapOper(vector<string>& tokens);
   void MapStat(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;
 };
 
@@ -60,8 +74,29 @@
 
 string hgrp = "SkyMapCmd";
-
-string kw = "map2double";
-string usage = "Convert a float map to a double map";
+string kw,usage;
+
+kw = "settypemap";
+usage = "Set le type de map par default";
+usage += "\n   type= H for Healpix (default)";
+usage += "\n         T for ThetaPhi";
+usage += "\n Usage: settypemap type";
+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";
+mpiac->RegisterCommand(kw, usage, this, hgrp);
+
+kw = "map2map";
+usage = "Convert a map into an other map type";
+usage += "\n   type= H for Healpix";
+usage += "\n         T for ThetaPhi";
+usage += "\n   if nside <=1 use nside from map";
+usage += "\n Usage: map2map map type [nside]";
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
@@ -89,4 +124,5 @@
 usage = "Generate SkyMap from Cl";
 usage += "\n Usage: cl2map clvec map nside [fwhm]";
+usage += "\n    (see command settypemap)";
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
@@ -101,4 +137,5 @@
 usage = "Generate SkyMap from Alm";
 usage += "\n Usage: alm2map almmat map nside";
+usage += "\n    (see command settypemap)";
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
@@ -111,4 +148,5 @@
 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);
 
@@ -165,4 +203,5 @@
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
+DefTypMap = HealPix;
 }
 
@@ -176,5 +215,13 @@
 {
 
-if(kw == "map2double") {
+if(kw == "settypemap") {
+  SetTypeMap(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;
@@ -182,4 +229,10 @@
   }
   Map2Double(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) {
@@ -266,4 +319,32 @@
 
 /* --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)
 {
@@ -274,17 +355,75 @@
    return;
  }
- SphereHEALPix<r_4>* map = dynamic_cast<SphereHEALPix<r_4>*>(obj);
+ SphericalMap<r_4>* map = dynamic_cast<SphericalMap<r_4>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_4>"<<endl;
-   return;
- }
+   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_4>"<<endl;
+   return;
+ }
+ uint_2 t = typemap(obj);
  int_4 nside = map->SizeIndex();
- SphereHEALPix<r_8>* map8 = new SphereHEALPix<r_8>(nside);
- cout<<"Map2Double: converting to double "<<tokens[0]<<endl;
+
+ 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::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)
 {
@@ -295,7 +434,7 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
+   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
    return;
  }
@@ -314,7 +453,7 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
+   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
    return;
  }
@@ -323,5 +462,7 @@
  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++)
@@ -355,9 +496,10 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
-   return;
- }
+   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
+   return;
+ }
+ uint_2 t = typemap(obj);
 
  SphericalTransformServer<r_8> ylmserver;
@@ -369,9 +511,13 @@
      <<" 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]);
+ 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;
+
+ TVector<r_8> cltmp = ylmserver.DecomposeToCl(*tmpmap,lmax,0.,niter);
+ delete tmpmap;
+
+ omg.AddObj(cltmp,tokens[1]);
 }
 
@@ -394,8 +540,14 @@
  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());
 
- SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>;
+ 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;
@@ -423,9 +575,10 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
-   return;
- }
+   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
+   return;
+ }
+ uint_2 t = typemap(obj);
 
  SphericalTransformServer<r_8> ylmserver;
@@ -437,7 +590,12 @@
      <<" lmax="<<lmax<<" niter="<<niter<<endl;
 
- SphereHEALPix<r_8> tmpmap(*map,false);
+ 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;
+
  Alm<r_8> almtmp;
- ylmserver.DecomposeToAlm(tmpmap,almtmp,lmax,0.,niter);
+ ylmserver.DecomposeToAlm(*tmpmap,almtmp,lmax,0.,niter);
+ delete tmpmap;
 
  int n = almtmp.rowNumber();
@@ -471,5 +629,11 @@
  int nside = 128;
  if(tokens.size()>2) nside = atoi(tokens[2].c_str());
- SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>;
+ 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;
@@ -487,7 +651,11 @@
 
  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);
+ 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;
@@ -526,20 +694,25 @@
  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(!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);
+ 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;
@@ -550,5 +723,5 @@
  if(tokens.size()>4) sscanf(tokens[4].c_str(),"%lf,%lf",&valmsk,&unvalmsk);
 
- cout<<"CrMaskFrMap "<<tokens[0]<<" nside"<<nside<<" with "<<tokens[2]<<endl
+ cout<<"CrMaskFrMap "<<tokens[0]<<" nside="<<nside<<" with "<<tokens[2]<<endl
      <<"   for value in ["<<v1<<","<<v2<<"]"
      <<" valmsk="<<valmsk<<" unvalmsk="<<unvalmsk<<endl;
@@ -576,7 +749,7 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
+   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
    return;
  }
@@ -587,7 +760,7 @@
    return;
  }
- SphereHEALPix<r_8>* msk = dynamic_cast<SphereHEALPix<r_8>*>(objm);
+ SphericalMap<r_8>* msk = dynamic_cast<SphericalMap<r_8>*>(objm);
  if(msk==NULL) {
-   cout<<"Error: "<<tokens[1]<<" is not a SphereHEALPix<r_8>"<<endl;
+   cout<<"Error: "<<tokens[1]<<" is not a SphericalMap<r_8>"<<endl;
    return;
  }
@@ -614,16 +787,24 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
+   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(!IsNSideGood(nside)) {cout<<"MapProj: Bad nside "<<nside<<endl; return;}
-
- SphereHEALPix<r_8>* mapproj = new SphereHEALPix<r_8>(nside);
+ 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
@@ -638,6 +819,8 @@
    }
  } else {  // mapproj<map
-   SphereHEALPix<int_4> nproj(nside);
-   nproj = 0; *mapproj = 0.;
+   SphericalMap<int_4>*  nproj= NULL;
+   if(t&HealPix)       nproj= new SphereHEALPix<int_4>(nside);
+   else if(t&ThetaPhi) nproj= new SphereThetaPhi<int_4>(nside);
+   for(int_4 i=0;i<nproj->NbPixels();i++) {(*nproj)(i)=0;(*mapproj)(i)=0.;}
    for(int_4 i=0;i<map->NbPixels();i++) {
      double theta,phi; map->PixThetaPhi(i,theta,phi);
@@ -645,8 +828,9 @@
      if(ip<0 || ip>=mapproj->NbPixels()) continue;
      (*mapproj)(ip) += (*map)(i);
-     nproj(ip)++;
+     (*nproj)(ip)++;
    }
    for(int_4 i=0;i<mapproj->NbPixels();i++)
-      (*mapproj)(i) /= (r_8) nproj(i);
+      (*mapproj)(i) /= (r_8) (*nproj)(i);
+   delete nproj;
  }
  
@@ -664,7 +848,7 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
+   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
    return;
  }
@@ -739,7 +923,7 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
+   cout<<"Error: "<<tokens[0]<<" is not a SphericalMap<r_8>"<<endl;
    return;
  }
@@ -753,7 +937,7 @@
      {cout<<"Error: Pas d'objet de nom "<<tokens[iop+1]<<endl;
       continue;}
-   SphereHEALPix<r_8>* map1 = dynamic_cast<SphereHEALPix<r_8>*>(obj1);
+   SphericalMap<r_8>* map1 = dynamic_cast<SphericalMap<r_8>*>(obj1);
    if(map1==NULL)
-     {cout<<"Error: "<<tokens[iop+1]<<" is not a SphereHEALPix<r_8>"<<endl;
+     {cout<<"Error: "<<tokens[iop+1]<<" is not a SphericalMap<r_8>"<<endl;
       continue;}
    cout<<" "<<op<<"= "<<tokens[iop+1];
@@ -783,11 +967,11 @@
    return;
  }
- SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ SphericalMap<r_8>* map = dynamic_cast<SphericalMap<r_8>*>(obj);
  if(map==NULL) {
-   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
-   return;
- }
-
- SphereHEALPix<r_8>* msk = 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]!='!') {
@@ -797,7 +981,7 @@
        return;
      }
-     msk = dynamic_cast<SphereHEALPix<r_8>*>(objm);
+     msk = dynamic_cast<SphericalMap<r_8>*>(objm);
      if(msk==NULL) {
-       cout<<"Error: "<<tokens[1]<<" is not a SphereHEALPix<r_8>"<<endl;
+       cout<<"Error: "<<tokens[1]<<" is not a SphericalMap<r_8>"<<endl;
        return;
      }
@@ -852,4 +1036,56 @@
 
 ////////////////////////////////////////////////////////////
+
+/* --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;
 
@@ -866,4 +1102,2 @@
 piaskmex = NULL;
 }
-
-////////////////////////////////////////////////////////////
Index: trunk/SophyaPI/ProgPI/skymapmodule.pic
===================================================================
--- trunk/SophyaPI/ProgPI/skymapmodule.pic	(revision 2175)
+++ trunk/SophyaPI/ProgPI/skymapmodule.pic	(revision 2175)
@@ -0,0 +1,272 @@
+#----------------------
+shell rm -f h4*.fits t4*.ppf h8*.fits t8*.ppf clh*.ppf clt*.ppf almh*.ppf almt*.ppf
+delobjs *
+c++import SkyMap
+c++exec \
+  double theta,phi; \
+  SphereHEALPix<r_4>  h4(256); \
+  for(int_4 i=0;i<h4.NbPixels();i++) { \
+    h4.PixThetaPhi(i,theta,phi); \
+    h4(i) = 0.; \
+    double v = (M_PI_2-theta)*(M_PI_2-theta)/(0.50*0.50); \
+    v += (phi-M_PI)*(phi-M_PI)/(M_PI*M_PI/2.); \
+    v = 100.*exp(-v/2.); \
+    v += 10.*NorRand(); \
+    if(phi<20.*M_PI/180.||phi>340.*M_PI/180.) continue; \
+    if(phi<190.*M_PI/180.&&phi>170.*M_PI/180.) continue; \
+    h4(i) = v; \
+  } \
+  KeepObj(h4); \
+  SphereThetaPhi<r_4> t4(400); \
+  for(int_4 i=0;i<t4.NbPixels();i++) { \
+    t4.PixThetaPhi(i,theta,phi); \
+    t4(i) = 0.; \
+    double v = (M_PI_2-theta)*(M_PI_2-theta)/(0.50*0.50); \
+    v += (phi-M_PI)*(phi-M_PI)/(M_PI*M_PI/2.); \
+    v = 100.*exp(-v/2.); \
+    v += 10.*NorRand(); \
+    if(phi<20.*M_PI/180.||phi>340.*M_PI/180.) continue; \
+    if(phi<190.*M_PI/180.&&phi>170.*M_PI/180.) continue; \
+    t4(i) = v; \
+  } \
+  KeepObj(t4); \
+  cout<<"End of Job"<<endl;
+
+savefits h4 !h4.fits
+saveppf t4 t4.ppf
+
+zone 1 2
+disp h4 zoom/4
+disp t4 zoom/4
+zone
+
+#----------------------
+delobjs *
+openfits h4.fits
+
+map2double h4
+rename h4 h8
+savefits h8 !h8.fits
+
+delobjs *
+openppf t4.ppf
+map2double t4
+rename t4 t8
+saveppf t8 t8.ppf
+
+#----------------------
+delobjs *
+openfits h4.fits
+openfits h8.fits
+openppf t4.ppf
+openppf t8.ppf
+typemap h4
+typemap h8
+typemap t4
+typemap t8
+
+#----------------------
+delobjs *
+openfits h8.fits
+openppf t8.ppf
+map2map h8 t
+map2map t8 h
+
+#----------------------
+delobjs *
+openfits h8.fits
+openppf t8.ppf
+mapmult h8 1000
+mapmult t8 1000
+
+#----------------------
+delobjs *
+openfits h8.fits
+maproj h8 h8p
+mapop h8p - h8
+disp h8p zoom/4
+
+maproj h8 h8p64 64
+disp h8p64 zoom/4
+savefits h8p64 !h8p64.fits
+
+delobjs *
+openppf t8.ppf
+maproj t8 t8p
+mapop t8p - t8
+disp t8p zoom/4
+
+maproj t8 t8p100 100
+disp t8p100 zoom/4
+saveppf t8p100 t8p100.ppf
+
+#----------------------
+delobjs *
+openfits h8p64.fits
+map2cl h8p64 clh64
+saveppf clh64 clh64.ppf
+n/plot clh64.n*(n+1)*val%n ! ! crossmarker5
+
+map2cl h8p64 clh64 50
+n/plot clh64.n*(n+1)*val%n ! ! crossmarker5
+
+delobjs *
+openppf t8p100.ppf
+map2cl t8p100 clt100 191
+saveppf clt100 clt100.ppf
+n/plot clt100.n*(n+1)*val%n ! ! crossmarker5
+
+#----------------------
+delobjs *
+openfits h8p64.fits
+map2alm h8p64 almh64
+saveppf almh64 almh64.ppf
+disp almh64
+
+delobjs *
+openppf t8p100.ppf
+map2alm t8p100 almt100 191
+saveppf almt100 almt100.ppf
+disp almt100
+
+#----------------------
+delobjs *
+openppf clh64.ppf
+settypemap h
+cl2map clh64 h8p64frcl 64
+savefits h8p64frcl !h8p64frcl.fits
+disp h8p64frcl
+
+delobjs *
+openppf clt100.ppf
+settypemap t
+cl2map clt100 t8p100frcl 100
+saveppf t8p100frcl t8p100frcl.ppf
+disp t8p100frcl
+
+#----------------------
+delobjs *
+openppf almh64.ppf
+settypemap h
+alm2map almh64 h8p64fralm 64
+savefits h8p64fralm !h8p64fralm.fits
+disp h8p64fralm
+
+delobjs *
+openppf almt100.ppf
+settypemap t
+alm2map almt100 t8p100fralm 100
+saveppf t8p100fralm t8p100fralm.ppf
+disp t8p100fralm
+
+#----------------------
+delobjs *
+settypemap h
+crmapmask h8m 256 -20,20 100,130 0,1
+savefits h8m !h8m.fits
+disp h8m zoom/4
+
+delobjs *
+settypemap t
+crmapmask t8m 400 -20,20 100,130 0,1
+saveppf t8m t8m.ppf
+disp t8m zoom/4
+
+#----------------------
+delobjs *
+openfits h8.fits
+crmaskfrmap h8fm 256 h8 -1.e-30,1.e-30 0,1
+savefits h8fm !h8fm.fits
+disp h8fm zoom/4
+
+delobjs *
+openppf t8.ppf
+crmaskfrmap t8fm 400 t8 -1.e-30,1.e-30 0,1
+saveppf t8fm t8fm.ppf
+disp t8fm zoom/4
+
+#----------------------
+delobjs *
+openfits h8.fits
+openfits h8m.fits
+maskmap h8 h8m
+disp h8 zoom/4
+
+delobjs *
+openppf t8.ppf
+openppf t8m.ppf
+maskmap t8 t8m
+disp t8 zoom/4
+
+delobjs *
+openfits h8.fits
+openppf t8m.ppf
+maskmap h8 t8m
+disp h8 zoom/4
+
+delobjs *
+openppf t8.ppf
+openfits h8m.fits
+maskmap t8 h8m
+disp t8 zoom/4
+
+#----------------------
+delobjs *
+openfits h8m.fits
+mapcover h8m 0.9,1 couvh8m
+echo $couvh8m
+
+delobjs *
+openppf t8m.ppf
+mapcover t8m 0.9,1 couvt8m
+echo $couvt8m
+
+#----------------------
+delobjs *
+openfits h8.fits
+map2local h8 h8loc 200,300 20,30 0,90 !
+disp h8loc
+map2local h8 h8loc 200,300 20,30 0,90 ! 45
+disp h8loc
+
+delobjs *
+openppf t8.ppf
+map2local t8 t8loc 200,300 20,30 0,90 !
+disp t8loc
+map2local t8 t8loc 200,300 20,30 0,90 ! 45
+disp t8loc
+
+#----------------------
+delobjs *
+openfits h8.fits
+rename h8 h8save
+openfits h8.fits
+mapop h8 + h8save - h8save * h8save / h8save
+disp h8 zoom/4
+
+openppf t8.ppf
+rename t8 t8save
+openppf t8.ppf
+mapop t8 + t8save - t8save * t8save / t8save
+disp t8 zoom/4
+
+mapop h8 - t8save
+disp h8 zoom/4
+mapop t8 - h8save
+disp t8 zoom/4
+
+#----------------------
+delobjs *
+openfits h8.fits
+openfits h8fm.fits
+mapstat h8 !    mean sig
+echo mean=$mean sig=$sig
+mapstat h8 h8fm mean sig
+echo mean=$mean sig=$sig
+
+openppf t8.ppf
+openppf t8fm.ppf
+mapstat t8 !    mean sig
+echo mean=$mean sig=$sig
+mapstat t8 t8fm mean sig
+echo mean=$mean sig=$sig
