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

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

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