Changeset 2177 in Sophya for trunk/SophyaPI

Timestamp:

Aug 11, 2002, 4:53:54 PM (23 years ago)

Author:

cmv

Message:

+ de fonctionnalites cmv 11/8/02

Location:

trunk/SophyaPI/ProgPI

Files:

• skymapmodule.cc (modified) (27 diffs)
• skymapmodule.pic (modified) (29 diffs)

trunk/SophyaPI/ProgPI/skymapmodule.cc

@@ -13 +13 @@
 #include "pistdimgapp.h"
 #include "servnobjm.h"
-#include "tvector.h"
 #include "pidrawer.h"
 #include "nomgadapter.h"
 
-#include "skymap.h"
+#include "tvector.h"
+#include "ntuple.h"
+#include "slinparbuff.h"
+#include "localmap.h"
+#include "spherehealpix.h"
+#include "spherethetaphi.h"
 #include "sphericaltransformserver.h"
 
@@ -43 +47 @@
   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);
@@ -57 +62 @@
   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);
+  void VarOper(vector<string>& tokens);
 protected:
   void SetTypeMap(vector<string>& tokens);
@@ -77 +89 @@
 
 kw = "settypemap";
-usage = "Set le type de map par default";
+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";
-usage += "\n Usage: settypemap type";
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
@@ -90 +102 @@
 kw = "map2double";
 usage = "Convert a float map to a double map";
-usage += "\n Usage: map2double 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 into an other map type";
+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";
-usage += "\n Usage: map2map map type [nside]";
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
 kw = "mapmult";
-usage = "Multiply a map by a constant";
-usage += "\n Usage: mapmult map cste";
+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";
-usage += "\n Usage: mapcover map v1,v2 [varname]";
+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): \"!\")";
@@ -115 +132 @@
 
 kw = "map2cl";
-usage = "SkyMap to Cl";
+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)";
@@ -122 +139 @@
 
 kw = "cl2map";
-usage = "Generate SkyMap from Cl";
-usage += "\n Usage: cl2map clvec map nside [fwhm]";
+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 to Alm";
+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)";
@@ -135 +152 @@
 
 kw = "alm2map";
-usage = "Generate SkyMap from Alm";
-usage += "\n Usage: alm2map almmat map nside";
+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 (nside) between latitudes lat1,lat2 longitudes lon1,lon2";
-usage += "\n Usage: crmapmask mapmsk nside lat1,lat2 lon1,lon2 [valmsk,valnomsk]";
+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)";
@@ -152 +169 @@
 
 kw = "crmaskfrmap";
-usage = "Create a map mask (nside) relative to an other map pixels values";
-usage += "\n Usage: crmaskfrmap mapmsk nside map(double) v1,v2 [valmsk,valnomsk]";
+usage = "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)";
@@ -160 +177 @@
 
 kw = "maskmap";
-usage = "Mask a map with a mask map";
+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)";
@@ -166 +183 @@
 
 kw = "maproj";
-usage = "Project a map into an other one";
-usage += "\n Usage: maproj map mapproj [nside]";
+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)";
@@ -173 +190 @@
 
 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 = "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";
@@ -185 +202 @@
 
 kw = "mapop";
-usage = "operation on a map";
-usage += "\n Usage: mapop map op1 map1 [op2 map2] [op2 map3] [...]";
+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 = +,-,*,/";
@@ -192 +209 @@
 
 kw = "mapstat";
-usage = "Statistics from a map";
-usage += "\n Usage: mapstat map [msk] [meanvarname] [sigmavarname]";
+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";
@@ -203 +220 @@
 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);
+
+kw = "varop";
+usage = "operation on variables";
+usage += "\n Usage: varop var1 op var2 [resultvarname]";
+usage += "\n   Do result = var1 op var2";
+usage += "\n      where op = +,-,*,/,^(pow),exp,sqrt,log,log10";
+usage += "\n   resultvarname = variable name to store the result";
+mpiac->RegisterCommand(kw, usage, this, hgrp);
+
 DefTypMap = HealPix;
 }
@@ -229 +290 @@
   }
   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) {
@@ -313 +380 @@
   }
   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);
