source: Sophya/trunk/SophyaProg/PMixer/skymixer.cc@ 2681

Last change on this file since 2681 was 2615, checked in by cmv, 21 years ago

using namespace sophya enleve de machdefs.h, nouveau sopnamsp.h cmv 10/09/2004
File size: 31.6 KB
Line 
1#include "sopnamsp.h"
2#include "pmixer.h"
3
4/*!
5 * \defgroup PMixer PMixer module
6 * This module contains programs which:
7 * <UL>
8 * <LI> add several sky components, taking into account their
9 * radiation spectra and convoluting them with a given filter
10 * response : skymixer
11 * <LI> create a map with point source : extractRS
12 * <LI> generate sky components, radiation spectra and spectral
13 * response (small generator of maps) : tgsky and tgrsr
14 * </UL>
15 * A detailed description may be found at:
16 */
17/*!
18 * \ingroup PMixer
19 * \file skymixer.cc
20 *\brief \b PROGRAM \b skyMixer <BR>
21 * add several sky components, taking into account their
22 *radiation spectra and convoluting them with a given filter
23 *response
24 */
25
26// -----------------------------------------------------------------
27// ------------- Function declaration ------------------------------
28int CheckCards(DataCards & dc, string & msg, FitsOutFile& fios);
29char * BuildFITSFileName(string const & fname);
30SpectralResponse * getSpectralResponse(DataCards & dc, FitsOutFile& fios);
31RadSpectra * getEmissionSpectra(DataCards & dc, int nc, FitsOutFile& fios);
32void SKM_MergeFITSKeywords(char * flnm);
33void SKM_MergeFITSKeywords2(DVList & dvl);
34void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
35void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
36
37// to add different sky components and corresponding tools
38//----------------------------------------------------------
39template <class T>
40void addComponent(SpectralResponse& sr,
41 PixelMap<T>& finalMap,
42 PixelMap<T>& mapToAdd,
43 RadSpectra& rs, double K=1.);
44//
45template <class T>
46void addComponentBeta(SphereHEALPix<T>& finalMap,
47 SphereHEALPix<T>& mapToAdd,SpectralResponse& sr,
48 SphereHEALPix<T>& betaMap, double normFreq, double K);
49//
50template <class T>
51void integratedMap(SpectralResponse& sr,
52 SphereHEALPix<T>& betaMap, double normFreq, SphereHEALPix<T>& intBetaMap);
53
54//
55template <class T>
56void addComponentBeta(SphereHEALPix<T>& finalMap,
57 SphereHEALPix<T>& mapToAdd,
58 SphereHEALPix<T>& intBetaMap, double K);
59//
60template <class T>
61void addDipole(SpectralResponse& sr, PixelMap<T>& finalMap,
62 double theta,double phi,double amp,double temp);
63//
64// -----------------------------------------------------------------
65
66// ----- Global (static) variables ------------
67static bool rdmap = false; // true -> Read map first
68static char mapPath[256]; // Path for input maps
69static int hp_nside = 32; // HealPix NSide
70static int nskycomp = 0; // Number of sky components
71static int debuglev = 0; // Debug Level
72static int printlev = 0; // Print Level
73static POutPersist * so = NULL; // Debug PPFOut file
74static DVList * dvl_fitskw = NULL; // Global DVList for all FITS Keywords
75bool checkpdmtype=false;
76bool dvliston= false;
77// --------- SkyMixer Version --------------
78static double skm_version = 1.5;
79
80// -------------------------------------------------------------------------
81// main program
82// -------------------------------------------------------------------------
83int main(int narg, char * arg[])
84{
85 if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
86 cout << " Usage: skymixer parameterFile outputfitsname [outppfname]" << endl;
87 exit(0);
88 }
89
90 InitTim();
91 FitsOutFile fios(arg[2]);
92
93 cout << " skymixer Version= " << skm_version << " SophyaVersion= "
94 << SophyaVersion() << endl;
95
96 string msg;
97 int rc = 0;
98 RadSpectra * es = NULL;
99 SpectralResponse * sr = NULL;
100 double moy, sig;
101
102 DataCards dc;
103 so = NULL;
104 // DVList for merging all FITS keywords and datacards
105 dvl_fitskw = new DVList;
106
107 try {
108 string dcard = arg[1];
109 if(printlev > 1) cout << " Decoding parameters from file " << dcard << endl;
110 dc.ReadFile(dcard);
111
112 rc = CheckCards(dc, msg,fios);
113 if (rc) {
114 cerr << " Error condition -> Rc= " << rc << endl;
115 cerr << " Msg= " << msg << endl;
116 return(rc);
117 }
118 }
119 catch (PException exc) {
120 msg = exc.Msg();
121 cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
122 return(90);
123 }
124
125
126 cout << " skymix/Info : NComp = " << nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
127 cout << " ... MapPath = " << (string)mapPath << " DebugLev= " << debuglev
128 << " PrintLev= " << printlev << endl;
129
130 // We create an output persist file for writing debug objects
131 if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
132 int countHead = 2;
133 SphereHEALPix<float> outgs(hp_nside);
134 try{
135 if (rdmap) { // Reading map from FITS file
136 char ifnm[256];
137 strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
138 ifnm[255] = '\0';
139 cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
140 {
141 FitsInFile fitsin(ifnm);
142 fitsin >> outgs;
143 // Getting FITS keywords in primary header
144 // SKM_MergeFITSKeywords(ifnm);
145 // Getting FITS keywords in data object header
146 if(checkpdmtype)
147 {
148 int status = 1;
149 string whatsup=fitsin.getStringKeyword(2,"PDMTYPE",status);
150 if(!(status==0&&whatsup=="COMPMAP"))
151 {
152 cerr << " !!!! skymixer - file" << ifnm << " is not a COMPMAP " << endl;
153 throw ParmError("Error");
154 }
155 else cout << "read the comment COMPMAP from " << ifnm << endl;
156 }
157 if(dvliston) SKM_MergeFITSKeywords2(outgs.Info());
158 else
159 {
160 fios.appendInputHeader(fitsin, countHead);
161 }
162 }
163 if(printlev>0)
164 cout << " Output HealPIx Map read - NbPixels= " <<
165 outgs.NbPixels() << endl;
166 if (printlev > 0) {
167 MeanSig(outgs.DataBlock(), moy, sig );
168 cout << " MeanSig for outpout map - Mean= " <<
169 moy << " Sigma= " << sig << endl;
170 }
171 }
172 else {
173 if(printlev>0)
174 cout << " Output HealPix Map created - NbPixels= " <<
175 outgs.NbPixels() << endl;
176 outgs.SetPixels(0.);
177 }
178
179 // Decoding detection pass-band filter
180 sr = getSpectralResponse(dc,fios);
181 PrtTim(" After FilterCreation ");
182
183 char * flnm, buff[90];
184 string key;
185
186 double K = 1.;
187 double freqOfMap = 1.;
188 // Loop over sky component
189 int sk;
190 for(sk = 0; sk<nskycomp; sk++) {
191 cout << "------------------------------------" << endl;
192 cout << " Processing sky component No " << sk+1 << endl;
193
194
195 sprintf(buff, "%d", sk+1);
196 key = (string)"DIPOLE" + buff;
197 // check for an eventual dipole
198 if(dc.HasKey(key))
199 {
200 if (es) { delete es; es = NULL; }
201 double temp = -10.;
202 double theta= dc.DParam(key,1,1.);
203 double phi = dc.DParam(key,2,1.);
204 double amp = dc.DParam(key,3,1.);
205 if(dc.NbParam(key)>3)
206 {
207 temp = dc.DParam(key,4,1.);
208 }
209 cout << " creating dipole " << temp << " " << amp << " " << phi << " " << theta << " " << endl;
210 addDipole(*sr, outgs,theta,phi,amp,temp);
211 }
212 else
213 {
214 sprintf(buff, "%d", sk+1);
215 key = (string)"MAPFITSFILE" + buff;
216 flnm = BuildFITSFileName(dc.SParam(key, 0));
217
218 K = dc.DParam(key, 1, 1.);
219
220 cout << " Reading Input FITS map " << (string)flnm << endl;
221 SphereHEALPix<float> ings(hp_nside);
222 {
223 FitsInFile fiosIn(flnm);
224 fiosIn >> ings;
225 // Getting FITS keywords in primary header
226 // SKM_MergeFITSKeywords(flnm);
227 // Getting FITS keywords in data object header
228 //if(fiosIn.hasKeyword(string("PDMTYPE"), 1)) cout << "read the comment COMPMAP from " << flnm << endl;
229 if(checkpdmtype)
230 {
231 int status = 1;
232 string ch = "PDMTYPE";
233 string whatsup=fiosIn.getStringKeyword(2,ch,status);
234 if(!(status==0&&whatsup=="COMPMAP"))
235 {
236 cerr << " !!!! skymixer - file" << flnm << " is not a COMPMAP " << endl;
237 throw ParmError("Error");
238 }
239 else cout << "read the comment COMPMAP from " << flnm << endl;
240 }
241 if(dvliston) SKM_MergeFITSKeywords2(ings.Info());
242 else
243 {
244 fios.appendInputHeader(fiosIn, countHead);
245 }
246 }
247 if (debuglev > 4) { // Writing the input map to the outppf
248 FIO_SphereHEALPix<float> fiog(ings);
249 fiog.Write(*so, key);
250 }
251 if (printlev > 2) {
252 MeanSig(ings.DataBlock(), moy, sig );
253 cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
254 }
255 bool mapDependentOfFreq = false;
256 key = (string)"BETAFITSFILE"+ buff;
257 if(dc.HasKey(key))
258 {
259 mapDependentOfFreq = true;
260 }
261
262 // getting Emission spectra
263 if(!mapDependentOfFreq)
264 {
265 if (es) { delete es; es = NULL; }
266 es = getEmissionSpectra(dc, sk,fios);
267 addComponent(*sr, outgs, ings, *es, K);
268 }
269 else
270 {
271 key = (string)"BETAFITSFILE"+ buff;
272 //SphereHEALPix<float> betaMap;
273 flnm = BuildFITSFileName(dc.SParam(key, 0));
274 double normFreq = dc.DParam(key, 1, 1.);
275 if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
276 int nSideForInt = dc.DParam(key, 2, 1.);
277 if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
278 cout << "....BetaFits... Reading Beta FITS map " << (string)flnm << endl;
279 SphereHEALPix<float> betaMap(hp_nside);
280 {
281 FitsInFile fiosBM(flnm);
282 fiosBM >> betaMap;
283 // Getting FITS keywords in primary header
284 // SKM_MergeFITSKeywords(flnm);
285 // Getting FITS keywords in data object header
286 //if(fiosBM.hasKeyword(string("COMPMAP"), 1)) cout << "read the comment COMPMAP from " << flnm << endl;
287 if(checkpdmtype)
288 {
289 int status = 1;
290 string whatsup=fiosBM.getStringKeyword(2,"PDMTYPE",status);
291 if(!(status==0&&whatsup=="INDEXMAP"))
292 {
293 cerr << " !!!! skymixer - file" << flnm << " is not an INDEXMAP " << endl;
294 throw ParmError("Error");
295 }
296 else cout << "read the comment INDEXMAP from " << flnm << endl;
297 }
298 if(dvliston) SKM_MergeFITSKeywords2(betaMap.Info());
299 else
300 {
301 fios.appendInputHeader(fiosBM, countHead);
302 }
303 }
304 if (printlev > 2) {
305 MeanSig(betaMap.DataBlock(), moy, sig );
306 cout << " MeanSig for Beta map - Mean= " << moy << " Sigma= " << sig << endl;
307 }
308 if (debuglev > 4) { // Writing the input map to the outppf
309 FIO_SphereHEALPix<float> fiogs(betaMap);
310 fiogs.Write(*so, key);
311 }
312
313 if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
314 bool bydefault = true;
315 if(!bydefault)
316 addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
317 else
318 {
319 // integrate the betamap over the SpectralResponse
320 SphereHEALPix<float> intBetaMap(nSideForInt);
321 integratedMap(*sr, betaMap, normFreq, intBetaMap);
322 if (debuglev > 4) { // Writing the input map to the outppf
323 FIO_SphereHEALPix<float> fiogs2(intBetaMap);
324 fiogs2.Write(*so, "INTBETAMAP");
325 }
326
327 betaMap.Resize(8);
328 if (printlev > 4)
329 {
330 MeanSig(intBetaMap.DataBlock(), moy, sig );
331 cout << "....BetaFits... MeanSig for intBetaMap - Mean= "
332 << moy << " Sigma= " << sig << endl;
333 }
334 // add the integrated beta map
335 addComponentBeta(outgs,ings,intBetaMap, K);
336 }
337 }
338
339 MeanSig(outgs.DataBlock(), moy, sig );
340 cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
341 cout << "-------------------------------------------------" << endl;
342
343 sprintf(buff, "End of Processing Component %d ", sk+1);
344 PrtTim(buff);
345 }
346 }
347 }
348 catch(PException exc)
349 {
350 cout << "catched PException" << endl;
351 msg = exc.Msg();
352 cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
353 rc = 50;
354 return(50);
355 }
356
357 try {
358 // Saving the output map in FITS format
359 cout << "Output Map (SphereHEALPix<float>) written to FITS file "
360 << (string)(arg[2]) << endl;
361 {
362 // FitsOutFile fios(arg[2]);
363 fios.firstImageOnPrimaryHeader(false); // Use secondary header
364 DVList& dvl = (*dvl_fitskw);
365 dvl["PDMTYPE"] = "COMPMAP";
366 dvl.SetComment("PDMTYPE", "Planck Data Model Type");
367 dvl["SOPHYVER"] = SophyaVersion();
368 dvl.SetComment("SOPHYVER", "Sophya Version number");
369 dvl["SKYMVER"] = skm_version;
370 dvl.SetComment("SKYMVER", "skymixer Version number");
371 // fios.DVListIntoPrimaryHeader(dvl);
372 dvl["CREATOR"] = "SkyMixer";
373 outgs.Info() = dvl;
374 fios << outgs ;
375 }
376 PrtTim("End of WriteFITS ");
377 // Saving the output map in PPF format
378 if (narg > 3) {
379 POutPersist s(arg[3]);
380 FIO_SphereHEALPix<float> fiog(&outgs) ;
381 fiog.Write(s);
382 cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "
383 << (string)(arg[3]) << endl;
384 PrtTim("End of WritePPF ");
385 }
386 }
387 catch(PException exc)
388 {
389 cout << "catched PException (2)" << endl;
390 msg = exc.Msg();
391 cerr << " !!!! skymixer(2) - Catched exception - Msg= " << exc.Msg() << endl;
392 rc = 55;
393 return(55);
394 }
395
396 if (so) delete so; // Closing the debug ppf file
397 return(rc);
398}
399
400/* Nouvelle-Fonction */
401int CheckCards(DataCards & dc, string & msg, FitsOutFile& fios)
402// Function to check datacards
403{
404rdmap = false;
405mapPath[0] = '\0';
406hp_nside = 32;
407nskycomp = 0;
408debuglev = 0;
409printlev = 0;
410
411int rc = 0;
412string key, key2,key3;
413
414 DVList & dvl = *dvl_fitskw;
415 // Cheking datacards
416 if (dc.NbParam("SKYMIX") < 2) {
417 rc = 71;
418 msg = "Invalid parameters - Check @SKYMIX card ";
419 return(rc);
420 }
421 if(dvliston)
422 dvl.SetS("SKYMIX", dc.GetParams("SKYMIX"));
423 else
424 {
425 fios.insertCommentLineOnHeader("--------------------");
426 fios.insertCommentLineOnHeader("SkyMixer DataCard");
427 fios.insertCommentLineOnHeader("--------------------");
428
429 string key = "SKYMIX"+ (string)dc.GetParams("SKYMIX");
430 fios.insertCommentLineOnHeader(key);
431 }
432
433 key = "CHECKPDMTYPE";
434 if (dc.HasKey(key)) {
435 checkpdmtype=true;
436 }
437
438
439 key = "READMAP";
440 if (dc.HasKey(key)) {
441 if (dc.NbParam(key) < 1) {
442 rc = 72;
443 msg = "Invalid parameters - Check @READMAP card ";
444 return(rc);
445 }
446 else rdmap = true;
447
448 if(dvliston) dvl.SetS(key, dc.GetParams(key));
449 else
450 {
451 string keynew = key + (string)dc.GetParams(key);
452 fios.insertCommentLineOnHeader(keynew);
453 }
454 }
455
456// Checking detection filter specification
457 key = "GAUSSFILTER";
458 key2 = "FILTERFITSFILE";
459 key3 = "DIPOLE";
460 if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3) && (dc.NbParam(key3) < 3)) {
461 msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE or @DIPOLE";
462 rc = 73; return(rc);
463 }
464
465 if (dc.HasKey(key))
466 {
467 if(dvliston) dvl.SetS("GAUSFILT", dc.GetParams(key));
468 else
469 {
470 string keynew = key + (string)dc.GetParams(key);
471 fios.insertCommentLineOnHeader(keynew);
472 }
473 }
474 if (dc.HasKey(key2))
475 {
476 if(dvliston) dvl.SetS("FILTFILE", dc.GetParams(key2));
477 else
478 {
479 string keynew = key2 + (string)dc.GetParams(key2);
480 fios.insertCommentLineOnHeader(keynew);
481 }
482 }
483 if (dc.HasKey(key3))
484 {
485 if(dvliston) dvl.SetS("DIPOLE", dc.GetParams(key3));
486 else
487 {
488 string keynew = key3 + (string)dc.GetParams(key3);
489 fios.insertCommentLineOnHeader(keynew);
490 }
491 }
492 // Decoding number of component and pixelisation parameter
493 int mg = 32;
494 int ncomp = 0;
495 ncomp = dc.IParam("SKYMIX", 0, 0);
496 mg = dc.IParam("SKYMIX", 1, 32);
497 if (ncomp < 1) {
498 msg = "Invalid parameters - Check datacards @SKYMIX ";
499 rc = 74;
500 return(rc);
501 }
502
503 // Checking detection filter specification
504 // Checking input FITS file specifications
505 int kc;
506 char buff[256];
507 bool pb = false;
508 string key4;
509 string key5;
510 string key6;
511 for(kc=0; kc<ncomp; kc++) {
512 sprintf(buff, "MAPFITSFILE%d", kc+1);
513 key = buff;
514 sprintf(buff, "DIPOLE%d", kc+1);
515 key3 = buff;
516 if (dc.NbParam(key) < 1 && (dc.NbParam(key3)<1))
517 {
518 msg = "Invalid card or Missing card" + key + " or " + key3 ;
519 pb = true; break;
520 }
521
522 sprintf(buff, "MAPFIL%d", kc+1);
523 if (dc.HasKey(key))
524 {
525 if(dvliston) dvl.SetS(buff, dc.GetParams(key));
526 else
527 {
528 string keynew = buff + (string)dc.GetParams(key);
529 fios.insertCommentLineOnHeader(keynew);
530 }
531 }
532 sprintf(buff, "DIPOLE%d", kc+1);
533 if (dc.HasKey(key3))
534 {
535 if(dvliston) dvl.SetS(buff, dc.GetParams(key3));
536 else
537 {
538 string keynew = buff + (string)dc.GetParams(key3);
539 fios.insertCommentLineOnHeader(keynew);
540 }
541 }
542 sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
543 key = buff;
544 sprintf(buff, "BLACKBODY%d", kc+1);
545 key2 = buff;
546 sprintf(buff, "POWERLAWSPECTRA%d", kc+1);
547 key3 = buff;
548 sprintf(buff, "BETAFITSFILE%d", kc+1);
549 key4 = buff;
550 sprintf(buff, "DIPOLE%d", kc+1);
551 key5 = buff;
552 sprintf(buff, "DERIVBB%d", kc+1);
553 key6 = buff;
554 if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
555 && (dc.NbParam(key6)<1) && (dc.NbParam(key5)<1)) {
556 msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
557 pb = true; break;
558 }
559
560 sprintf(buff, "SPECTF%d", kc+1);
561 if (dc.HasKey(key)) {
562 if(dvliston) dvl.SetS(buff, dc.GetParams(key));
563 else {
564 string keynew = buff + (string)dc.GetParams(key);
565 fios.insertCommentLineOnHeader(keynew);}
566 }
567 sprintf(buff, "BLACKB%d", kc+1);
568 if (dc.HasKey(key2)) {
569 if(dvliston) dvl.SetS(buff, dc.GetParams(key2));
570 else {
571 string keynew = buff + (string)dc.GetParams(key2);
572 fios.insertCommentLineOnHeader(keynew);}
573 }
574
575 sprintf(buff, "PLAWSP%d", kc+1);
576 if (dc.HasKey(key3)) {
577 if(dvliston) dvl.SetS(buff, dc.GetParams(key3));
578 else {
579 string keynew = buff + (string)dc.GetParams(key3);
580 fios.insertCommentLineOnHeader(keynew);}
581 }
582
583 sprintf(buff, "BETAFI%d", kc+1);
584 if (dc.HasKey(key4)) {
585 if(dvliston) dvl.SetS(buff, dc.GetParams(key4));
586 else {
587 string keynew = buff + (string)dc.GetParams(key4);
588 fios.insertCommentLineOnHeader(keynew);}
589 }
590
591 sprintf(buff, "DIPOLE%d", kc+1);
592 if (dc.HasKey(key5)) {
593 if(dvliston) dvl.SetS(buff, dc.GetParams(key5));
594 else {
595 string keynew = buff + (string)dc.GetParams(key5);
596 fios.insertCommentLineOnHeader(keynew);}
597 }
598
599 sprintf(buff, "DERIBB%d", kc+1);
600 if (dc.HasKey(key6)) {
601 if(dvliston) dvl.SetS(buff, dc.GetParams(key6));
602 else {
603 string keynew = buff + (string)dc.GetParams(key6);
604 fios.insertCommentLineOnHeader(keynew);}
605 }
606
607 }
608
609 if (pb) {
610 rc = 75;
611 return(75);
612 }
613
614
615// Initialiazing parameters
616 rc = 0;
617 msg = "OK";
618 nskycomp = ncomp;
619 hp_nside = mg;
620
621
622// Checking for PATH definition card
623 key = "MAPPATH";
624 if (dc.NbParam(key) < 3) strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
625 mapPath[255] = '\0';
626 key = "DEBUGLEVEL";
627 debuglev = dc.IParam(key, 0, 0);
628 key = "PRINTLEVEL";
629 printlev = dc.IParam(key, 0, 0);
630 return(rc);
631}
632
633static char buff_flnm[1024]; // Mal protege !
634/* Nouvelle-Fonction */
635char* BuildFITSFileName(string const & fname)
636{
637if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
638else sprintf(buff_flnm, "%s", fname.c_str());
639return(buff_flnm);
640}
641
642/* Nouvelle-Fonction */
643SpectralResponse * getSpectralResponse(DataCards & dc, FitsOutFile& fios)
644{
645 SpectralResponse * filt = NULL;
646
647 string key = "FILTERFITSFILE";
648 string key2 = "GAUSSFILTER";
649 string ppfname = "filter";
650
651 if (dc.HasKey(key) ) { // Reading FITS filter file
652 char ifnm[256];
653 strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
654 ifnm[255] = '\0';
655 Matrix mtx;
656 FitsInFile fiis(ifnm);
657 fiis.firstImageOnPrimaryHeader(false); // Use secondary header HDU=2
658 fiis >> mtx ;
659 // Getting FITS keywords in primary header
660 // SKM_MergeFITSKeywords(ifnm);
661 // Getting FITS keywords in data object header
662 if(checkpdmtype)
663 {
664 int status = 1;
665 string whatsup=fiis.getStringKeyword(2,"PDMTYPE",status);
666 if(!(status==0&&whatsup=="COMPMAP"))
667 {
668 cerr << " !!!! skymixer - file" << ifnm << " is not a COMPMAP " << endl;
669 throw ParmError("Error");
670 }
671 else cout << "read the comment COMPMAP from " << ifnm << endl;
672 }
673 if(dvliston) SKM_MergeFITSKeywords2(mtx.Info());
674 else
675 {
676 fios.appendInputHeader(fiis, 2);
677 }
678 double numin = dc.DParam(key, 2, 1.);
679 double numax = dc.DParam(key, 3, 9999.);
680 Vector nu(mtx.NCols());
681 Vector fnu(mtx.NCols());
682 for(int k=0; k<mtx.NCols(); k++) {
683 nu(k) = mtx(0, k);
684 fnu(k) = mtx(1, k);
685 }
686 filt = new SpecRespVec(nu, fnu, numin, numax);
687 ppfname = key;
688 }
689 else if (dc.HasKey(key2) ) { // creating GaussianFilter
690 double nu0 = dc.DParam(key2, 0, 10.);
691 double s = dc.DParam(key2, 1, 1.);
692 double a = dc.DParam(key2, 2, 1.);
693 double numin = dc.DParam(key2, 3, 0.1);
694 double numax = dc.DParam(key2, 4, 9999);
695 filt = new GaussianFilter(nu0, s, a, numin, numax);
696 ppfname = key2;
697 }
698 if (filt == NULL) throw ParmError("datacard error ! No detection filter");
699 if(printlev>0)
700 {
701 cout << endl;
702 cout << " Filter decoded - Created " << endl;
703 cout << *filt << endl;
704 }
705 // for debug
706 if (debuglev > 1) SpectralResponse2Nt(*filt, *so, ppfname);
707 return(filt);
708}
709
710/* Nouvelle-Fonction */
711RadSpectra * getEmissionSpectra(DataCards & dc, int nc, FitsOutFile& fios)
712{
713 char numb[16];
714 sprintf(numb, "%d", nc+1);
715
716 string key = (string)"SPECTRAFITSFILE" + numb;
717 string key2 = (string)"BLACKBODY" + numb;
718 string key5 = (string)"DERIVBB" + numb;
719 string key3 = (string)"POWERLAWSPECTRA" + numb;
720 string ppfname = "espectra";
721
722 RadSpectra * rs = NULL;
723 if (dc.HasKey(key) ) { // Reading emission spectra from file
724 char * ifnm = BuildFITSFileName(dc.SParam(key, 0));
725 cout << " Reading Input FITS spectra file " << (string)ifnm << endl;
726 Matrix mtx;
727 FitsInFile fiis(ifnm);
728 fiis.firstImageOnPrimaryHeader(false); // Use secondary header HDU=2
729 fiis >> mtx ;
730 // Getting FITS keywords in primary header
731 // SKM_MergeFITSKeywords(ifnm);
732 // Getting FITS keywords in data object header
733 //if(fiis.hasKeyword(string("PDMTYPE"), 1)) cout << "read the comment COMPMAP from " << ifnm << endl;
734 if(checkpdmtype)
735 {
736 int status = 1;
737 string whatsup=fiis.getStringKeyword(2,"PDMTYPE",status);
738 if(!(status==0&&whatsup=="FGRSPEC"))
739 {
740 cerr << " !!!! skymixer - file" << ifnm << " is not a FGRSPEC" << endl;
741 throw ParmError("Error");
742 }
743 else cout << "read the comment FGRSPEC from " << ifnm << endl;
744 }
745 if(dvliston) SKM_MergeFITSKeywords2(mtx.Info());
746 else
747 {
748 fios.appendInputHeader(fiis, 2);
749 }
750 double numin = dc.DParam(key, 2, 1.);
751 double numax = dc.DParam(key, 3, 9999.);
752 Vector nu(mtx.NCols());
753 Vector tnu(mtx.NCols());
754 for(int k=0; k<mtx.NCols(); k++) {
755 nu(k) = mtx(0, k);
756 tnu(k) = mtx(1, k);
757 }
758 rs = new RadSpectraVec(nu, tnu, numin, numax);
759 ppfname = key;
760 }
761 else if (dc.HasKey(key2) ) { // Creating BlackBody emission spectra
762 rs = new BlackBody(dc.DParam(key2, 0, 2.726));
763 ppfname = key2;
764 }
765 else if (dc.HasKey(key5) ) { // Creating Dipole
766 rs = new DerivBlackBody(dc.DParam(key5, 0, 3.E-3));
767 ppfname = key5;
768 }
769 else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
770 double a = dc.DParam(key3, 0, 1.);
771 double nu0 = dc.DParam(key3, 1, 100.);
772 double dnu = dc.DParam(key3, 2, 10.);
773 double b = dc.DParam(key3, 3, 0.);
774 double numin = dc.DParam(key3, 4, 0.1);
775 double numax = dc.DParam(key3, 5, 9999);
776 rs = new PowerLawSpectra(a, b, nu0, dnu, numin, numax);
777 ppfname = key3;
778 }
779 if (rs == NULL) throw ParmError("datacard error ! missing Emission spectra");
780 cout << " Emission spectra decoded - Created (" << ppfname << ")" << endl;
781 cout << *rs << endl;
782// for debug
783 if (debuglev > 2) RadSpec2Nt(*rs, *so, ppfname);
784 return(rs);
785}
786
787
788
789
790template <class T>
791void addDipole(SpectralResponse& sr, PixelMap<T>& finalMap,
792 double theta,double phi,double amp,double temp)
793{
794 DerivBlackBody dbb;
795 if(temp>0) dbb.setTemperature(temp);
796 double coeff = dbb.filteredIntegratedFlux(sr) * amp;
797 UnitVector vd(theta,phi);
798 UnitVector vc(theta,phi);
799
800 for(int i=0; i<finalMap.NbPixels(); i++)
801 {
802 double thetar,phir;
803 finalMap.PixThetaPhi(i,thetar,phir);
804 vc.SetThetaPhi(thetar, phir);
805 finalMap(i) += vd.Psc(vc)*coeff;
806 }
807 if (debuglev > 4) { // Writing the input map to the outppf
808 SphereHEALPix<float> ings(sqrt((double)finalMap.NbPixels()/12));
809 for(int i=0; i<finalMap.NbPixels(); i++)
810 {
811 double thetar,phir;
812 finalMap.PixThetaPhi(i,thetar,phir);
813 vc.SetThetaPhi(thetar, phir);
814 ings(i) = vd.Psc(vc);
815 }
816 FIO_SphereHEALPix<float> fiog(ings);
817 fiog.Write(*so, "dipole");
818 cout << "Debug the dipole map....saved in debug file !" << endl;
819 }
820}
821/* Nouvelle-Fonction */
822template <class T>
823void addComponent(SpectralResponse& sr, PixelMap<T>& finalMap,
824 PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
825{
826 // finalMap = finalMap + coeff* mapToAdd
827 // coeff = convolution of sr and rs
828 // compute the coefficient corresponding to mapToAdd
829 if (finalMap.NbPixels() != mapToAdd.NbPixels())
830 throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
831 double coeff = rs.filteredIntegratedFlux(sr) * K;
832 if (printlev > 1)
833 cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl;
834 for(int i=0; i<finalMap.NbPixels(); i++)
835 {
836 finalMap(i) += coeff * mapToAdd(i);
837 }
838}
839/* Nouvelle-Fonction */
840template <class T>
841void addComponentBeta(SphereHEALPix<T>& finalMap,
842 SphereHEALPix<T>& mapToAdd,SpectralResponse& sr,
843 SphereHEALPix<T>& betaMap, double normFreq, double K)
844{
845 // finalMap = finalMap + coeff* mapToAdd
846 // coeff = convolution of sr and rs
847 // compute the coefficient corresponding to mapToAdd
848 // betaMap is the map of (beta(theta,phi))
849
850 int nbpix = finalMap.NbPixels();
851 if (nbpix != mapToAdd.NbPixels())
852 throw SzMismatchError("addComponentBeta()/Error: Unequal number of Input/Output map pixels");
853 if (printlev > 1)
854 {
855 cout << "addComponentBeta - Coeff= " << K << endl;
856 cout << "nb pixels: " << finalMap.NbPixels() << endl;
857 }
858 SphereHEALPix<T> bigBetaMap(sqrt((double)nbpix/12));
859 if(nbpix != betaMap.NbPixels())
860 {
861 Sph2Sph(betaMap,bigBetaMap);
862 }
863 for(int i=0; i<finalMap.NbPixels(); i++)
864 {
865 // coeff = integration of (nu/normFreq)^(-beta(theta,phi)) in each pixels
866 RadSpectra* rs = new PowerLawSpectra
867 (1.,-bigBetaMap(i), 0., normFreq, 0.01, 800.);
868 double coeff = rs->filteredIntegratedFlux(sr);
869 finalMap(i) += coeff*K*mapToAdd(i);
870 }
871}
872
873template <class T>
874void integratedMap(SpectralResponse& sr,
875 SphereHEALPix<T>& betaMap,
876 double normFreq,
877 SphereHEALPix<T>& intBetaMap)
878{
879 PowerLawSpectra rs(1.,-2., 0., normFreq);
880
881 if(betaMap.NbPixels()!=intBetaMap.NbPixels())
882 {
883 Sph2Sph(betaMap,intBetaMap);
884 for(int i=0; i<intBetaMap.NbPixels(); i++)
885 {
886 rs.setExp(-intBetaMap(i));
887 double coeff = rs.filteredIntegratedFlux(sr);
888 intBetaMap(i) = coeff;
889 }
890 }
891 else
892 {
893 for(int i=0; i<intBetaMap.NbPixels(); i++)
894 {
895 rs.setExp(-betaMap(i));
896 double coeff = rs.filteredIntegratedFlux(sr);
897 intBetaMap(i) = coeff;
898 }
899 }
900}
901
902template <class T>
903void addComponentBeta(SphereHEALPix<T>& finalMap,
904 SphereHEALPix<T>& mapToAdd,SphereHEALPix<T>& intBetaMap, double K)
905{
906 // finalMap = finalMap + coeff* mapToAdd
907 // coeff = convolution of sr and rs
908 // compute the coefficient corresponding to mapToAdd
909 // integBetaMap is the map of the integration (nu/normFreq)^(-beta(theta,phi)) over
910 // the spectralResponse
911 // different from addComponentBeta(PixelMap<T>& finalMap,
912 // PixelMap<T>& mapToAdd,SpectralResponse& sr, PixelMap<T>& betaMap, double normFreq, double K)
913 // since it permits to use a intBetaMap with a different number of pixels than
914 // the other maps
915
916 int nbpix = finalMap.NbPixels();
917 if (nbpix != mapToAdd.NbPixels())
918 throw SzMismatchError("addComponentBeta(PixelMap<T>&,PixelMap<T>&,PixelMap<T>&,double)/Error: Unequal number of Input/Output map pixels");
919 double coeff = K;
920
921 if(nbpix != intBetaMap.NbPixels())
922 {
923 for(int i=0; i<finalMap.NbPixels();i++)
924 {
925 double teta,phi,val;
926 finalMap.PixThetaPhi(i, teta, phi);
927 int pixel = intBetaMap.PixIndexSph(teta,phi);
928 val = intBetaMap.PixVal(pixel);
929 finalMap(i) += coeff*mapToAdd(i)*val;
930 }
931 }
932 else
933 {
934 for(int i=0; i<finalMap.NbPixels();i++)
935 {
936 finalMap(i) += coeff*mapToAdd(i)*intBetaMap(i);
937 }
938 }
939 if (printlev > 1)
940 {
941 cout << "addComponentBeta(SG<T>,SG<T>,SG<T>,double) - Coeff= " << K << endl;
942 }
943}
944
945
946
947
948
949/* Nouvelle-Fonction */
950void SKM_MergeFITSKeywords(char * flnm)
951{
952 DVList dvl;
953 FitsFile::FitsExtensionType typeOfExtension;
954 int naxis;
955 vector<int> naxisn;
956 FitsFile::FitsDataType dataType;
957 FitsInFile::GetBlockType(flnm, 1, typeOfExtension, naxis, naxisn, dataType, dvl);
958 // Cleaning the keywords
959#define SZexlst 21
960 char *exlst[SZexlst]=
961 {"SIMPLE","BITPIX" ,"NAXIS" ,"NAXIS#" ,"PCOUNT","GCOUNT",
962 "EXTEND","ORIGIN" ,"DATE*" ,"TFIELDS","TTYPE#","TFORM#",
963 "TUNIT#","EXTNAME","CTYPE#","CRVAL#" ,"CRPIX#","CDELT#",
964 "XTENSION","INSTRUME","TELESCOP"};
965 char kwex[32];
966 int i,l;
967 for (i=0; i<SZexlst; i++) {
968 strncpy(kwex, exlst[i], 32);
969 l = strlen(kwex)-1;
970 if ((kwex[l] != '*') && (kwex[l] != '#')) {
971 dvl.DeleteKey(kwex);
972 }
973 else {
974 bool fgd = (kwex[l] == '#') ? true : false;
975 list<string> lstsup;
976 kwex[l] = '\0';
977 DVList::ValList::const_iterator it;
978 for (it = dvl.Begin(); it != dvl.End(); it++) {
979 if ((*it).first.substr(0,l) != kwex) continue;
980 if (fgd && !isdigit((*it).first[l])) continue;
981 lstsup.push_back((*it).first);
982 }
983 list<string>::iterator it2;
984 for (it2 = lstsup.begin(); it2 != lstsup.end(); it2++)
985 dvl.DeleteKey(*it2);
986 }
987 }
988 dvl_fitskw->Merge(dvl);
989}
990
991/* Nouvelle-Fonction */
992void SKM_MergeFITSKeywords2(DVList & dvl)
993{
994 // Cleaning the keywords
995#define SZexlst 21
996 char *exlst[SZexlst]=
997 {"SIMPLE","BITPIX" ,"NAXIS" ,"NAXIS#" ,"PCOUNT","GCOUNT",
998 "EXTEND","ORIGIN" ,"DATE*" ,"TFIELDS","TTYPE#","TFORM#",
999 "TUNIT#","EXTNAME","CTYPE#","CRVAL#" ,"CRPIX#","CDELT#",
1000 "XTENSION","INSTRUME","TELESCOP"};
1001 char kwex[32];
1002 int i,l;
1003 for (i=0; i<SZexlst; i++) {
1004 strncpy(kwex, exlst[i], 32);
1005 l = strlen(kwex)-1;
1006 if ((kwex[l] != '*') && (kwex[l] != '#')) {
1007 dvl.DeleteKey(kwex);
1008 }
1009 else {
1010 bool fgd = (kwex[l] == '#') ? true : false;
1011 list<string> lstsup;
1012 kwex[l] = '\0';
1013 DVList::ValList::const_iterator it;
1014 for (it = dvl.Begin(); it != dvl.End(); it++) {
1015 if ((*it).first.substr(0,l) != kwex) continue;
1016 if (fgd && !isdigit((*it).first[l])) continue;
1017 lstsup.push_back((*it).first);
1018 }
1019 list<string>::iterator it2;
1020 for (it2 = lstsup.begin(); it2 != lstsup.end(); it2++)
1021 dvl.DeleteKey(*it2);
1022 }
1023 }
1024 dvl_fitskw->Merge(dvl);
1025}
1026
1027/* Nouvelle-Fonction */
1028void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
1029{
1030 char *ntn[2] = {"nu","fnu"};
1031 NTuple nt(2,ntn); // Creation NTuple (AVEC new )
1032 float xnt[2];
1033 double nu;
1034 double numin = rs.minFreq();
1035 double numax = rs.maxFreq();
1036 int nmax = 500;
1037 double dnu = (numax-numin)/nmax;
1038 for(int k=0; k<nmax; k++) {
1039 nu = numin+k*dnu;
1040 xnt[0] = nu;
1041 xnt[1] = rs.flux(nu);
1042 nt.Fill(xnt);
1043 }
1044 ObjFileIO<NTuple> oiont(nt);
1045 oiont.Write(so, name);
1046 return;
1047}
1048
1049/* Nouvelle-Fonction */
1050void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name)
1051{
1052 char *ntn[2] = {"nu","tnu"};
1053 NTuple nt(2,ntn); // Creation NTuple (AVEC new )
1054 float xnt[2];
1055 double nu;
1056 double numin = sr.minFreq();
1057 double numax = sr.maxFreq();
1058 int nmax = 500;
1059 double dnu = (numax-numin)/nmax;
1060 for(int k=0; k<nmax; k++) {
1061 nu = numin+k*dnu;
1062 xnt[0] = nu;
1063 xnt[1] = sr.transmission(nu);
1064 nt.Fill(xnt);
1065 }
1066 ObjFileIO<NTuple> oiont(nt);
1067 oiont.Write(so, name);
1068 return;
1069}
1070
Note: See TracBrowser for help on using the repository browser.