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

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

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