1 | #include "machdefs.h"
|
---|
2 | #include <stdio.h>
|
---|
3 | #include <stdlib.h>
|
---|
4 | #include <iostream.h>
|
---|
5 | #include <math.h>
|
---|
6 |
|
---|
7 | #include <typeinfo>
|
---|
8 | #include <vector>
|
---|
9 | #include <string>
|
---|
10 |
|
---|
11 | #include "piacmd.h"
|
---|
12 | #include "nobjmgr.h"
|
---|
13 | #include "pistdimgapp.h"
|
---|
14 | #include "servnobjm.h"
|
---|
15 | #include "tvector.h"
|
---|
16 | #include "pidrawer.h"
|
---|
17 | #include "nomgadapter.h"
|
---|
18 |
|
---|
19 | #include "skymap.h"
|
---|
20 | #include "sphericaltransformserver.h"
|
---|
21 |
|
---|
22 | extern "C" {
|
---|
23 | void skymapmodule_init();
|
---|
24 | void skymapmodule_end();
|
---|
25 | }
|
---|
26 |
|
---|
27 | class skymapmoduleExecutor : public CmdExecutor {
|
---|
28 | public:
|
---|
29 | skymapmoduleExecutor();
|
---|
30 | virtual ~skymapmoduleExecutor();
|
---|
31 | virtual int Execute(string& keyw, vector<string>& args, string& toks);
|
---|
32 | void Map2Double(vector<string>& tokens);
|
---|
33 | void MapMult(vector<string>& tokens);
|
---|
34 | void MapCover(vector<string>& tokens);
|
---|
35 | void Map2Cl(vector<string>& tokens);
|
---|
36 | void Map2Alm(vector<string>& tokens);
|
---|
37 | void CrMapMask(vector<string>& tokens);
|
---|
38 | void CrMaskFrMap(vector<string>& tokens);
|
---|
39 | void MaskMap(vector<string>& tokens);
|
---|
40 | void MapProj(vector<string>& tokens);
|
---|
41 | void MapOper(vector<string>& tokens);
|
---|
42 | };
|
---|
43 |
|
---|
44 | /* --Methode-- */
|
---|
45 | skymapmoduleExecutor::skymapmoduleExecutor()
|
---|
46 | {
|
---|
47 | PIACmd * mpiac;
|
---|
48 | NamedObjMgr omg;
|
---|
49 | mpiac = omg.GetImgApp()->CmdInterpreter();
|
---|
50 |
|
---|
51 | string hgrp = "SkyMapCmd";
|
---|
52 |
|
---|
53 | string kw = "map2double";
|
---|
54 | string usage = "Convert a float map to a double map";
|
---|
55 | usage += "\n Usage: map2double map";
|
---|
56 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
57 |
|
---|
58 | kw = "mapmult";
|
---|
59 | usage = "Multiply a map by a constant";
|
---|
60 | usage += "\n Usage: mapmult map cste";
|
---|
61 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
62 |
|
---|
63 | kw = "mapcover";
|
---|
64 | usage = "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";
|
---|
67 | usage += "\n $varname: percent of sky covering";
|
---|
68 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
69 |
|
---|
70 | kw = "map2cl";
|
---|
71 | usage = "SkyMap to Cl";
|
---|
72 | usage += "\n Usage: map2cl map(double) clvec [lmax] [niter]";
|
---|
73 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
74 |
|
---|
75 | kw = "map2alm";
|
---|
76 | usage = "SkyMap to Alm";
|
---|
77 | usage += "\n Usage: map2alm map(double) almmat [lmax] [niter]";
|
---|
78 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
79 |
|
---|
80 | kw = "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]";
|
---|
83 | usage += "\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";
|
---|
86 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
87 |
|
---|
88 | kw = "crmaskfrmap";
|
---|
89 | usage = "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]";
|
---|
91 | usage += "\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";
|
---|
94 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
95 |
|
---|
96 | kw = "maskmap";
|
---|
97 | usage = "Mask a map with a mask map";
|
---|
98 | usage += "\n Usage: maskmap map msk";
|
---|
99 | usage += "\n operation is map(i) *= msk(theta,phi)";
|
---|
100 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
101 |
|
---|
102 | kw = "maproj";
|
---|
103 | usage = "Project a map into an other one";
|
---|
104 | usage += "\n Usage: maproj map mapproj [nside]";
|
---|
105 | usage += "\n Create mapproj(nside) and project map into mapproj";
|
---|
106 | usage += "\n (use \"maproj map mapproj\" to make a copy)";
|
---|
107 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
108 |
|
---|
109 | kw = "mapop";
|
---|
110 | usage = "operation on a map";
|
---|
111 | usage += "\n Usage: mapop map op1 map1 [op2 map2] [op2 map3] [...]";
|
---|
112 | usage += "\n Do \"map op1= map1\", \"map op2=map2\", ...";
|
---|
113 | usage += "\n where op = +,-,*,/";
|
---|
114 | mpiac->RegisterCommand(kw, usage, this, hgrp);
|
---|
115 |
|
---|
116 | }
|
---|
117 |
|
---|
118 | /* --Methode-- */
|
---|
119 | skymapmoduleExecutor::~skymapmoduleExecutor()
|
---|
120 | {
|
---|
121 | }
|
---|
122 |
|
---|
123 | /* --Methode-- */
|
---|
124 | int skymapmoduleExecutor::Execute(string& kw, vector<string>& tokens, string& toks)
|
---|
125 | {
|
---|
126 |
|
---|
127 | if(kw == "map2double") {
|
---|
128 | if(tokens.size()<1) {
|
---|
129 | cout<<"Usage: map2double map"<<endl;
|
---|
130 | return(0);
|
---|
131 | }
|
---|
132 | Map2Double(tokens);
|
---|
133 | } else if(kw == "mapmult") {
|
---|
134 | if(tokens.size()<2) {
|
---|
135 | cout<<"Usage: mapmult map cste"<<endl;
|
---|
136 | return(0);
|
---|
137 | }
|
---|
138 | 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);
|
---|
145 | } else if(kw == "map2cl") {
|
---|
146 | if(tokens.size()<2) {
|
---|
147 | cout<<"Usage: map2cl map clvec [lmax] [niter]"<<endl;
|
---|
148 | return(0);
|
---|
149 | }
|
---|
150 | 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);
|
---|
157 | } else if(kw == "crmapmask") {
|
---|
158 | if(tokens.size()<2) {
|
---|
159 | cout<<"Usage: crmapmask mapmsk nside lat1 lat2 [valmsk=0] [valnomsk=1]"<<endl;
|
---|
160 | return(0);
|
---|
161 | }
|
---|
162 | CrMapMask(tokens);
|
---|
163 | } else if(kw == "crmaskfrmap") {
|
---|
164 | if(tokens.size()<3) {
|
---|
165 | cout<<"Usage: crmaskfrmap mapmsk nside map(double) v1 v2 [valmsk=0] [valnomsk=1]"<<endl;
|
---|
166 | return(0);
|
---|
167 | }
|
---|
168 | CrMaskFrMap(tokens);
|
---|
169 | } else if(kw == "maskmap") {
|
---|
170 | if(tokens.size()<2) {
|
---|
171 | cout<<"Usage: maskmap map msk"<<endl;
|
---|
172 | return(0);
|
---|
173 | }
|
---|
174 | MaskMap(tokens);
|
---|
175 | } else if(kw == "maproj") {
|
---|
176 | if(tokens.size()<2) {
|
---|
177 | cout<<"Usage: maproj map mapproj [nside]"<<endl;
|
---|
178 | return(0);
|
---|
179 | }
|
---|
180 | MapProj(tokens);
|
---|
181 | } else if(kw == "mapop") {
|
---|
182 | if(tokens.size()<3) {
|
---|
183 | cout<<"Usage: mapop map op map1 [op map2] [op map3] [...]"<<endl;
|
---|
184 | return(0);
|
---|
185 | }
|
---|
186 | MapOper(tokens);
|
---|
187 | }
|
---|
188 |
|
---|
189 | return(0);
|
---|
190 | }
|
---|
191 |
|
---|
192 | /* --Methode-- */
|
---|
193 | void skymapmoduleExecutor::Map2Double(vector<string>& tokens)
|
---|
194 | {
|
---|
195 | NamedObjMgr omg;
|
---|
196 | AnyDataObj* obj=omg.GetObj(tokens[0]);
|
---|
197 | if(obj==NULL) {
|
---|
198 | cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
|
---|
199 | return;
|
---|
200 | }
|
---|
201 | SphereHEALPix<r_4>* map = dynamic_cast<SphereHEALPix<r_4>*>(obj);
|
---|
202 | if(map==NULL) {
|
---|
203 | cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_4>"<<endl;
|
---|
204 | return;
|
---|
205 | }
|
---|
206 | int_4 nside = map->SizeIndex();
|
---|
207 | SphereHEALPix<r_8>* map8 = new SphereHEALPix<r_8>(nside);
|
---|
208 | cout<<"Map2Double: converting to double "<<tokens[0]<<endl;
|
---|
209 | for(int_4 i=0;i<map8->NbPixels();i++) (*map8)(i) = (r_8) (*map)(i);
|
---|
210 | omg.AddObj(map8,tokens[0]);
|
---|
211 | }
|
---|
212 |
|
---|
213 | /* --Methode-- */
|
---|
214 | void skymapmoduleExecutor::MapMult(vector<string>& tokens)
|
---|
215 | {
|
---|
216 | NamedObjMgr omg;
|
---|
217 | AnyDataObj* obj=omg.GetObj(tokens[0]);
|
---|
218 | if(obj==NULL) {
|
---|
219 | cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
|
---|
220 | return;
|
---|
221 | }
|
---|
222 | SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
|
---|
223 | if(map==NULL) {
|
---|
224 | cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
|
---|
225 | return;
|
---|
226 | }
|
---|
227 | double v = atof(tokens[1].c_str());
|
---|
228 | cout<<"MapMult: Multiply "<<tokens[0]<<" by "<<v<<endl;
|
---|
229 | for(int_4 i=0;i<map->NbPixels();i++) (*map)(i) *= v;
|
---|
230 | }
|
---|
231 |
|
---|
232 | /* --Methode-- */
|
---|
233 | void 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-- */
|
---|
267 | void skymapmoduleExecutor::Map2Cl(vector<string>& tokens)
|
---|
268 | {
|
---|
269 | NamedObjMgr omg;
|
---|
270 |
|
---|
271 | int_4 lmax=0, niter=3;
|
---|
272 | if(tokens.size()>2) lmax = atoi(tokens[2].c_str());
|
---|
273 | if(tokens.size()>3) niter = atoi(tokens[3].c_str());
|
---|
274 | if(niter<=0) niter = 3;
|
---|
275 |
|
---|
276 | AnyDataObj* obj=omg.GetObj(tokens[0]);
|
---|
277 | if(obj==NULL) {
|
---|
278 | cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
|
---|
279 | return;
|
---|
280 | }
|
---|
281 | SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
|
---|
282 | if(map==NULL) {
|
---|
283 | cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
|
---|
284 | return;
|
---|
285 | }
|
---|
286 |
|
---|
287 | SphericalTransformServer<r_8> ylmserver;
|
---|
288 |
|
---|
289 | int_4 nside = map->SizeIndex();
|
---|
290 | if(lmax<=0) lmax = 3*nside-1;
|
---|
291 |
|
---|
292 | cout<<"Map2Cl: "<<tokens[0]<<" -> "<<tokens[1]
|
---|
293 | <<" lmax="<<lmax<<" niter="<<niter<<endl;
|
---|
294 |
|
---|
295 | SphereHEALPix<r_8> tmpmap(*map,false);
|
---|
296 | TVector<r_8> cltmp = ylmserver.DecomposeToCl(tmpmap,lmax,0.,niter);
|
---|
297 |
|
---|
298 | TVector<r_8>* cl = new TVector<r_8>(cltmp,false);
|
---|
299 | omg.AddObj(cl,tokens[1]);
|
---|
300 | }
|
---|
301 |
|
---|
302 | /* --Methode-- */
|
---|
303 | void 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]);
|
---|
340 | }
|
---|
341 |
|
---|
342 | /* --Methode-- */
|
---|
343 | void skymapmoduleExecutor::CrMapMask(vector<string>& tokens)
|
---|
344 | {
|
---|
345 | NamedObjMgr omg;
|
---|
346 |
|
---|
347 | int_4 nside = atoi(tokens[1].c_str());
|
---|
348 | 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());
|
---|
358 |
|
---|
359 | 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());
|
---|
409 |
|
---|
410 | cout<<"CrMaskFrMap "<<tokens[0]<<" nside"<<nside<<" with "<<tokens[2]<<endl
|
---|
411 | <<" for value in ["<<v1<<","<<v2<<"]"
|
---|
412 | <<" valmsk="<<valmsk<<" unvalmsk="<<unvalmsk<<endl;
|
---|
413 |
|
---|
414 | for(int_4 i=0;i<msk->NbPixels();i++) {
|
---|
415 | double theta,phi;
|
---|
416 | msk->PixThetaPhi(i,theta,phi);
|
---|
417 | int_4 ip = map->PixIndexSph(theta,phi);
|
---|
418 | if(ip<0 || ip>=map->NbPixels()) continue;
|
---|
419 | double v = (*map)(ip);
|
---|
420 | if(v>v1 && v<v2) (*msk)(i)=valmsk;
|
---|
421 | else (*msk)(i)=unvalmsk;
|
---|
422 | }
|
---|
423 |
|
---|
424 | omg.AddObj(msk,tokens[0]);
|
---|
425 | }
|
---|
426 |
|
---|
427 | /* --Methode-- */
|
---|
428 | void skymapmoduleExecutor::MaskMap(vector<string>& tokens)
|
---|
429 | {
|
---|
430 | NamedObjMgr omg;
|
---|
431 |
|
---|
432 | AnyDataObj* obj=omg.GetObj(tokens[0]);
|
---|
433 | if(obj==NULL) {
|
---|
434 | cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
|
---|
435 | return;
|
---|
436 | }
|
---|
437 | SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
|
---|
438 | if(map==NULL) {
|
---|
439 | cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
|
---|
440 | return;
|
---|
441 | }
|
---|
442 |
|
---|
443 | AnyDataObj* objm=omg.GetObj(tokens[1]);
|
---|
444 | if(obj==NULL) {
|
---|
445 | cout<<"Error: Pas d'objet de nom "<<tokens[1]<<endl;
|
---|
446 | return;
|
---|
447 | }
|
---|
448 | SphereHEALPix<r_8>* msk = dynamic_cast<SphereHEALPix<r_8>*>(objm);
|
---|
449 | if(msk==NULL) {
|
---|
450 | cout<<"Error: "<<tokens[1]<<" is not a SphereHEALPix<r_8>"<<endl;
|
---|
451 | return;
|
---|
452 | }
|
---|
453 |
|
---|
454 | cout<<"MskMap: "<<tokens[0]<<" with mask "<<tokens[1]<<endl;
|
---|
455 |
|
---|
456 | for(int_4 i=0;i<map->NbPixels();i++) {
|
---|
457 | double theta,phi;
|
---|
458 | map->PixThetaPhi(i,theta,phi);
|
---|
459 | int_4 ip = msk->PixIndexSph(theta,phi);
|
---|
460 | if(ip<0 || ip>=msk->NbPixels()) continue;
|
---|
461 | (*map)(i) *= (*msk)(ip);
|
---|
462 | }
|
---|
463 |
|
---|
464 | }
|
---|
465 |
|
---|
466 | /* --Methode-- */
|
---|
467 | void skymapmoduleExecutor::MapProj(vector<string>& tokens)
|
---|
468 | {
|
---|
469 | NamedObjMgr omg;
|
---|
470 |
|
---|
471 | AnyDataObj* obj=omg.GetObj(tokens[0]);
|
---|
472 | if(obj==NULL) {
|
---|
473 | cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
|
---|
474 | return;
|
---|
475 | }
|
---|
476 | SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
|
---|
477 | if(map==NULL) {
|
---|
478 | cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
|
---|
479 | return;
|
---|
480 | }
|
---|
481 | int_4 nsidemap = map->SizeIndex();
|
---|
482 |
|
---|
483 | int_4 nside = nsidemap;
|
---|
484 | 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;}
|
---|
489 |
|
---|
490 | SphereHEALPix<r_8>* mapproj = new SphereHEALPix<r_8>(nside);
|
---|
491 |
|
---|
492 | if(nsidemap==nside) { // tailles egales
|
---|
493 | for(int_4 i=0;i<mapproj->NbPixels();i++)
|
---|
494 | (*mapproj)(i) = (*map)(i);
|
---|
495 | } else if(nsidemap<nside) { // mapproj>map
|
---|
496 | for(int_4 i=0;i<mapproj->NbPixels();i++) {
|
---|
497 | double theta,phi;
|
---|
498 | mapproj->PixThetaPhi(i,theta,phi);
|
---|
499 | int_4 ip = map->PixIndexSph(theta,phi);
|
---|
500 | if(ip<0 || ip>=map->NbPixels()) continue;
|
---|
501 | (*mapproj)(i) = (*map)(ip);
|
---|
502 | }
|
---|
503 | } else { // mapproj<map
|
---|
504 | SphereHEALPix<int_4> nproj(nside);
|
---|
505 | nproj = 0; *mapproj = 0.;
|
---|
506 | for(int_4 i=0;i<map->NbPixels();i++) {
|
---|
507 | double theta,phi;
|
---|
508 | map->PixThetaPhi(i,theta,phi);
|
---|
509 | int_4 ip = mapproj->PixIndexSph(theta,phi);
|
---|
510 | if(ip<0 || ip>=mapproj->NbPixels()) continue;
|
---|
511 | (*mapproj)(ip) += (*map)(i);
|
---|
512 | nproj(ip)++;
|
---|
513 | }
|
---|
514 | for(int_4 i=0;i<mapproj->NbPixels();i++)
|
---|
515 | (*mapproj)(i) /= (r_8) nproj(i);
|
---|
516 | }
|
---|
517 |
|
---|
518 | omg.AddObj(mapproj,tokens[1]);
|
---|
519 | }
|
---|
520 |
|
---|
521 | /* --Methode-- */
|
---|
522 | void skymapmoduleExecutor::MapOper(vector<string>& tokens)
|
---|
523 | {
|
---|
524 | NamedObjMgr omg;
|
---|
525 |
|
---|
526 | AnyDataObj* obj=omg.GetObj(tokens[0]);
|
---|
527 | if(obj==NULL) {
|
---|
528 | cout<<"Error: Pas d'objet de nom "<<tokens[0]<<endl;
|
---|
529 | return;
|
---|
530 | }
|
---|
531 | SphereHEALPix<r_8>* map = dynamic_cast<SphereHEALPix<r_8>*>(obj);
|
---|
532 | if(map==NULL) {
|
---|
533 | cout<<"Error: "<<tokens[0]<<" is not a SphereHEALPix<r_8>"<<endl;
|
---|
534 | return;
|
---|
535 | }
|
---|
536 |
|
---|
537 | cout<<"MapOper: "<<tokens[0];
|
---|
538 | for(unsigned short iop=1;iop<tokens.size()-1;iop+=2) {
|
---|
539 | char op = tokens[iop][0];
|
---|
540 | if(op!='+' && op!='-' && op!='*' && op!='/') continue;
|
---|
541 | AnyDataObj* obj1=omg.GetObj(tokens[iop+1]);
|
---|
542 | if(obj1==NULL)
|
---|
543 | {cout<<"Error: Pas d'objet de nom "<<tokens[iop+1]<<endl;
|
---|
544 | continue;}
|
---|
545 | SphereHEALPix<r_8>* map1 = dynamic_cast<SphereHEALPix<r_8>*>(obj1);
|
---|
546 | if(map1==NULL)
|
---|
547 | {cout<<"Error: "<<tokens[iop+1]<<" is not a SphereHEALPix<r_8>"<<endl;
|
---|
548 | continue;}
|
---|
549 | cout<<" "<<op<<"= "<<tokens[iop+1];
|
---|
550 | for(int_4 i=0;i<map->NbPixels();i++) {
|
---|
551 | double theta,phi;
|
---|
552 | map->PixThetaPhi(i,theta,phi);
|
---|
553 | int_4 ip = map1->PixIndexSph(theta,phi);
|
---|
554 | if(ip<0 || ip>=map1->NbPixels()) continue;
|
---|
555 | if(op=='+') (*map)(i) += (*map1)(ip);
|
---|
556 | else if(op=='-') (*map)(i) -= (*map1)(ip);
|
---|
557 | else if(op=='*') (*map)(i) *= (*map1)(ip);
|
---|
558 | else if(op=='/')
|
---|
559 | {if((*map1)(ip)!=0.) (*map)(i) /= (*map1)(ip);
|
---|
560 | else (*map)(i) = 0.;}
|
---|
561 | }
|
---|
562 | }
|
---|
563 | cout<<endl;
|
---|
564 | }
|
---|
565 |
|
---|
566 | ////////////////////////////////////////////////////////////
|
---|
567 | static skymapmoduleExecutor * piaskmex = NULL;
|
---|
568 |
|
---|
569 | void skymapmodule_init()
|
---|
570 | {
|
---|
571 | if (piaskmex) delete piaskmex;
|
---|
572 | piaskmex = new skymapmoduleExecutor;
|
---|
573 | }
|
---|
574 |
|
---|
575 | void skymapmodule_end()
|
---|
576 | {
|
---|
577 | // Desactivation du module
|
---|
578 | if (piaskmex) delete piaskmex;
|
---|
579 | piaskmex = NULL;
|
---|
580 | }
|
---|
581 |
|
---|
582 | ////////////////////////////////////////////////////////////
|
---|