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

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

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