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

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

Sophie: adapt to new class MapOperations
File size: 21.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 RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
32void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
33
34// to add different sky components and corresponding tools
35//----------------------------------------------------------
36template <class T>
37void addComponent(SpectralResponse& sr,
38 PixelMap<T>& finalMap,
39 PixelMap<T>& mapToAdd,
40 RadSpectra& rs, double K=1.);
41//
42template <class T>
43void addComponentBeta(SphereHEALPix<T>& finalMap,
44 SphereHEALPix<T>& mapToAdd,SpectralResponse& sr,
45 SphereHEALPix<T>& betaMap, double normFreq, double K);
46//
47template <class T>
48void integratedMap(SpectralResponse& sr,
49 SphereHEALPix<T>& betaMap, double normFreq, SphereHEALPix<T>& intBetaMap);
50
51//
52template <class T>
53void addComponentBeta(SphereHEALPix<T>& finalMap,
54 SphereHEALPix<T>& mapToAdd,
55 SphereHEALPix<T>& intBetaMap, double K);
56//
57template <class T>
58void addDipole(SpectralResponse& sr, PixelMap<T>& finalMap,
59 double theta,double phi,double amp,double temp);
60//
61// -----------------------------------------------------------------
62
63// ----- Global (static) variables ------------
64static bool rdmap = false; // true -> Read map first
65static char mapPath[256]; // Path for input maps
66static int hp_nside = 32; // HealPix NSide
67static int nskycomp = 0; // Number of sky components
68static int debuglev = 0; // Debug Level
69static int printlev = 0; // Print Level
70static POutPersist * so = NULL; // Debug PPFOut file
71
72// -------------------------------------------------------------------------
73// main program
74// -------------------------------------------------------------------------
75int main(int narg, char * arg[])
76{
77 if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
78 cout << " Usage: skymixer parameterFile outputfitsname [outppfname]" << endl;
79 exit(0);
80 }
81
82 InitTim();
83
84 string msg;
85 int rc = 0;
86 RadSpectra * es = NULL;
87 SpectralResponse * sr = NULL;
88 double moy, sig;
89
90 DataCards dc;
91 so = NULL;
92
93 try {
94 string dcard = arg[1];
95 if(printlev > 1) cout << " Decoding parameters from file " << dcard << endl;
96 dc.ReadFile(dcard);
97
98 rc = CheckCards(dc, msg);
99 if (rc) {
100 cerr << " Error condition -> Rc= " << rc << endl;
101 cerr << " Msg= " << msg << endl;
102 return(rc);
103 }
104 }
105 catch (PException exc) {
106 msg = exc.Msg();
107 cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
108 return(90);
109 }
110
111
112 cout << " skymix/Info : NComp = " << nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
113 cout << " ... MapPath = " << (string)mapPath << " DebugLev= " << debuglev
114 << " PrintLev= " << printlev << endl;
115
116 // We create an output persist file for writing debug objects
117 if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
118
119 SphereHEALPix<float> outgs(hp_nside);
120 try{
121 if (rdmap) { // Reading map from FITS file
122 char ifnm[256];
123 strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
124 ifnm[255] = '\0';
125 cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
126 {
127 FITS_SphereHEALPix<float> fios(ifnm);
128 outgs = (SphereHEALPix<float>)fios;
129 }
130 if(printlev>0)
131 cout << " Output HealPIx Map read - NbPixels= " <<
132 outgs.NbPixels() << endl;
133 if (printlev > 0) {
134 MeanSig(outgs.DataBlock(), moy, sig );
135 cout << " MeanSig for outpout map - Mean= " <<
136 moy << " Sigma= " << sig << endl;
137 }
138 }
139 else {
140 if(printlev>0)
141 cout << " Output HealPix Map created - NbPixels= " <<
142 outgs.NbPixels() << endl;
143 outgs.SetPixels(0.);
144 }
145
146 // Decoding detection pass-band filter
147 sr = getSpectralResponse(dc);
148 PrtTim(" After FilterCreation ");
149
150 char * flnm, buff[90];
151 string key;
152
153 double K = 1.;
154 double freqOfMap = 1.;
155 // Loop over sky component
156 int sk;
157 for(sk = 0; sk<nskycomp; sk++) {
158 cout << "------------------------------------" << endl;
159 cout << " Processing sky component No " << sk+1 << endl;
160
161
162 sprintf(buff, "%d", sk+1);
163 key = (string)"DIPOLE" + buff;
164 // check for an eventual dipole
165 if(dc.HasKey(key))
166 {
167 if (es) { delete es; es = NULL; }
168 double temp = -10.;
169 double theta= dc.DParam(key,1,1.);
170 double phi = dc.DParam(key,2,1.);
171 double amp = dc.DParam(key,3,1.);
172 if(dc.NbParam(key)>3)
173 {
174 temp = dc.DParam(key,4,1.);
175 }
176 cout << " creating dipole " << temp << " " << amp << " " << phi << " " << theta << " " << endl;
177 addDipole(*sr, outgs,theta,phi,amp,temp);
178 }
179 else
180 {
181 sprintf(buff, "%d", sk+1);
182 key = (string)"MAPFITSFILE" + buff;
183 flnm = BuildFITSFileName(dc.SParam(key, 0));
184
185 K = dc.DParam(key, 1, 1.);
186
187 cout << " Reading Input FITS map " << (string)flnm << endl;
188 SphereHEALPix<float> ings(hp_nside);
189 {
190 FITS_SphereHEALPix<float> fiosIn(flnm);
191 ings = (SphereHEALPix<float>)fiosIn;
192 }
193 if (debuglev > 4) { // Writing the input map to the outppf
194 FIO_SphereHEALPix<float> fiog(ings);
195 fiog.Write(*so, key);
196 }
197 if (printlev > 2) {
198 MeanSig(ings.DataBlock(), moy, sig );
199 cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
200 }
201 bool mapDependentOfFreq = false;
202 key = (string)"BETAFITSFILE"+ buff;
203 if(dc.HasKey(key))
204 {
205 mapDependentOfFreq = true;
206 }
207
208 // getting Emission spectra
209 if(!mapDependentOfFreq)
210 {
211 if (es) { delete es; es = NULL; }
212 es = getEmissionSpectra(dc, sk);
213 addComponent(*sr, outgs, ings, *es, K);
214 }
215 else
216 {
217 key = (string)"BETAFITSFILE"+ buff;
218 //SphereHEALPix<float> betaMap;
219 flnm = BuildFITSFileName(dc.SParam(key, 0));
220 double normFreq = dc.DParam(key, 1, 1.);
221 if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
222 int nSideForInt = dc.DParam(key, 2, 1.);
223 if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
224 cout << "....BetaFits... Reading Beta FITS map " << (string)flnm << endl;
225 SphereHEALPix<float> betaMap(hp_nside);
226 {
227 FITS_SphereHEALPix<float> fiosBM(flnm);
228 betaMap = (SphereHEALPix<float>)fiosBM;
229 }
230 if (debuglev > 4) { // Writing the input map to the outppf
231 FIO_SphereHEALPix<float> fiogs(betaMap);
232 fiogs.Write(*so, key);
233 }
234
235 if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
236 bool bydefault = true;
237 if(!bydefault)
238 addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
239 else
240 {
241 // integrate the betamap over the SpectralResponse
242 SphereHEALPix<float> intBetaMap(nSideForInt);
243 integratedMap(*sr, betaMap, normFreq, intBetaMap);
244 if (debuglev > 4) { // Writing the input map to the outppf
245 FIO_SphereHEALPix<float> fiogs2(intBetaMap);
246 fiogs2.Write(*so, "INTBETAMAP");
247 }
248
249 betaMap.Resize(8);
250 MeanSig(intBetaMap.DataBlock(), moy, sig );
251 if (printlev > 4)
252 cout << "....BetaFits... MeanSig for intBetaMap - Mean= "
253 << moy << " Sigma= " << sig << endl;
254 // add the integrated beta map
255 addComponentBeta(outgs,ings,intBetaMap, K);
256 }
257 }
258
259 MeanSig(outgs.DataBlock(), moy, sig );
260 cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
261 cout << "-------------------------------------------------" << endl;
262
263 sprintf(buff, "End of Processing Component %d ", sk+1);
264 PrtTim(buff);
265 }
266 }
267 }
268 catch(PException exc)
269 {
270 cout << "catched PException" << endl;
271 msg = exc.Msg();
272 cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
273 rc = 50;
274 return(50);
275 }
276
277 // Saving the output map in FITS format
278 cout << "Output Map (SphereHEALPix<float>) written to FITS file "
279 << (string)(arg[2]) << endl;
280 {
281 FITS_SphereHEALPix<float> fios2(outgs);
282 fios2.Write(arg[2]);
283 }
284 PrtTim("End of WriteFITS ");
285 // Saving the output map in PPF format
286 if (narg > 3) {
287 POutPersist s(arg[3]);
288 FIO_SphereHEALPix<float> fiog(&outgs) ;
289 fiog.Write(s);
290 cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "
291 << (string)(arg[3]) << endl;
292 PrtTim("End of WritePPF ");
293 }
294 if (so) delete so; // Closing the debug ppf file
295 return(rc);
296}
297
298/* Nouvelle-Fonction */
299int CheckCards(DataCards & dc, string & msg)
300// Function to check datacards
301{
302rdmap = false;
303mapPath[0] = '\0';
304hp_nside = 32;
305nskycomp = 0;
306debuglev = 0;
307printlev = 0;
308
309int rc = 0;
310string key, key2,key3;
311
312 // Cheking datacards
313 if (dc.NbParam("SKYMIX") < 2) {
314 rc = 71;
315 msg = "Invalid parameters - Check @SKYMIX card ";
316 return(rc);
317 }
318 key = "READMAP";
319 if (dc.HasKey(key)) {
320 if (dc.NbParam(key) < 1) {
321 rc = 72;
322 msg = "Invalid parameters - Check @READMAP card ";
323 return(rc);
324 }
325 else rdmap = true;
326 }
327
328// Checking detection filter specification
329 key = "GAUSSFILTER";
330 key2 = "FILTERFITSFILE";
331 key3 = "DIPOLE";
332 if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3) && (dc.NbParam(key3) < 3)) {
333 msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE or @DIPOLE";
334 rc = 73; return(rc);
335 }
336
337 // Decoding number of component and pixelisation parameter
338 int mg = 32;
339 int ncomp = 0;
340 ncomp = dc.IParam("SKYMIX", 0, 0);
341 mg = dc.IParam("SKYMIX", 1, 32);
342 if (ncomp < 1) {
343 msg = "Invalid parameters - Check datacards @SKYMIX ";
344 rc = 74;
345 return(rc);
346 }
347
348 // Checking detection filter specification
349 // Checking input FITS file specifications
350 int kc;
351 char buff[256];
352 bool pb = false;
353 string key4;
354 string key5;
355 string key6;
356 for(kc=0; kc<ncomp; kc++) {
357 sprintf(buff, "MAPFITSFILE%d", kc+1);
358 key = buff;
359 sprintf(buff, "DIPOLE%d", kc+1);
360 key3 = buff;
361 if (dc.NbParam(key) < 1 && dc.NbParam(key3)<1) {
362 msg = "Missing or invalid card : " + key + " " + key2 + " " + key3;
363 pb = true; break;
364 }
365 sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
366 key = buff;
367 sprintf(buff, "BLACKBODY%d", kc+1);
368 key2 = buff;
369 sprintf(buff, "POWERLAWSPECTRA%d", kc+1);
370 key3 = buff;
371 sprintf(buff, "BETAFITSFILE%d", kc+1);
372 key4 = buff;
373 sprintf(buff, "DIPOLE%d", kc+1);
374 key5 = buff;
375 sprintf(buff, "DERIVBB%d", kc+1);
376 key6 = buff;
377 if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
378 && (dc.NbParam(key6)<1) && (dc.NbParam(key5)<3)) {
379 msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
380 pb = true; break;
381 }
382
383 }
384
385 if (pb) {
386 rc = 75;
387 return(75);
388 }
389
390
391// Initialiazing parameters
392 rc = 0;
393 msg = "OK";
394 nskycomp = ncomp;
395 hp_nside = mg;
396
397// Checking for PATH definition card
398 key = "MAPPATH";
399 if (dc.NbParam(key) < 3) strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
400 mapPath[255] = '\0';
401 key = "DEBUGLEVEL";
402 debuglev = dc.IParam(key, 0, 0);
403 key = "PRINTLEVEL";
404 printlev = dc.IParam(key, 0, 0);
405 return(rc);
406}
407
408static char buff_flnm[1024]; // Mal protege !
409/* Nouvelle-Fonction */
410char* BuildFITSFileName(string const & fname)
411{
412if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
413else sprintf(buff_flnm, "%s", fname.c_str());
414return(buff_flnm);
415}
416
417/* Nouvelle-Fonction */
418SpectralResponse * getSpectralResponse(DataCards & dc)
419{
420 SpectralResponse * filt = NULL;
421
422 string key = "FILTERFITSFILE";
423 string key2 = "GAUSSFILTER";
424 string ppfname = "filter";
425
426 if (dc.HasKey(key) ) { // Reading FITS filter file
427 FitsIoServer fios;
428 char ifnm[256];
429 strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
430 ifnm[255] = '\0';
431 Matrix mtx(2,10);
432 fios.load(mtx, ifnm);
433 double numin = dc.DParam(key, 2, 1.);
434 double numax = dc.DParam(key, 3, 9999.);
435 Vector nu(mtx.NCols());
436 Vector fnu(mtx.NCols());
437 for(int k=0; k<mtx.NCols(); k++) {
438 nu(k) = mtx(0, k);
439 fnu(k) = mtx(1, k);
440 }
441 filt = new SpecRespVec(nu, fnu, numin, numax);
442 ppfname = key;
443 }
444 else if (dc.HasKey(key2) ) { // creating GaussianFilter
445 double nu0 = dc.DParam(key2, 0, 10.);
446 double s = dc.DParam(key2, 1, 1.);
447 double a = dc.DParam(key2, 2, 1.);
448 double numin = dc.DParam(key2, 3, 0.1);
449 double numax = dc.DParam(key2, 4, 9999);
450 filt = new GaussianFilter(nu0, s, a, numin, numax);
451 ppfname = key2;
452 }
453 if (filt == NULL) throw ParmError("datacard error ! No detection filter");
454 if(printlev>0)
455 {
456 cout << endl;
457 cout << " Filter decoded - Created " << endl;
458 cout << *filt << endl;
459 }
460 // for debug
461 if (debuglev > 1) SpectralResponse2Nt(*filt, *so, ppfname);
462 return(filt);
463}
464
465/* Nouvelle-Fonction */
466RadSpectra * getEmissionSpectra(DataCards & dc, int nc)
467{
468 char numb[16];
469 sprintf(numb, "%d", nc+1);
470
471 string key = (string)"SPECTRAFITSFILE" + numb;
472 string key2 = (string)"BLACKBODY" + numb;
473 string key5 = (string)"DERIVBB" + numb;
474 string key3 = (string)"POWERLAWSPECTRA" + numb;
475 string ppfname = "espectra";
476
477 RadSpectra * rs = NULL;
478 if (dc.HasKey(key) ) { // Reading emission spectra from file
479 char * ifnm = BuildFITSFileName(dc.SParam(key, 0));
480 cout << " Reading Input FITS spectra file " << (string)ifnm << endl;
481 FitsIoServer fios;
482 Matrix mtx(2,10);
483 fios.load(mtx, ifnm);
484 double numin = dc.DParam(key, 2, 1.);
485 double numax = dc.DParam(key, 3, 9999.);
486 Vector nu(mtx.NCols());
487 Vector tnu(mtx.NCols());
488 for(int k=0; k<mtx.NCols(); k++) {
489 nu(k) = mtx(0, k);
490 tnu(k) = mtx(1, k);
491 }
492 rs = new RadSpectraVec(nu, tnu, numin, numax);
493 ppfname = key;
494 }
495 else if (dc.HasKey(key2) ) { // Creating BlackBody emission spectra
496 rs = new BlackBody(dc.DParam(key2, 0, 2.726));
497 ppfname = key2;
498 }
499 else if (dc.HasKey(key5) ) { // Creating Dipole
500 rs = new DerivBlackBody(dc.DParam(key5, 0, 3.E-3));
501 ppfname = key5;
502 }
503 else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
504 double a = dc.DParam(key3, 0, 1.);
505 double nu0 = dc.DParam(key3, 1, 100.);
506 double dnu = dc.DParam(key3, 2, 10.);
507 double b = dc.DParam(key3, 3, 0.);
508 double numin = dc.DParam(key3, 4, 0.1);
509 double numax = dc.DParam(key3, 5, 9999);
510 rs = new PowerLawSpectra(a, b, nu0, dnu, numin, numax);
511 ppfname = key3;
512 }
513 if (rs == NULL) throw ParmError("datacard error ! missing Emission spectra");
514 cout << " Emission spectra decoded - Created (" << ppfname << ")" << endl;
515 cout << *rs << endl;
516// for debug
517 if (debuglev > 2) RadSpec2Nt(*rs, *so, ppfname);
518 return(rs);
519}
520
521
522
523
524template <class T>
525void addDipole(SpectralResponse& sr, PixelMap<T>& finalMap,
526 double theta,double phi,double amp,double temp)
527{
528 DerivBlackBody dbb;
529 if(temp>0) dbb.setTemperature(temp);
530 double coeff = dbb.filteredIntegratedFlux(sr) * amp;
531 UnitVector vd(theta,phi);
532 UnitVector vc(theta,phi);
533
534 for(int i=0; i<finalMap.NbPixels(); i++)
535 {
536 double thetar,phir;
537 finalMap.PixThetaPhi(i,thetar,phir);
538 vc.SetThetaPhi(thetar, phir);
539 finalMap(i) += vd.Psc(vc)*coeff;
540 }
541 if (debuglev > 4) { // Writing the input map to the outppf
542 SphereHEALPix<float> ings(sqrt((double)finalMap.NbPixels()/12));
543 for(int i=0; i<finalMap.NbPixels(); i++)
544 {
545 double thetar,phir;
546 finalMap.PixThetaPhi(i,thetar,phir);
547 vc.SetThetaPhi(thetar, phir);
548 ings(i) = vd.Psc(vc);
549 }
550 FIO_SphereHEALPix<float> fiog(ings);
551 fiog.Write(*so, "dipole");
552 cout << "Debug the dipole map....saved in debug file !" << endl;
553 }
554}
555/* Nouvelle-Fonction */
556template <class T>
557void addComponent(SpectralResponse& sr, PixelMap<T>& finalMap,
558 PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
559{
560 // finalMap = finalMap + coeff* mapToAdd
561 // coeff = convolution of sr and rs
562 // compute the coefficient corresponding to mapToAdd
563 if (finalMap.NbPixels() != mapToAdd.NbPixels())
564 throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
565 double coeff = rs.filteredIntegratedFlux(sr) * K;
566 if (printlev > 1)
567 cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl;
568 for(int i=0; i<finalMap.NbPixels(); i++)
569 {
570 finalMap(i) += coeff * mapToAdd(i);
571 }
572}
573/* Nouvelle-Fonction */
574template <class T>
575void addComponentBeta(SphereHEALPix<T>& finalMap,
576 SphereHEALPix<T>& mapToAdd,SpectralResponse& sr,
577 SphereHEALPix<T>& betaMap, double normFreq, double K)
578{
579 // finalMap = finalMap + coeff* mapToAdd
580 // coeff = convolution of sr and rs
581 // compute the coefficient corresponding to mapToAdd
582 // betaMap is the map of (beta(theta,phi))
583
584 int nbpix = finalMap.NbPixels();
585 if (nbpix != mapToAdd.NbPixels())
586 throw SzMismatchError("addComponentBeta()/Error: Unequal number of Input/Output map pixels");
587 if (printlev > 1)
588 {
589 cout << "addComponentBeta - Coeff= " << K << endl;
590 cout << "nb pixels: " << finalMap.NbPixels() << endl;
591 }
592 SphereHEALPix<T> bigBetaMap(sqrt((double)nbpix/12));
593 if(nbpix != betaMap.NbPixels())
594 {
595 Sph2Sph(betaMap,bigBetaMap);
596 }
597 for(int i=0; i<finalMap.NbPixels(); i++)
598 {
599 // coeff = integration of (nu/normFreq)^(-beta(theta,phi)) in each pixels
600 RadSpectra* rs = new PowerLawSpectra
601 (1.,-bigBetaMap(i), 0., normFreq, 0.01, 800.);
602 double coeff = rs->filteredIntegratedFlux(sr);
603 finalMap(i) += coeff*K*mapToAdd(i);
604 }
605}
606
607template <class T>
608void integratedMap(SpectralResponse& sr,
609 SphereHEALPix<T>& betaMap,
610 double normFreq,
611 SphereHEALPix<T>& intBetaMap)
612{
613 PowerLawSpectra rs(1.,-2., 0., normFreq);
614
615 if(betaMap.NbPixels()!=intBetaMap.NbPixels())
616 {
617 Sph2Sph(betaMap,intBetaMap);
618 for(int i=0; i<intBetaMap.NbPixels(); i++)
619 {
620 rs.setExp(-intBetaMap(i));
621 double coeff = rs.filteredIntegratedFlux(sr);
622 intBetaMap(i) = coeff;
623 }
624 }
625 else
626 {
627 for(int i=0; i<intBetaMap.NbPixels(); i++)
628 {
629 rs.setExp(-betaMap(i));
630 double coeff = rs.filteredIntegratedFlux(sr);
631 intBetaMap(i) = coeff;
632 }
633 }
634}
635
636template <class T>
637void addComponentBeta(SphereHEALPix<T>& finalMap,
638 SphereHEALPix<T>& mapToAdd,SphereHEALPix<T>& intBetaMap, double K)
639{
640 // finalMap = finalMap + coeff* mapToAdd
641 // coeff = convolution of sr and rs
642 // compute the coefficient corresponding to mapToAdd
643 // integBetaMap is the map of the integration (nu/normFreq)^(-beta(theta,phi)) over
644 // the spectralResponse
645 // different from addComponentBeta(PixelMap<T>& finalMap,
646 // PixelMap<T>& mapToAdd,SpectralResponse& sr, PixelMap<T>& betaMap, double normFreq, double K)
647 // since it permits to use a intBetaMap with a different number of pixels than
648 // the other maps
649
650 int nbpix = finalMap.NbPixels();
651 if (nbpix != mapToAdd.NbPixels())
652 throw SzMismatchError("addComponentBeta(PixelMap<T>&,PixelMap<T>&,PixelMap<T>&,double)/Error: Unequal number of Input/Output map pixels");
653 double coeff = K;
654
655 if(nbpix != intBetaMap.NbPixels())
656 {
657 for(int i=0; i<finalMap.NbPixels();i++)
658 {
659 double teta,phi,val;
660 finalMap.PixThetaPhi(i, teta, phi);
661 int pixel = intBetaMap.PixIndexSph(teta,phi);
662 val = intBetaMap.PixVal(pixel);
663 finalMap(i) += coeff*mapToAdd(i)*val;
664 }
665 }
666 else
667 {
668 for(int i=0; i<finalMap.NbPixels();i++)
669 {
670 finalMap(i) += coeff*mapToAdd(i)*intBetaMap(i);
671 }
672 }
673 if (printlev > 1)
674 {
675 cout << "addComponentBeta(SG<T>,SG<T>,SG<T>,double) - Coeff= " << K << endl;
676 }
677}
678
679
680
681
682
683/* Nouvelle-Fonction */
684
685/* Nouvelle-Fonction */
686void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
687{
688 char *ntn[2] = {"nu","fnu"};
689 NTuple nt(2,ntn); // Creation NTuple (AVEC new )
690 float xnt[2];
691 double nu;
692 double numin = rs.minFreq();
693 double numax = rs.maxFreq();
694 int nmax = 500;
695 double dnu = (numax-numin)/nmax;
696 for(int k=0; k<nmax; k++) {
697 nu = numin+k*dnu;
698 xnt[0] = nu;
699 xnt[1] = rs.flux(nu);
700 nt.Fill(xnt);
701 }
702 ObjFileIO<NTuple> oiont(nt);
703 oiont.Write(so, name);
704 return;
705}
706
707/* Nouvelle-Fonction */
708void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name)
709{
710 char *ntn[2] = {"nu","tnu"};
711 NTuple nt(2,ntn); // Creation NTuple (AVEC new )
712 float xnt[2];
713 double nu;
714 double numin = sr.minFreq();
715 double numax = sr.maxFreq();
716 int nmax = 500;
717 double dnu = (numax-numin)/nmax;
718 for(int k=0; k<nmax; k++) {
719 nu = numin+k*dnu;
720 xnt[0] = nu;
721 xnt[1] = sr.transmission(nu);
722 nt.Fill(xnt);
723 }
724 ObjFileIO<NTuple> oiont(nt);
725 oiont.Write(so, name);
726 return;
727}
728
Note: See TracBrowser for help on using the repository browser.