source: Sophya/trunk/SophyaPI/ProgPI/skymapmodule.cc@ 2158

#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>

#include <typeinfo>
#include <vector>
#include <string>

#include "piacmd.h"
#include "nobjmgr.h"
#include "pistdimgapp.h"
#include "servnobjm.h"
#include "tvector.h"
#include "pidrawer.h"
#include "nomgadapter.h"

#include "skymap.h"
#include "sphericaltransformserver.h"

extern "C" {
  void skymapmodule_init();
  void skymapmodule_end();
}

class skymapmoduleExecutor : public CmdExecutor {
public:
                skymapmoduleExecutor();
  virtual       ~skymapmoduleExecutor();
  virtual int   Execute(string& keyw, vector<string>& args, string& toks);
  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);
  void MaskMap(vector<string>& tokens);
  void MapProj(vector<string>& tokens);
  void MapOper(vector<string>& tokens);
};

/* --Methode-- */
skymapmoduleExecutor::skymapmoduleExecutor()
{
PIACmd * mpiac;
NamedObjMgr omg;
mpiac = omg.GetImgApp()->CmdInterpreter();

string hgrp = "SkyMapCmd";

string kw = "map2double";
string usage = "Convert a float map to a double map";
usage += "\n Usage: map2double map";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "mapmult";
usage = "Multiply a map by a constant";
usage += "\n Usage: mapmult map 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 += "\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);

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 += "\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";
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    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";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "maskmap";
usage = "Mask a map with a mask map";
usage += "\n Usage: maskmap map msk";
usage += "\n    operation is map(i) *= msk(theta,phi)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "maproj";
usage = "Project a map into an other one";
usage += "\n Usage: maproj map mapproj [nside]";
usage += "\n   Create mapproj(nside) and project map into mapproj";
usage += "\n   (use \"maproj map mapproj\" to make a copy)";
mpiac->RegisterCommand(kw, usage, this, hgrp);

kw = "mapop";
usage = "operation on a map";
usage += "\n Usage: mapop map op1 map1 [op2 map2] [op2 map3] [...]";
usage += "\n   Do \"map op1= map1\", \"map op2=map2\", ...";
usage += "\n      where op = +,-,*,/";
mpiac->RegisterCommand(kw, usage, this, hgrp);

}

/* --Methode-- */
skymapmoduleExecutor::~skymapmoduleExecutor()
{
}

/* --Methode-- */
int skymapmoduleExecutor::Execute(string& kw, vector<string>& tokens, string& toks)
{

if(kw == "map2double") {
  if(tokens.size()<1) {
    cout<<"Usage: map2double map"<<endl;
    return(0);
  }
  Map2Double(tokens);
} else if(kw == "mapmult") {
  if(tokens.size()<2) {
    cout<<"Usage: mapmult map cste"<<endl;
    return(0);
  }
  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) {
    cout<<"Usage: map2cl map clvec [lmax] [niter]"<<endl;
    return(0);
  }
  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) {
    cout<<"Usage: crmapmask mapmsk nside lat1 lat2 [valmsk=0] [valnomsk=1]"<<endl;
    return(0);
  }
  CrMapMask(tokens);
} else if(kw == "crmaskfrmap") {
  if(tokens.size()<3) {
    cout<<"Usage: crmaskfrmap mapmsk nside map(double) v1 v2 [valmsk=0] [valnomsk=1]"<<endl;
    return(0);
  }
  CrMaskFrMap(tokens);
} else if(kw == "maskmap") {
  if(tokens.size()<2) {
    cout<<"Usage: maskmap map msk"<<endl;
    return(0);
  }
  MaskMap(tokens);
} else if(kw == "maproj") {
  if(tokens.size()<2) {
    cout<<"Usage: maproj map mapproj [nside]"<<endl;
    return(0);
  }
  MapProj(tokens);
} else if(kw == "mapop") {
  if(tokens.size()<3) {
    cout<<"Usage: mapop map op map1 [op map2] [op map3] [...]"<<endl;
    return(0);
  }
  MapOper(tokens);
}

return(0);
}

/* --Methode-- */
void skymapmoduleExecutor::Map2Double(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_4>* map = dynamic_cast<SphereHEALPix<r_4>*>(obj);
 if(map==NULL) {
   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_4>"<<endl;
   return;
 }
 int_4 nside = map->SizeIndex();
 SphereHEALPix<r_8>* map8 = new SphereHEALPix<r_8>(nside);
 cout<<"Map2Double: converting to double "<<tokens[0]<<endl;
 for(int_4 i=0;i<map8->NbPixels();i++) (*map8)(i) = (r_8) (*map)(i);
 omg.AddObj(map8,tokens[0]);
}

/* --Methode-- */
void skymapmoduleExecutor::MapMult(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 v = atof(tokens[1].c_str());
 cout<<"MapMult: Multiply "<<tokens[0]<<" by "<<v<<endl;
 for(int_4 i=0;i<map->NbPixels();i++) (*map)(i) *= v;
}

