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

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

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