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

Last change on this file since 1282 was 1282, checked in by ansari, 25 years ago

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