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

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

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