Changeset 2156 in Sophya for trunk/SophyaPI/ProgPI/skymapmodule.cc

Timestamp:

Aug 4, 2002, 11:28:16 AM (23 years ago)

Author:

cmv

Message:

add methods cmv 4/8/2002

File:

• trunk/SophyaPI/ProgPI/skymapmodule.cc (modified) (6 diffs)

trunk/SophyaPI/ProgPI/skymapmodule.cc

@@ -32 +32 @@
   void Map2Double(vector<string>& tokens);
   void MapMult(vector<string>& tokens);
+  void MapCover(vector<string>& tokens);
   void Map2Cl(vector<string>& tokens);
+  void Map2Alm(vector<string>& tokens);
   void CrMapMask(vector<string>& tokens);
   void CrMaskFrMap(vector<string>& tokens);
@@ -59 +61 @@
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
+kw = "mapcover";
+usage = "Print the covering of a map";
+usage += "\n Usage: mapcover map v1 v2 [varname]";
+usage += "\n    v1 v2: good pixels have v1<=val<=v2";
+usage += "\n    $varname: percent of sky covering";
+mpiac->RegisterCommand(kw, usage, this, hgrp);
+
 kw = "map2cl";
 usage = "SkyMap to Cl";
 usage += "\n Usage: map2cl map(double) clvec [lmax] [niter]";
+mpiac->RegisterCommand(kw, usage, this, hgrp);
+
+kw = "map2alm";
+usage = "SkyMap to Alm";
+usage += "\n Usage: map2alm map(double) almmat [lmax] [niter]";
 mpiac->RegisterCommand(kw, usage, this, hgrp);
 
@@ -123 +137 @@
   }
   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) {
@@ -129 +149 @@
   }
   Map2Cl(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 == "crmapmask") {
   if(tokens.size()<2) {
@@ -205 +231 @@
 
 /* --Methode-- */
+void skymapmoduleExecutor::MapCover(vector<string>& tokens)
+{
+ NamedObjMgr omg;
+ AnyDataObj* obj=omg.GetObj(tokens[0]);
+ if(obj==NULL) {
+   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
+   return;
+ }
+ SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ if(map==NULL) {
+   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
+   return;
+ }
+
+ double v1=0.99, v2=1.01;
+ if(tokens.size()>1) v1=atof(tokens[1].c_str());
+ if(tokens.size()>2) v2=atof(tokens[2].c_str());
+ 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()>3) {
+   if(omg.HasVar(tokens[3])) omg.DeleteVar(tokens[3]);
+   char str[128]; sprintf(str,"%g",per);
+   omg.SetVar(tokens[3],(string)str);
+   cout<<" stored into $"<<tokens[3];
+ }
+ cout<<endl;
+}
+
+/* --Methode-- */
 void skymapmoduleExecutor::Map2Cl(vector<string>& tokens)
 {
@@ -238 +298 @@
  TVector<r_8>* cl = new TVector<r_8>(cltmp,false);
  omg.AddObj(cl,tokens[1]);
+}
+
+/* --Methode-- */
+void skymapmoduleExecutor::Map2Alm(vector<string>& tokens)
+{
+ NamedObjMgr omg;
+
+ int_4 lmax=0, niter=3;
+ if(tokens.size()>2) lmax = atoi(tokens[2].c_str());
+ if(tokens.size()>3) niter = atoi(tokens[3].c_str());
+ if(niter<=0) niter = 3;
+
+ AnyDataObj* obj=omg.GetObj(tokens[0]);
+ if(obj==NULL) {
+   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
+   return;
+ }
+ SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
+ if(map==NULL) {
+   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
+   return;
+ }
+
+ SphericalTransformServer<r_8> ylmserver;
+
+ int_4 nside = map->SizeIndex();
+ if(lmax<=0) lmax = 3*nside-1;
+
+ cout<<"Map2Alm: "<<tokens[0]<<" -> "<<tokens[1]
+     <<" lmax="<<lmax<<" niter="<<niter<<endl;
+
+ SphereHEALPix<r_8> tmpmap(*map,false);
+ Alm<r_8> almtmp;
+ ylmserver.DecomposeToAlm(tmpmap,almtmp,lmax,0.,niter);
+
+ int n = almtmp.rowNumber();
+ TMatrix< complex<r_8> >* alm = new TMatrix< complex<r_8> >(n,n);
+ *alm = (complex<r_8>)0.;
+ for(int i=0;i<n;i++) for(int j=0;j<=i;j++) (*alm)(i,j)=almtmp(i,j);
+ omg.AddObj(alm,tokens[1]);
 }