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

Last change on this file since 887 was 886, checked in by ansari, 26 years ago

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