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

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

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