Ignore:
Timestamp:
Aug 4, 2002, 11:28:16 AM (23 years ago)
Author:
cmv
Message:

add methods cmv 4/8/2002

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/SophyaPI/ProgPI/skymapmodule.cc

    r2155 r2156  
    3232  void Map2Double(vector<string>& tokens);
    3333  void MapMult(vector<string>& tokens);
     34  void MapCover(vector<string>& tokens);
    3435  void Map2Cl(vector<string>& tokens);
     36  void Map2Alm(vector<string>& tokens);
    3537  void CrMapMask(vector<string>& tokens);
    3638  void CrMaskFrMap(vector<string>& tokens);
     
    5961mpiac->RegisterCommand(kw, usage, this, hgrp);
    6062
     63kw = "mapcover";
     64usage = "Print the covering of a map";
     65usage += "\n Usage: mapcover map v1 v2 [varname]";
     66usage += "\n    v1 v2: good pixels have v1<=val<=v2";
     67usage += "\n    $varname: percent of sky covering";
     68mpiac->RegisterCommand(kw, usage, this, hgrp);
     69
    6170kw = "map2cl";
    6271usage = "SkyMap to Cl";
    6372usage += "\n Usage: map2cl map(double) clvec [lmax] [niter]";
     73mpiac->RegisterCommand(kw, usage, this, hgrp);
     74
     75kw = "map2alm";
     76usage = "SkyMap to Alm";
     77usage += "\n Usage: map2alm map(double) almmat [lmax] [niter]";
    6478mpiac->RegisterCommand(kw, usage, this, hgrp);
    6579
     
    123137  }
    124138  MapMult(tokens);
     139} else if(kw == "mapcover") {
     140  if(tokens.size()<1) {
     141    cout<<"Usage: mapcover map v1 v2 [varname]"<<endl;
     142    return(0);
     143  }
     144  MapCover(tokens);
    125145} else if(kw == "map2cl") {
    126146  if(tokens.size()<2) {
     
    129149  }
    130150  Map2Cl(tokens);
     151} else if(kw == "map2alm") {
     152  if(tokens.size()<2) {
     153    cout<<"Usage: map2alm map(double) almmat [lmax] [niter]"<<endl;
     154    return(0);
     155  }
     156  Map2Alm(tokens);
    131157} else if(kw == "crmapmask") {
    132158  if(tokens.size()<2) {
     
    205231
    206232/* --Methode-- */
     233void skymapmoduleExecutor::MapCover(vector<string>& tokens)
     234{
     235 NamedObjMgr omg;
     236 AnyDataObj* obj=omg.GetObj(tokens[0]);
     237 if(obj==NULL) {
     238   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
     239   return;
     240 }
     241 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
     242 if(map==NULL) {
     243   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
     244   return;
     245 }
     246
     247 double v1=0.99, v2=1.01;
     248 if(tokens.size()>1) v1=atof(tokens[1].c_str());
     249 if(tokens.size()>2) v2=atof(tokens[2].c_str());
     250 cout<<"MapCover: "<<tokens[0]<<" good px=["<<v1<<","<<v2<<"]"<<endl;
     251 int_4 ngood=0;
     252 for(int_4 i=0;i<map->NbPixels();i++)
     253   if(v1<=(*map)(i) && (*map)(i)<=v2) ngood++;
     254 double per = (double)ngood/(double)map->NbPixels();
     255 cout<<"NGood="<<ngood<<" NTot="<<map->NbPixels()
     256     <<" -> "<<100.*per<<" %";
     257 if(tokens.size()>3) {
     258   if(omg.HasVar(tokens[3])) omg.DeleteVar(tokens[3]);
     259   char str[128]; sprintf(str,"%g",per);
     260   omg.SetVar(tokens[3],(string)str);
     261   cout<<" stored into $"<<tokens[3];
     262 }
     263 cout<<endl;
     264}
     265
     266/* --Methode-- */
    207267void skymapmoduleExecutor::Map2Cl(vector<string>& tokens)
    208268{
     
    238298 TVector<r_8>* cl = new TVector<r_8>(cltmp,false);
    239299 omg.AddObj(cl,tokens[1]);
     300}
     301
     302/* --Methode-- */
     303void skymapmoduleExecutor::Map2Alm(vector<string>& tokens)
     304{
     305 NamedObjMgr omg;
     306
     307 int_4 lmax=0, niter=3;
     308 if(tokens.size()>2) lmax = atoi(tokens[2].c_str());
     309 if(tokens.size()>3) niter = atoi(tokens[3].c_str());
     310 if(niter<=0) niter = 3;
     311
     312 AnyDataObj* obj=omg.GetObj(tokens[0]);
     313 if(obj==NULL) {
     314   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
     315   return;
     316 }
     317 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
     318 if(map==NULL) {
     319   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
     320   return;
     321 }
     322
     323 SphericalTransformServer<r_8> ylmserver;
     324
     325 int_4 nside = map->SizeIndex();
     326 if(lmax<=0) lmax = 3*nside-1;
     327
     328 cout<<"Map2Alm: "<<tokens[0]<<" -> "<<tokens[1]
     329     <<" lmax="<<lmax<<" niter="<<niter<<endl;
     330
     331 SphereHEALPix<r_8> tmpmap(*map,false);
     332 Alm<r_8> almtmp;
     333 ylmserver.DecomposeToAlm(tmpmap,almtmp,lmax,0.,niter);
     334
     335 int n = almtmp.rowNumber();
     336 TMatrix< complex<r_8> >* alm = new TMatrix< complex<r_8> >(n,n);
     337 *alm = (complex<r_8>)0.;
     338 for(int i=0;i<n;i++) for(int j=0;j<=i;j++) (*alm)(i,j)=almtmp(i,j);
     339 omg.AddObj(alm,tokens[1]);
    240340}
    241341
Note: See TracChangeset for help on using the changeset viewer.