Changeset 2162 in Sophya for trunk/SophyaPI


Ignore:
Timestamp:
Aug 6, 2002, 4:14:30 PM (23 years ago)
Author:
cmv
Message:

more methods cmv 6/8/02

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaPI/ProgPI/skymapmodule.cc

    r2156 r2162  
    1919#include "skymap.h"
    2020#include "sphericaltransformserver.h"
     21
     22static bool IsNSideGood(int_4 nside) {
     23 if(nside<=1) return false;
     24 while(nside>1) {if(nside%2!=0) return false; else nside/=2;}
     25 return true;
     26}
    2127
    2228extern "C" {
    … …  
    3440  void MapCover(vector<string>& tokens);
    3541  void Map2Cl(vector<string>& tokens);
     42  void Cl2Map(vector<string>& tokens);
    3643  void Map2Alm(vector<string>& tokens);
     44  void Alm2Map(vector<string>& tokens);
    3745  void CrMapMask(vector<string>& tokens);
    3846  void CrMaskFrMap(vector<string>& tokens);
    3947  void MaskMap(vector<string>& tokens);
    4048  void MapProj(vector<string>& tokens);
     49  void Map2Local(vector<string>& tokens);
    4150  void MapOper(vector<string>& tokens);
     51  void MapStat(vector<string>& tokens);
    4252};
    4353
    … …  
    6373kw = "mapcover";
    6474usage = "Print the covering of a map";
    65 usage += "\n Usage: mapcover map v1 v2 [varname]";
    66 usage += "\n    v1 v2: good pixels have v1<=val<=v2";
     75usage += "\n Usage: mapcover map v1,v2 [varname]";
     76usage += "\n    v1,v2: good pixels have v1<=val<=v2 (def=0.99,1.01)";
     77usage += "\n          (use default (0.99,1.01): \"!\")";
    6778usage += "\n    $varname: percent of sky covering";
    6879mpiac->RegisterCommand(kw, usage, this, hgrp);
    … …  
    7182usage = "SkyMap to Cl";
    7283usage += "\n Usage: map2cl map(double) clvec [lmax] [niter]";
     84usage += "\n   lmax: if <=0 or \"!\" use default (3*nside-1)";
     85usage += "\n   niter: number of iterations (<=0 or def=3)";
     86mpiac->RegisterCommand(kw, usage, this, hgrp);
     87
     88kw = "cl2map";
     89usage = "Generate SkyMap from Cl";
     90usage += "\n Usage: cl2map clvec map nside [fwhm]";
    7391mpiac->RegisterCommand(kw, usage, this, hgrp);
    7492
    … …  
    7694usage = "SkyMap to Alm";
    7795usage += "\n Usage: map2alm map(double) almmat [lmax] [niter]";
     96usage += "\n   lmax: if <=0 or \"!\" use default (3*nside-1)";
     97usage += "\n   niter: number of iterations (<=0 or def=3)";
     98mpiac->RegisterCommand(kw, usage, this, hgrp);
     99
     100kw = "alm2map";
     101usage = "Generate SkyMap from Alm";
     102usage += "\n Usage: alm2map almmat map nside";
    78103mpiac->RegisterCommand(kw, usage, this, hgrp);
    79104
    80105kw = "crmapmask";
    81 usage = "Create a map mask (nside) between latitudes lat1 lat2";
    82 usage += "\n Usage: crmapmask mapmsk nside lat1 lat2 [valmsk=0] [valnomsk=1]";
     106usage = "Create a map mask (nside) between latitudes lat1,lat2 longitudes lon1,lon2";
     107usage += "\n Usage: crmapmask mapmsk nside lat1,lat2 lon1,lon2 [valmsk,valnomsk]";
    83108usage += "\n    mask pixels between [lat1,lat2] ([-90,90] in degrees)";
    84 usage += "\n    valmsk=value to be written into masked pixels";
    85 usage += "\n    valnomsk=value to be written into unmasked pixels";
     109usage += "\n                and    [lon1,lon2] ([0,360[ in degrees)";
     110usage += "\n                (for no mask \"!\")";
     111usage += "\n    valmsk=value to be written into masked pixels (def=0)";
     112usage += "\n    valnomsk=value to be written into unmasked pixels (def=1)";
    86113mpiac->RegisterCommand(kw, usage, this, hgrp);
    87114
    88115kw = "crmaskfrmap";
    89116usage = "Create a map mask (nside) relative to an other map pixels values";
    90 usage += "\n Usage: crmaskfrmap mapmsk nside map(double) v1 v2 [valmsk=0] [valnomsk=1]";
     117usage += "\n Usage: crmaskfrmap mapmsk nside map(double) v1,v2 [valmsk,valnomsk]";
    91118usage += "\n    mask if v1<mapmsk(i)<v2";
    92 usage += "\n    valmsk=value to be written into masked pixels";
    93 usage += "\n    valnomsk=value to be written into unmasked pixels";
     119usage += "\n    valmsk=value to be written into masked pixels (def=0)";
     120usage += "\n    valnomsk=value to be written into unmasked pixels (def=1)";
    94121mpiac->RegisterCommand(kw, usage, this, hgrp);
    95122
    … …  
    107134mpiac->RegisterCommand(kw, usage, this, hgrp);
    108135
     136kw = "map2local";
     137usage = "Project a map into a local map";
     138usage += "\n Usage: map2local map localmap nx,ny angleX,angleY phi0,theta0 [x0,y0] [angle]";
     139usage += "\n   nx,ny: number of pixels in x(col),y(row) direction";
     140usage += "\n          X-axis==phi, Y-axis==theta";
     141usage += "\n   angleX,angleY: total angle aperture in x,y direction (degrees)";
     142usage += "\n   phi0,theta0: define origin in phi,theta (degrees)";
     143usage += "\n   x0,y0: map phi0,theta0 to pixel x0,y0 (\"!\" or def=middle of the localmap)";
     144usage += "\n   angle: angle (degrees) of rotation between x-axis of  map-coordinate system";
     145usage += "\n           and the tangent to parallel on the sphere (def: 0.)";
     146mpiac->RegisterCommand(kw, usage, this, hgrp);
     147
    109148kw = "mapop";
    110149usage = "operation on a map";
    … …  
    114153mpiac->RegisterCommand(kw, usage, this, hgrp);
    115154
     155kw = "mapstat";
     156usage = "Statistics from a map";
     157usage += "\n Usage: mapstat map [msk] [meanvarname] [sigmavarname]";
     158usage += "\n   msk = mask sphere used as weight";
     159usage += "\n     if msk(i)<=0 do not use that pixel";
     160usage += "\n     if msk(i)>0  use that pixel with weight msk(i)";
     161usage += "\n   msk = \"!\" means no mask sphere";
     162usage += "\n   meanvarname  = variable name to store mean";
     163usage += "\n   sigmavarname = variable name to store sigma";
     164usage += "\n   (\"!\" means do not store)";
     165mpiac->RegisterCommand(kw, usage, this, hgrp);
     166
    116167}
    117168
    … …  
    139190} else if(kw == "mapcover") {
    140191  if(tokens.size()<1) {
    141     cout<<"Usage: mapcover map v1 v2 [varname]"<<endl;
     192    cout<<"Usage: mapcover map v1,v2 [varname]"<<endl;
    142193    return(0);
    143194  }
    … …  
    149200  }
    150201  Map2Cl(tokens);
     202} else if(kw == "cl2map") {
     203  if(tokens.size()<2) {
     204    cout<<"Usage: cl2map clvec map nside [fwhm]"<<endl;
     205    return(0);
     206  }
     207  Cl2Map(tokens);
    151208} else if(kw == "map2alm") {
    152209  if(tokens.size()<2) {
    … …  
    155212  }
    156213  Map2Alm(tokens);
     214} else if(kw == "alm2map") {
     215  if(tokens.size()<3) {
     216    cout<<"Usage: alm2map almmat map nside"<<endl;
     217    return(0);
     218  }
     219  Alm2Map(tokens);
    157220} else if(kw == "crmapmask") {
    158221  if(tokens.size()<2) {
    159     cout<<"Usage: crmapmask mapmsk nside lat1 lat2 [valmsk=0] [valnomsk=1]"<<endl;
     222    cout<<"Usage: crmapmask mapmsk nside lat1,lat2 lon1,lon2 [valmsk,valnomsk]"<<endl;
    160223    return(0);
    161224  }
    … …  
    163226} else if(kw == "crmaskfrmap") {
    164227  if(tokens.size()<3) {
    165     cout<<"Usage: crmaskfrmap mapmsk nside map(double) v1 v2 [valmsk=0] [valnomsk=1]"<<endl;
     228    cout<<"Usage: crmaskfrmap mapmsk nside map(double) v1,v2 [valmsk,valnomsk]"<<endl;
    166229    return(0);
    167230  }
    … …  
    179242  }
    180243  MapProj(tokens);
     244} else if(kw == "map2local") {
     245  if(tokens.size()<5) {
     246    cout<<"Usage: map2local map localmap nx,ny angleX,angleY phi0,theta0 [x0,y0] [angle]"<<endl;
     247    return(0);
     248  }
     249  Map2Local(tokens);
    181250} else if(kw == "mapop") {
    182251  if(tokens.size()<3) {
    … …  
    185254  }
    186255  MapOper(tokens);
     256} else if(kw == "mapstat") {
     257  if(tokens.size()<1) {
     258    cout<<"Usage: Usage: mapstat map [msk] [meanvarname] [sigmavarname]"<<endl;
     259    return(0);
     260  }
     261  MapStat(tokens);
    187262}
    188263
    … …  
    246321
    247322 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());
     323 if(tokens.size()>1) if(tokens[1][0]!='!')
     324   sscanf(tokens[1].c_str(),"%lf,%lf",&v1,&v2);
    250325 cout<<"MapCover: "<<tokens[0]<<" good px=["<<v1<<","<<v2<<"]"<<endl;
    251326 int_4 ngood=0;
    … …  
    255330 cout<<"NGood="<<ngood<<" NTot="<<map->NbPixels()
    256331     <<" -> "<<100.*per<<" %";
    257  if(tokens.size()>3) {
    258    if(omg.HasVar(tokens[3])) omg.DeleteVar(tokens[3]);
     332 if(tokens.size()>2) {
     333   if(omg.HasVar(tokens[2])) omg.DeleteVar(tokens[2]);
    259334   char str[128]; sprintf(str,"%g",per);
    260    omg.SetVar(tokens[3],(string)str);
    261    cout<<" stored into $"<<tokens[3];
     335   omg.SetVar(tokens[2],(string)str);
     336   cout<<" stored into $"<<tokens[2];
    262337 }
    263338 cout<<endl;
    … …  
    270345
    271346 int_4 lmax=0, niter=3;
    272  if(tokens.size()>2) lmax = atoi(tokens[2].c_str());
     347 if(tokens.size()>2) if(tokens[2][0]!='!')
     348   lmax = atoi(tokens[2].c_str());
    273349 if(tokens.size()>3) niter = atoi(tokens[3].c_str());
    274350 if(niter<=0) niter = 3;
    … …  
    301377
    302378/* --Methode-- */
     379void skymapmoduleExecutor::Cl2Map(vector<string>& tokens)
     380{
     381 NamedObjMgr omg;
     382
     383 AnyDataObj* obj=omg.GetObj(tokens[0]);
     384 if(obj==NULL) {
     385   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
     386   return;
     387 }
     388 TVector<r_8>* cl = dynamic_cast<TVector<r_8>*>(obj);
     389 if(cl==NULL) {
     390   cout<<"Error: "<<tokens[0]<<" is not a TVector<r_8>"<<endl;
     391   return;
     392 }
     393
     394 int nside = 128;
     395 if(tokens.size()>2) nside = atoi(tokens[2].c_str());
     396 double fwhm = 0.;
     397 if(tokens.size()>3) fwhm = atof(tokens[3].c_str());
     398
     399 SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>;
     400
     401 SphericalTransformServer<r_8> ylmserver;
     402 cout<<"Cl2Map: Generate map "<<tokens[1]<<" with nside="<<nside
     403     <<" from cl "<<tokens[0]<<" with fwhm="<<fwhm<<endl;
     404 ylmserver.GenerateFromCl(*map,nside,*cl,fwhm);
     405
     406 omg.AddObj(map,tokens[1]);
     407}
     408
     409/* --Methode-- */
    303410void skymapmoduleExecutor::Map2Alm(vector<string>& tokens)
    304411{
    … …  
    306413
    307414 int_4 lmax=0, niter=3;
    308  if(tokens.size()>2) lmax = atoi(tokens[2].c_str());
     415 if(tokens.size()>2) if(tokens[2][0]!='!')
     416   lmax = atoi(tokens[2].c_str());
    309417 if(tokens.size()>3) niter = atoi(tokens[3].c_str());
    310418 if(niter<=0) niter = 3;
    … …  
    341449
    342450/* --Methode-- */
     451void skymapmoduleExecutor::Alm2Map(vector<string>& tokens)
     452{
     453 NamedObjMgr omg;
     454
     455 AnyDataObj* obj=omg.GetObj(tokens[0]);
     456 if(obj==NULL) {
     457   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
     458   return;
     459 }
     460 TMatrix< complex<r_8> >* almmat = dynamic_cast<TMatrix< complex<r_8> >*>(obj);
     461 if(almmat==NULL) {
     462   cout<<"Error: "<<tokens[0]<<" is not a TMatrix< complex<r_8> >"<<endl;
     463   return;
     464 }
     465
     466 int lmax = almmat->NRows()-1;
     467 Alm<r_8> alm(lmax);
     468 alm = (complex<r_8>)0.;
     469 for(int i=0;i<lmax+1;i++) for(int j=0;j<=i;j++) alm(i,j)=(*almmat)(i,j);
     470
     471 int nside = 128;
     472 if(tokens.size()>2) nside = atoi(tokens[2].c_str());
     473 SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>;
     474
     475 SphericalTransformServer<r_8> ylmserver;
     476 cout<<"Alm2Map: Generate map "<<tokens[1]<<" with nside="<<nside
     477     <<" from alm "<<tokens[0]<<endl;
     478 ylmserver.GenerateFromAlm(*map,nside,alm);
     479
     480 omg.AddObj(map,tokens[1]);
     481}
     482
     483/* --Methode-- */
    343484void skymapmoduleExecutor::CrMapMask(vector<string>& tokens)
     485{
     486 NamedObjMgr omg;
     487
     488 int_4 nside = atoi(tokens[1].c_str());
     489 if(!IsNSideGood(nside)) {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}
     490
     491 SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>(nside);
     492
     493 bool lat=false, lon=false;
     494 double lat1=-99., lat2=99., lon1=-999., lon2=999.;
     495 if(tokens.size()>2) if(tokens[2][0]!='!')
     496   {lat=true; sscanf(tokens[2].c_str(),"%lf,%lf",&lat1,&lat2);}
     497 if(tokens.size()>3) if(tokens[3][0]!='!')
     498   {lon=true; sscanf(tokens[3].c_str(),"%lf,%lf",&lon1,&lon2);}
     499
     500 double valmsk=0., unvalmsk=1.;
     501 if(tokens.size()>4) sscanf(tokens[4].c_str(),"%lf,%lf",&valmsk,&unvalmsk);
     502
     503 cout<<"CrMapMask "<<tokens[0]<<" nside="<<nside;
     504 if(lat) cout<<" for latitude in ["<<lat1<<","<<lat2<<"]";
     505 if(lon) cout<<" for longitude in ["<<lon1<<","<<lon2<<"]";
     506 cout<<" valmsk="<<valmsk<<" unvalmsk="<<unvalmsk<<endl;
     507
     508 //Attention: conversion en co-latitude ==> inversion: [lat2,lat1]
     509 lat1 = (90.-lat1)*M_PI/180.; lat2 = (90.-lat2)*M_PI/180.;
     510 lon1 *= M_PI/180.; lon2 *= M_PI/180.;
     511 for(int_4 i=0;i<map->NbPixels();i++) {
     512   (*map)(i)=unvalmsk;
     513   if(!lat && !lon) continue;
     514   double theta,phi; map->PixThetaPhi(i,theta,phi);
     515   if(lat && (theta<lat2 || lat1<theta)) continue;
     516   if(lon && (  phi<lon1 || lon2<phi  )) continue;
     517   (*map)(i)=valmsk;
     518 }
     519
     520 omg.AddObj(map,tokens[0]);
     521}
     522
     523/* --Methode-- */
     524void skymapmoduleExecutor::CrMaskFrMap(vector<string>& tokens)
    344525{
    345526 NamedObjMgr omg;
    … …  
    347528 int_4 nside = atoi(tokens[1].c_str());
    348529 if(nside<=1) return;
    349  int_4 n=nside; bool ok=true;
    350  while(n>1) {if(n%2!=0) {ok=false; break;} else n/=2;}
    351  if(!ok) {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}
    352 
    353  SphereHEALPix<r_8>* map = new SphereHEALPix<r_8>(nside);
    354 
    355  double lat1=-1.e30, lat2=1.e30;
    356  if(tokens.size()>2) lat1=atof(tokens[2].c_str());
    357  if(tokens.size()>3) lat2=atof(tokens[3].c_str());
     530 if(!IsNSideGood(nside)) {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}
     531
     532 AnyDataObj* obj=omg.GetObj(tokens[2]);
     533 if(obj==NULL) {
     534   cout<<"Error: Pas d'objet de nom "<<tokens[2]<<endl;
     535   return;
     536 }
     537 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
     538 if(map==NULL) {
     539   cout<<"Error: "<<tokens[2]<<" is not a SphereHEALPix<r_8>"<<endl;
     540   return;
     541 }
     542
     543 SphereHEALPix<r_8>* msk = new SphereHEALPix<r_8>(nside);
     544
     545 double v1=-1.e100, v2=1.e100;
     546 if(tokens.size()>3) if(tokens[3][0]!='!')
     547   sscanf(tokens[3].c_str(),"%lf,%lf",&v1,&v2);
    358548
    359549 double valmsk=0., unvalmsk=1.;
    360  if(tokens.size()>4) valmsk=atof(tokens[4].c_str());
    361  if(tokens.size()>5) unvalmsk=atof(tokens[5].c_str());
    362 
    363  cout<<"CrMapMask "<<tokens[0]<<" nside="<<nside
    364      <<" for latitude in ["<<lat1<<","<<lat2<<"]"
    365      <<" valmsk="<<valmsk<<" unvalmsk="<<unvalmsk<<endl;
    366 
    367  lat1 = (90.-lat1)*M_PI/180.; lat2 = (90.-lat2)*M_PI/180.;
    368  for(int_4 i=0;i<map->NbPixels();i++) {
    369    double theta,phi;
    370    map->PixThetaPhi(i,theta,phi);
    371    if(theta<lat1 && theta>lat2) (*map)(i)=valmsk;
    372      else                       (*map)(i)=unvalmsk;
    373  }
    374 
    375  omg.AddObj(map,tokens[0]);
    376 }
    377 
    378 /* --Methode-- */
    379 void skymapmoduleExecutor::CrMaskFrMap(vector<string>& tokens)
    380 {
    381  NamedObjMgr omg;
    382 
    383  int_4 nside = atoi(tokens[1].c_str());
    384  if(nside<=1) return;
    385  int_4 n=nside; bool ok=true;
    386  while(n>1) {if(n%2!=0) {ok=false; break;} else n/=2;}
    387  if(!ok) {cout<<"CrMapMask: Bad nside "<<nside<<endl; return;}
    388 
    389  AnyDataObj* obj=omg.GetObj(tokens[2]);
    390  if(obj==NULL) {
    391    cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
    392    return;
    393  }
    394  SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
    395  if(map==NULL) {
    396    cout<<"Error: "<<tokens[1]<<" is not a SphereHEALPix<r_8>"<<endl;
    397    return;
    398  }
    399 
    400  SphereHEALPix<r_8>* msk = new SphereHEALPix<r_8>(nside);
    401 
    402  double v1=-1.e50, v2=1.e50;
    403  if(tokens.size()>3) v1 = atof(tokens[3].c_str());
    404  if(tokens.size()>4) v2 = atof(tokens[4].c_str());
    405 
    406  double valmsk=0., unvalmsk=1.;
    407  if(tokens.size()>5) valmsk=atof(tokens[5].c_str());
    408  if(tokens.size()>6) unvalmsk=atof(tokens[6].c_str());
     550 if(tokens.size()>4) sscanf(tokens[4].c_str(),"%lf,%lf",&valmsk,&unvalmsk);
    409551
    410552 cout<<"CrMaskFrMap "<<tokens[0]<<" nside"<<nside<<" with "<<tokens[2]<<endl
    … …  
    413555
    414556 for(int_4 i=0;i<msk->NbPixels();i++) {
    415    double theta,phi;
    416    msk->PixThetaPhi(i,theta,phi);
     557   double theta,phi; msk->PixThetaPhi(i,theta,phi);
    417558   int_4 ip = map->PixIndexSph(theta,phi);
    418559   if(ip<0 || ip>=map->NbPixels()) continue;
    … …  
    455596
    456597 for(int_4 i=0;i<map->NbPixels();i++) {
    457    double theta,phi;
    458    map->PixThetaPhi(i,theta,phi);
     598   double theta,phi; map->PixThetaPhi(i,theta,phi);
    459599   int_4 ip = msk->PixIndexSph(theta,phi);
    460600   if(ip<0 || ip>=msk->NbPixels()) continue;
    … …  
    483623 int_4 nside = nsidemap;
    484624 if(tokens.size()>2) nside = atoi(tokens[2].c_str());
    485  if(nside<=1) return;
    486  int_4 n=nside; bool ok=true;
    487  while(n>1) {if(n%2!=0) {ok=false; break;} else n/=2;}
    488  if(!ok) {cout<<"MapProj: Bad nside "<<nside<<endl; return;}
     625 if(!IsNSideGood(nside)) {cout<<"MapProj: Bad nside "<<nside<<endl; return;}
    489626
    490627 SphereHEALPix<r_8>* mapproj = new SphereHEALPix<r_8>(nside);
    … …  
    495632  } else if(nsidemap<nside) {  // mapproj>map
    496633   for(int_4 i=0;i<mapproj->NbPixels();i++) {
    497      double theta,phi;
    498      mapproj->PixThetaPhi(i,theta,phi);
     634     double theta,phi; mapproj->PixThetaPhi(i,theta,phi);
    499635     int_4 ip = map->PixIndexSph(theta,phi);
    500636     if(ip<0 || ip>=map->NbPixels()) continue;
    … …  
    505641   nproj = 0; *mapproj = 0.;
    506642   for(int_4 i=0;i<map->NbPixels();i++) {
    507      double theta,phi;
    508      map->PixThetaPhi(i,theta,phi);
     643     double theta,phi; map->PixThetaPhi(i,theta,phi);
    509644     int_4 ip = mapproj->PixIndexSph(theta,phi);
    510645     if(ip<0 || ip>=mapproj->NbPixels()) continue;
    … …  
    517652 
    518653 omg.AddObj(mapproj,tokens[1]);
     654}
     655
     656/* --Methode-- */
     657void skymapmoduleExecutor::Map2Local(vector<string>& tokens)
     658{
     659 NamedObjMgr omg;
     660
     661 AnyDataObj* obj=omg.GetObj(tokens[0]);
     662 if(obj==NULL) {
     663   cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
     664   return;
     665 }
     666 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
     667 if(map==NULL) {
     668   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
     669   return;
     670 }
     671
     672 int_4 nx=10,ny=10;
     673 if(tokens.size()>2) sscanf(tokens[2].c_str(),"%d,%d",&nx,&ny);
     674 if(nx<=0)
     675   {cout<<"Error: nx="<<nx<<" value >0"<<endl; return;}
     676 if(ny<=0)
     677   {cout<<"Error: ny="<<ny<<" value >0"<<endl; return;}
     678
     679 double angleX=1., angleY=1.;
     680 if(tokens.size()>3) sscanf(tokens[3].c_str(),"%lf,%lf",&angleX,&angleY);
     681 if(angleX<=0. || angleX>180.)
     682   {cout<<"Error: bad angleX="<<angleX<<" value #]0,180]"<<endl; return;}
     683 if(angleY<=0. || angleY>180.)
     684   {cout<<"Error: bad angleY="<<angleY<<" value #]0,180]"<<endl; return;}
     685
     686 double phi0=0., theta0=0.;
     687 if(tokens.size()>4) sscanf(tokens[4].c_str(),"%lf,%lf",&phi0,&theta0);
     688 if(phi0<0. || phi0>=360.)
     689   {cout<<"Error: bad phi0="<<phi0<<" value #[0,360["<<endl; return;}
     690 if(theta0<-90. || theta0>90.)
     691   {cout<<"Error: bad theta0="<<theta0<<" value #[-90,90]"<<endl; return;}
     692
     693 int_4 x0=nx/2, y0=ny/2;
     694 if(tokens.size()>5) if(tokens[5][0]!='!')
     695   sscanf(tokens[5].c_str(),"%d,%d",&x0,&y0);
     696
     697 double angle=0.;
     698 if(tokens.size()>6) sscanf(tokens[6].c_str(),"%lf",&angle);
     699
     700 cout<<"Map2Local: proj "<<tokens[0]<<" to local "<<tokens[1]<<endl
     701     <<" nx="<<nx<<" ny="<<ny<<" angleX="<<angleX<<" angleY="<<angleY
     702     <<" phi0="<<phi0<<" theta0="<<theta0<<" x0="<<x0<<" y0="<<y0
     703     <<" angle="<<angle<<endl;
     704 if(angle<-90. || angle>90.)
     705   {cout<<"Error: bad angle="<<angle<<" value #[-90,90]"<<endl; return;}
     706
     707 theta0 = (90.-theta0);
     708 LocalMap<r_8>* lmap = new LocalMap<r_8>(nx,ny,angleX,angleY,theta0,phi0,x0,y0,angle);
     709 *lmap = 0.;  //lmap->print(cout);
     710
     711 int_4 nbad=0;
     712 double tmax=-9999., pmax=-9999., tmin=9999., pmin=9999.;
     713 for(int_4 i=0;i<lmap->NbPixels();i++) {
     714   double theta,phi;
     715   lmap->PixThetaPhi(i,theta,phi);
     716   int_4 ip = map->PixIndexSph(theta,phi);
     717   if(ip<0 || ip>=map->NbPixels()) {nbad++; continue;}
     718   if(theta<tmin) tmin=theta; if(theta>tmax) tmax=theta;
     719   if(phi<pmin) pmin=phi; if(phi>pmax) pmax=phi;
     720   (*lmap)(i) = (*map)(ip);
     721 }
     722 if(nbad) cout<<"Warning: "<<nbad<<" bad pixels numbers"<<endl;
     723 cout<<" phi=["<<pmin/M_PI*180.<<","<<pmax/M_PI*180.<<"]="
     724               <<(pmax-pmin)/M_PI*180.
     725     <<" theta=["<<tmin/M_PI*180.<<","<<tmax/M_PI*180.<<"]="
     726               <<(tmax-tmin)/M_PI*180.<<endl;
     727
     728 omg.AddObj(lmap,tokens[1]);
    519729}
    520730
    … …  
    549759   cout<<" "<<op<<"= "<<tokens[iop+1];
    550760   for(int_4 i=0;i<map->NbPixels();i++) {
    551      double theta,phi;
    552      map->PixThetaPhi(i,theta,phi);
     761     double theta,phi; map->PixThetaPhi(i,theta,phi);
    553762     int_4 ip = map1->PixIndexSph(theta,phi);
    554763     if(ip<0 || ip>=map1->NbPixels()) continue;
    … …  
    564773}
    565774
     775/* --Methode-- */
     776void skymapmoduleExecutor::MapStat(vector<string>& tokens)
     777{
     778 NamedObjMgr omg;
     779
     780 AnyDataObj* obj=omg.GetObj(tokens[0]);
     781 if(obj==NULL) {
     782   cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
     783   return;
     784 }
     785 SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
     786 if(map==NULL) {
     787   cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
     788   return;
     789 }
     790
     791 SphereHEALPix<r_8>* msk = NULL;
     792 if(tokens.size()>1) {
     793   if(tokens[1][0]!='!') {
     794     AnyDataObj* objm=omg.GetObj(tokens[1]);
     795     if(objm==NULL) {
     796       cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
     797       return;
     798     }
     799     msk = dynamic_cast<SphereHEALPix<r_8>*>(objm);
     800     if(msk==NULL) {
     801       cout<<"Error: "<<tokens[1]<<" is not a SphereHEALPix<r_8>"<<endl;
     802       return;
     803     }
     804   }
     805 }
     806
     807 cout<<"MapStat for map "<<tokens[0];
     808 if(msk) cout<<" using mask "<<tokens[1];
     809 cout<<endl;
     810
     811 double sum=0., sum2=0., sumw=0.;
     812 int npix = map->NbPixels(), npixgood=0;
     813 for(int_4 i=0;i<npix;i++) {
     814   double w = 1.;
     815   if(msk) {
     816     double theta,phi; map->PixThetaPhi(i,theta,phi);
     817     int_4 ip = msk->PixIndexSph(theta,phi);
     818     if(ip<0 || ip>=msk->NbPixels()) w=0.;
     819       else                          w=(*msk)(ip);
     820   }
     821   if(w<=0.) continue;
     822   npixgood++; sumw += w;
     823   sum += (*map)(i)*w;
     824   sum2 += (*map)(i)*(*map)(i)*w*w;
     825 }
     826 if(sumw<=0. || npixgood==0) {
     827   sum = sum2 = sumw=0.;
     828 } else {
     829   sum /= sumw;
     830   sum2 = sum2/sumw - sum*sum;
     831   if(sum2<=0.) sum2=0.; else sum2=sqrt(sum2);
     832 }
     833 cout<<"Number of good px="<<npixgood<<" / "<<npix
     834     <<" ("<<100.*npixgood/npix<<" %)"<<endl;
     835 cout<<" mean="<<sum<<" sigma="<<sum2<<endl;
     836
     837 if(tokens.size()>2) if(tokens[2][0]!='!') {
     838   if(omg.HasVar(tokens[2])) omg.DeleteVar(tokens[2]);
     839   char str[128]; sprintf(str,"%.8f",sum);
     840   omg.SetVar(tokens[2],(string)str);
     841   cout<<" mean stored into $"<<tokens[2];
     842 }
     843 if(tokens.size()>3) if(tokens[3][0]!='!') {
     844   if(omg.HasVar(tokens[3])) omg.DeleteVar(tokens[3]);
     845   char str[128]; sprintf(str,"%.8f",sum2);
     846   omg.SetVar(tokens[3],(string)str);
     847   cout<<" sigma stored into $"<<tokens[3];
     848 }
     849 if(tokens.size()>2) cout<<endl;
     850
     851}
     852
    566853////////////////////////////////////////////////////////////
    567854static skymapmoduleExecutor * piaskmex = NULL;
Note: See TracChangeset for help on using the changeset viewer.