+} else if(kw == "varop") {
+  if(tokens.size()<2) {
+    cout<<"Usage: varop var1 op var2 [resultvarname]"<<endl;
+    return(0);
+  }
+  VarOper(tokens);
 }
 
@@ -372 +481 @@
 
  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]);
 }
 
@@ -508 +645 @@
  if(lmax<=0) lmax = 3*nside-1;
 
- cout<<"Map2Cl: "<<tokens[0]<<" -> "<<tokens[1]
-     <<" lmax="<<lmax<<" niter="<<niter<<endl;
-
  SphericalMap<r_8>* tmpmap = NULL;
  if(t&HealPix)       tmpmap = new SphereHEALPix<r_8>(*((SphereHEALPix<r_8>*)map),false);
@@ -516 +650 @@
  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;
 
@@ -553 +691 @@
  SphericalTransformServer<r_8> ylmserver;
  cout<<"Cl2Map: Generate map "<<tokens[1]<<" with nside="<<nside
-     <<" from cl "<<tokens[0]<<" with fwhm="<<fwhm<<endl;
+     <<" from cl "<<tokens[0]<<"("<<cl->Size()<<")"
+     <<" with fwhm="<<fwhm<<endl;
  ylmserver.GenerateFromCl(*map,nside,*cl,fwhm);
 
@@ -587 +726 @@
  if(lmax<=0) lmax = 3*nside-1;
 
- cout<<"Map2Alm: "<<tokens[0]<<" -> "<<tokens[1]
-     <<" lmax="<<lmax<<" niter="<<niter<<endl;
-
  SphericalMap<r_8>* tmpmap = NULL;
  if(t&HealPix)       tmpmap = new SphereHEALPix<r_8>(*((SphereHEALPix<r_8>*)map),false);
@@ -595 +731 @@
  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;
 
@@ -639 +779 @@
  SphericalTransformServer<r_8> ylmserver;
  cout<<"Alm2Map: Generate map "<<tokens[1]<<" with nside="<<nside
-     <<" from alm "<<tokens[0]<<endl;
+     <<" from alm "<<tokens[0]<<"("<<alm.rowNumber()<<")"<<endl;
  ylmserver.GenerateFromAlm(*map,nside,alm);
 
@@ -931 +1071 @@
  cout<<"MapOper: "<<tokens[0];
  for(unsigned short iop=1;iop<tokens.size()-1;iop+=2) {
-   char op = tokens[iop][0];
-   if(op!='+' && op!='-' && op!='*' && op!='/') continue;
+   string op = tokens[iop];
+   if(op!="+" && op!="-" && op!="*" && op!="/") continue;
    AnyDataObj* obj1=omg.GetObj(tokens[iop+1]);
    if(obj1==NULL)
@@ -946 +1086 @@
      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.;}
+     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.;}
    }
  }
@@ -1035 +1174 @@
 }
 