/* --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)
{
 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<<"Map2Cl: "<<tokens[0]<<" -> "<<tokens[1]
     <<" 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]);
}

/* --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]);
}

/* --Methode-- */
void skymapmoduleExecutor::CrMapMask(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;}

 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());

 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());

 cout<<"CrMaskFrMap "<<tokens[0]<<" nside"<<nside<<" with "<<tokens[2]<<endl
     <<"   for value in ["<<v1<<","<<v2<<"]"
     <<" valmsk="<<valmsk<<" unvalmsk="<<unvalmsk<<endl;

 for(int_4 i=0;i<msk->NbPixels();i++) {
   double theta,phi;
   msk->PixThetaPhi(i,theta,phi);
   int_4 ip = map->PixIndexSph(theta,phi);
   if(ip<0 || ip>=map->NbPixels()) continue;
   double v = (*map)(ip);
   if(v>v1 && v<v2) (*msk)(i)=valmsk;
     else           (*msk)(i)=unvalmsk;
 }

 omg.AddObj(msk,tokens[0]);
}

/* --Methode-- */
void skymapmoduleExecutor::MaskMap(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;
 }

 AnyDataObj* objm=omg.GetObj(tokens[1]);
 if(obj==NULL) {
   cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
   return;
 }
 SphereHEALPix<r_8>* msk = dynamic_cast<SphereHEALPix<r_8>*>(objm);
 if(msk==NULL) {
   cout<<"Error: "<<tokens[1]<<" is not a SphereHEALPix<r_8>"<<endl;
   return;
 }

 cout<<"MskMap: "<<tokens[0]<<" with mask "<<tokens[1]<<endl;

 for(int_4 i=0;i<map->NbPixels();i++) {
   double theta,phi;
   map->PixThetaPhi(i,theta,phi);
   int_4 ip = msk->PixIndexSph(theta,phi);
   if(ip<0 || ip>=msk->NbPixels()) continue;
   (*map)(i) *= (*msk)(ip);
 }

}

/* --Methode-- */
void skymapmoduleExecutor::MapProj(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 nsidemap = map->SizeIndex();

 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;}

 SphereHEALPix<r_8>* mapproj = new SphereHEALPix<r_8>(nside);

 if(nsidemap==nside) {  // tailles egales
   for(int_4 i=0;i<mapproj->NbPixels();i++)
      (*mapproj)(i) = (*map)(i);
  } else if(nsidemap<nside) {  // mapproj>map
   for(int_4 i=0;i<mapproj->NbPixels();i++) {
     double theta,phi;
     mapproj->PixThetaPhi(i,theta,phi);
     int_4 ip = map->PixIndexSph(theta,phi);
     if(ip<0 || ip>=map->NbPixels()) continue;
     (*mapproj)(i) = (*map)(ip);
   }
 } else {  // mapproj<map
   SphereHEALPix<int_4> nproj(nside);
   nproj = 0; *mapproj = 0.;
   for(int_4 i=0;i<map->NbPixels();i++) {
     double theta,phi;
     map->PixThetaPhi(i,theta,phi);
     int_4 ip = mapproj->PixIndexSph(theta,phi);
     if(ip<0 || ip>=mapproj->NbPixels()) continue;
     (*mapproj)(ip) += (*map)(i);
     nproj(ip)++;
   }
   for(int_4 i=0;i<mapproj->NbPixels();i++)
      (*mapproj)(i) /= (r_8) nproj(i);
 }
 
 omg.AddObj(mapproj,tokens[1]);
}

/* --Methode-- */
void skymapmoduleExecutor::MapOper(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;
 }

 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;
   AnyDataObj* obj1=omg.GetObj(tokens[iop+1]);
   if(obj1==NULL)
     {cout<<"Error: Pas d'objet de nom "<<tokens[iop+1]<<endl;
      continue;}
   SphereHEALPix<r_8>* map1 = dynamic_cast<SphereHEALPix<r_8>*>(obj1);
   if(map1==NULL)
     {cout<<"Error: "<<tokens[iop+1]<<" is not a SphereHEALPix<r_8>"<<endl;
      continue;}
   cout<<" "<<op<<"= "<<tokens[iop+1];
   for(int_4 i=0;i<map->NbPixels();i++) {
     double theta,phi;
     map->PixThetaPhi(i,theta,phi);
     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.;}
   }
 }
 cout<<endl;
}

////////////////////////////////////////////////////////////
static skymapmoduleExecutor * piaskmex = NULL;

void skymapmodule_init()
{
if (piaskmex) delete piaskmex;
piaskmex = new skymapmoduleExecutor;
}

void skymapmodule_end()
{
// Desactivation du module
if (piaskmex) delete piaskmex;
piaskmex = NULL;
}

////////////////////////////////////////////////////////////