+/* --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[128]; 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();
+   for(int_4 j=0;j<ndim;j++) xnt[j]=0.;
+   for(int_4 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;
+   xnt[0] = ll / (r_8) slp.NPoints();
+   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::VarOper(vector<string>& tokens)
+{
+ NamedObjMgr omg;
+
+ double var1=0., var2=0.; string op = "?";
+ if(tokens.size()>0) var1 = atof(tokens[0].c_str());
+ if(tokens.size()>1) op = tokens[1];
+ if(tokens.size()>2) if(tokens[2][0]!='!') var2 = atof(tokens[2].c_str());
+
+ double result = 0.;
+ try {
+   if(op=="+")                     result = var1 + var2;
+   else if(op=="-")                result = var1 - var2;
+   else if(op=="*")                result = var1 * var2;
+   else if(op=="/" && var2!=0.)    result = var1 / var2;   else if(op=="^" && var2>0.)     result = pow(var1,var2);
+   else if(op=="pow" && var2>0.)   result = pow(var1,var2);
+   else if(op=="exp")              result = exp(var1);
+   else if(op=="sqrt" && var1>0.)  result = sqrt(var1);
+   else if(op=="log" && var1>0.)   result = log(var1);
+   else if(op=="log10" && var1>0.) result = log10(var1);
+ } catch ( ... ) {
+   cout<<"An operation exception occurs"<<endl;
+   return;
+ }
+
+ cout<<"VarOper: "<<var1<<" "<<op<<" "<<var2<<" = "<<result;
+
+ if(tokens.size()>3) {
+   if(omg.HasVar(tokens[3])) omg.DeleteVar(tokens[3]);
+   char str[128]; sprintf(str,"%.8f",result);
+   omg.SetVar(tokens[3],(string)str);
+   cout<<" stored into $"<<tokens[3];
+ }
+ cout<<endl;
+
+}
+
 ////////////////////////////////////////////////////////////
 

trunk/SophyaPI/ProgPI/skymapmodule.pic

@@ -1 @@
-#----------------------
+#---------------------- Create Maps and Save
 shell rm -f h4*.fits t4*.ppf h8*.fits t8*.ppf clh*.ppf clt*.ppf almh*.ppf almt*.ppf
 delobjs *
@@ -39 +39 @@
 disp h4 zoom/4
 disp t4 zoom/4
-zone
-
-#----------------------
-delobjs *
+
+#---------------------- test Map2Double
+delobjs *
+zone
 openfits h4.fits
 
@@ -50 +50 @@
 
 delobjs *
+zone
 openppf t4.ppf
 map2double t4
@@ -55 +56 @@
 saveppf t8 t8.ppf
 
-#----------------------
-delobjs *
+#---------------------- test Map2Float
+delobjs *
+zone
+openfits h8.fits
+openppf t8.ppf
+
+map2float h8
+map2float t8
+typemap h8
+typemap t8
+
+#---------------------- test TypeMap
+delobjs *
+zone
 openfits h4.fits
 openfits h8.fits
@@ -66 +79 @@
 typemap t8
 
-#----------------------
-delobjs *
+#---------------------- test Map2Map
+delobjs *
+zone
 openfits h8.fits
 openppf t8.ppf
@@ -73 +87 @@
 map2map t8 h
 
-#----------------------
-delobjs *
+#---------------------- test MapMult
+delobjs *
+zone
 openfits h8.fits
 openppf t8.ppf
@@ -80 +95 @@
 mapmult t8 1000
 
-#----------------------
-delobjs *
+#---------------------- test MapProj
+delobjs *
+zone
 openfits h8.fits
 maproj h8 h8p
@@ -92 +108 @@
 
 delobjs *
+zone
 openppf t8.ppf
 maproj t8 t8p
@@ -101 +118 @@
 saveppf t8p100 t8p100.ppf
 
-#----------------------
-delobjs *
+#---------------------- test Map2Cl
+delobjs *
+zone
 openfits h8p64.fits
 map2cl h8p64 clh64
@@ -112 +130 @@
 
 delobjs *
+zone
 openppf t8p100.ppf
 map2cl t8p100 clt100 191
@@ -117 +136 @@
 n/plot clt100.n*(n+1)*val%n ! ! crossmarker5
 
-#----------------------
-delobjs *
+#---------------------- test Map2Alm
+delobjs *
+zone
 openfits h8p64.fits
 map2alm h8p64 almh64
@@ -125 +145 @@
 
 delobjs *
+zone
 openppf t8p100.ppf
 map2alm t8p100 almt100 191
@@ -130 +151 @@
 disp almt100
 
-#----------------------
-delobjs *
+#---------------------- test Alm2Cl
+delobjs *
+zone
+openppf clh64.ppf
+openppf almh64.ppf
+alm2cl almh64 clh64fralm
+n/plot clh64.n*(n+1)*val%n ! ! crossmarker5
+n/plot clh64fralm.n*(n+1)*val%n ! ! "red same circlemarker5"
+c++exec clh64fralm -= clh64;
+n/plot clh64fralm.n*(n+1)*val%n ! ! crossmarker5
+
+delobjs *
+zone
+openppf clt100.ppf
+openppf almt100.ppf
+alm2cl almt100 clt100fralm
+n/plot clt100.n*(n+1)*val%n ! ! crossmarker5
+n/plot clt100fralm.n*(n+1)*val%n ! ! "red same circlemarker5"
+c++exec clt100fralm -= clt100;
+n/plot clt100fralm.n*(n+1)*val%n ! ! crossmarker5
+
+#---------------------- test Cl2Map et SetTypeMap
+delobjs *
+zone
 openppf clh64.ppf
 settypemap h
@@ -139 +182 @@
 
 delobjs *
+zone
 openppf clt100.ppf
 settypemap t
@@ -145 +189 @@
 disp t8p100frcl
 
-#----------------------
-delobjs *
+#---------------------- test Alm2Map et SetTypeMap
+delobjs *
+zone
 openppf almh64.ppf
 settypemap h
@@ -154 +199 @@
 
 delobjs *
+zone
 openppf almt100.ppf
 settypemap t
@@ -160 +206 @@
 disp t8p100fralm
 
-#----------------------
-delobjs *
+#---------------------- test Cl2llCl
+delobjs *
+zone
+openppf clh64.ppf
+cl2llcl clh64 llclh64
+n/plot clh64.n*(n+1)*val%n ! ! crossmarker5
+n/plot llclh64.val%n ! ! "red same circlemarker5"
+
+delobjs *
+zone
+openppf clt100.ppf
+cl2llcl clt100 llclt100
+n/plot clt100.n*(n+1)*val%n ! ! crossmarker5
+n/plot llclt100.val%n ! ! "red same circlemarker5"
+
+#---------------------- test ClMean
+delobjs *
+zone
+openppf clh64.ppf
+cl2llcl clh64 llclh64
+n/plot llclh64.val%n ! ! circlemarker5
+clmean llclh64
+clmean llclh64 !
+clmean llclh64 100,9999 mean
+echo $mean
+
+delobjs *
+zone
+openppf clt100.ppf
+cl2llcl clt100 llclt100
+n/plot llclt100.val%n ! ! circlemarker5
+clmean llclt100
+clmean llclt100 !
+clmean llclt100 100,9999 mean
+echo $mean
+
+#---------------------- test ClMult
+delobjs *
+openppf clh64.ppf
+zone 1 2
+n/plot clh64.val%n ! ! circlemarker5
+clmult clh64 1000.
+n/plot clh64.val%n ! ! circlemarker5
+
+delobjs *
+openppf clt100.ppf
+zone 1 2
+n/plot clt100.val%n ! ! circlemarker5
+clmult clt100 1000.
+n/plot clt100.val%n ! ! circlemarker5
+
+#---------------------- test ClOper
+delobjs *
+zone
+openppf clh64.ppf
+cp clh64 dum
+clop dum + clh64 - clh64 * clh64 / clh64
+n/plot clh64.val%n ! ! crossmarker5
+n/plot dum.val%n ! ! "same red circlemarker5"
+
+delobjs *
+zone
+openppf clt100.ppf
+cp clt100 dum
+clop dum + clt100 - clt100 * clt100 / clt100
+n/plot clt100.val%n ! ! crossmarker5
+n/plot dum.val%n ! ! "same red circlemarker5"
+
+delobjs *
+zone
+openppf clh64.ppf
+openppf clt100.ppf
+clop clh64 - clt100
+n/plot clh64.val%n ! ! crossmarker5
+
+#---------------------- test ClRebin
+delobjs *
+zone
+openppf clh64.ppf
+cl2llcl clh64 llclh64
+clrebin llclh64 clntu 10,0
+
+n/plot clntu.n%l ! ! crossmarker5
+
+n/plot llclh64.val%n  ! ! "crossmarker3"
+n/plot clntu.clmean%l ! ! "same circlemarker5 red"
+n/plot clntu.cllin%l  ! ! "same boxmarker5 blue"
+n/plot clntu.clpar%l  ! ! "same trianglemarker5 violet"
+
+n/plot llclh64.val%n  ! ! "crossmarker3"
+nt2d clntu l clmean 0 sclmean 1 " " "same circlemarker5 red thinline"
+
+n/plot llclh64.val%n  ! ! "crossmarker3"
+nt2d clntu l cllin 0 scllin 1 " " "same boxmarker5 blue thinline"
+
+n/plot llclh64.val%n  ! ! "crossmarker3"
+nt2d clntu l clpar 0 sclpar 1 " " "same trianglemarker5 violet thinline"
+
+#---------------------- test CrMaskMap et SetTypeMap
+delobjs *
+zone
 settypemap h
 crmapmask h8m 256 -20,20 100,130 0,1
@@ -168 +313 @@
 
 delobjs *
+zone
 settypemap t
 crmapmask t8m 400 -20,20 100,130 0,1
@@ -173 +319 @@
 disp t8m zoom/4
 
-#----------------------
-delobjs *
+#---------------------- test CrMaskFrMap
+delobjs *
+zone
 openfits h8.fits
 crmaskfrmap h8fm 256 h8 -1.e-30,1.e-30 0,1
@@ -181 +328 @@
 
 delobjs *
+zone
 openppf t8.ppf
 crmaskfrmap t8fm 400 t8 -1.e-30,1.e-30 0,1
@@ -186 +334 @@
 disp t8fm zoom/4
 
-#----------------------
-delobjs *
+#---------------------- test MaskMap
+delobjs *
+zone
 openfits h8.fits
 openfits h8m.fits
@@ -194 +343 @@
 
 delobjs *
+zone
 openppf t8.ppf
 openppf t8m.ppf
@@ -200 +350 @@
 
 delobjs *
+zone
 openfits h8.fits
 openppf t8m.ppf
@@ -206 +357 @@
 
 delobjs *
+zone
 openppf t8.ppf
 openfits h8m.fits
@@ -211 +363 @@
 disp t8 zoom/4
 
-#----------------------
-delobjs *
+#---------------------- test MapCover
+delobjs *
+zone
 openfits h8m.fits
 mapcover h8m 0.9,1 couvh8m
@@ -218 +371 @@
 
 delobjs *
+zone
 openppf t8m.ppf
 mapcover t8m 0.9,1 couvt8m
 echo $couvt8m
 
-#----------------------
-delobjs *
+#---------------------- test Map2Local
+delobjs *
+zone
 openfits h8.fits
 map2local h8 h8loc 200,300 20,30 0,90 !
@@ -231 +386 @@
 
 delobjs *
+zone
 openppf t8.ppf
 map2local t8 t8loc 200,300 20,30 0,90 !
@@ -237 +393 @@
 disp t8loc
 
-#----------------------
-delobjs *
-openfits h8.fits
-rename h8 h8save
-openfits h8.fits
-mapop h8 + h8save - h8save * h8save / h8save
+#---------------------- test MapOper
+delobjs *
+zone
+openfits h8.fits
+cp h8 h8save
+mapop h8save + h8 - h8 * h8 / h8
+disp h8save zoom/4
+mapop h8save - h8
+disp h8save zoom/4
+
+delobjs *
+zone
+openppf t8.ppf
+cp t8 t8save
+mapop t8save + t8 - t8 * t8 / t8
+disp t8save zoom/4
+mapop t8save - t8
+disp t8save zoom/4
+
+delobjs *
+zone
+openfits h8.fits
+openppf t8.ppf
+mapop h8 - t8
 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 *
+#---------------------- test MapStat
+delobjs *
+zone
 openfits h8.fits
 openfits h8fm.fits
@@ -271 +435 @@
 mapstat t8 t8fm mean sig
 echo mean=$mean sig=$sig
+
+#---------------------- test VarOper
+varop 1000 + 100 result
+  echo result is $result
+varop 1000 - 100 result
+  echo result is $result
+varop 1000 * 100 result
+  echo result is $result
+varop 1000 / 100 result
+  echo result is $result
+varop 1000 / 0 result
+  echo result is $result
+
+varop 1000 sqrt ! result
+  echo result is $result
+varop -1000 sqrt ! result
+  echo result is $result
+varop 10 exp ! result
+  echo result is $result
+varop 1000 log ! result
+  echo result is $result
+varop -1000 log ! result
+  echo result is $result
+varop 1000 log10 ! result
+  echo result is $result
+varop -1000 log10 ! result
+  echo result is $result
+
+varop 1000 pow 2 result
+  echo result is $result
+varop 1000 pow 1000 result
+  echo result is $result