source: Sophya/trunk/SigPredictor/TestSources.cc@ 3736

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

mise a jour
File size: 15.8 KB
Line 
1 // Dominique YVON, CEA/DAPNIA/SPP 02/2000
2
3#include <stdio.h>
4#include <stdlib.h>
5#include <math.h>
6#include <iostream>
7#include <fstream>
8//#include <SIOUX.h>
9
10#ifdef __MWERKS__
11 #include "unixmac.h"
12 #include "macenvvariables.h"
13#endif
14
15#include "squarefilt.h"
16
17#include "alllobe.h"
18#include "alllightsources.h"
19#include "sigcalctools.h"
20
21#define NLatTestMap (256)
22
23int CompareMapResults();
24int TestSommeMapsInBand(char file1[], long nlat);
25
26int CompareMapResults()
27{
28 cout<< "Test du fonctionnement par la comparaison systematique";
29 cout<< " Avec les resultats de supertango"<<endl;
30
31
32
33#ifndef __MWERKS__
34 char* PATHResults=getenv("PATHResults");
35 char* PATHDataLScr=getenv("PATHDataLScr");
36 char* PATHCarteLobe=getenv("PATHCarteLobe");
37 char* PATHSTangoRes=getenv("PATHSTangoRes");
38#endif
39
40// double Min, Max, Moy, sigma;
41 char filename [150];
42
43 sprintf(filename,"CarteCieltot_res0256.fits");
44
45 TestSommeMapsInBand(filename,NLatTestMap);
46
47 return 0;
48}
49
50int TestSommeMapsInBand(char fileTotPower[], long nlat) // Pas fini Ne pas utiliser!
51{
52 double Min, Max, Moy, sigma;
53 char filename[150];
54 // On lit les donnees brutes
55 SphereHEALPix<r_4> PowerMap(nlat); // Pour se mettre a l'unite nW/m2/St
56
57 // Definition channel C
58 LobeGaussien LobeGaussChanC(28.2/60.,30.e9,50.e09);
59 // LobeConique ConeLobe(1.,10.e9,3000.e9); // pour des test simples
60 SquareFilter FlatFilter(36.e9,44.e9); // Hz
61
62 // Test lobes
63 // On somme les trois contributions intŽgrŽes sur la bande passante
64/*
65 // OK pour le synchrotron
66 { LightSynchro Synchro(nlat);
67 SigCalcTool ChanCTool(&Synchro,&LobeGaussChanC,&FlatFilter);
68 addInInBandPowerMap(PowerMap,ChanCTool); // OK pour le synchrotron
69 }
70
71 // On ajoute la poussiere diffuse
72 { LightDiffDust DiffDust(nlat);
73 SigCalcTool ChanCTool(&DiffDust,&LobeGaussChanC,&FlatFilter);
74 addInInBandPowerMap(PowerMap,ChanCTool);
75 }
76
77 // On ajoute le CMB
78 { LightCMBPrim CMB(nlat);
79 SigCalcTool ChanCTool(&CMB,&LobeGaussChanC,&FlatFilter);
80 addInInBandPowerMap(PowerMap,ChanCTool);
81 }
82
83
84 // On ajoute du bruit
85 { SphereHEALPix<r_4> NoiseChanCMap(nlat);
86 sprintf(filename, "%schannelC_noise_res%04i.fits",PATHSTangoRes,nlat);
87 FitsServer.load(NoiseChanCMap,filename);
88 scaleMap(1.e-9,NoiseChanCMap); // Pour se mettre W/m2/St
89 addMap(PowerMap,NoiseChanCMap);
90 }
91
92 // Bilan carte sommee
93 cout<<"Carte Totale"<<endl;
94 MinMaxSigMap(PowerMap,Min,Max,Moy,sigma);
95 cout <<"PowerMap: "<<"Min="<<Min<<" Max="<<Max;
96 cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
97
98 FitsServer.save(PowerMap,file1);
99*/
100/*
101 // Test-Validation Ptotale si ca se passe bien
102 {
103 SphereHEALPix<r_4> ObsTotSTango(256);
104 sprintf(filename, "%schannelC_obs_res0256.fits",PATHSTangoRes);
105 FitsServer.load(ObsTotSTango,filename);
106 scaleMap(1.e-9,ObsTotSTango); // Pour se mettre W/m2/St
107
108 MinMaxSigMap(ObsTotSTango,Min,Max,Moy,sigma);
109 cout <<"ObsTotSTango W/m2/sr: "<<"Min="<<Min<<" Max="<<Max;
110 cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
111
112 divMap1WithMap2(PowerMap,ObsTotSTango);
113 MinMaxSigMap(PowerMap,Min,Max,Moy,sigma);
114 cout <<"Test: "<<"Min="<<Min<<" Max="<<Max;
115 cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
116
117 // On sauve la carte sur disque
118 sprintf(filename,"SinusDivCieltot_res0256.fits");
119 FitsServer.sinus_picture_projection(PowerMap,filename);
120
121 // On regarde la carte obtenue
122 sprintf(filename,"SinusCieltot_res%04i.fits",nlat);
123 FitsServer.sinus_picture_projection(PowerMap,filename);
124 }
125*/
126
127/*
128// On veut convoluer avec un lobe gaussien
129 {
130
131 LightSrcMapPowerInband LightMapInCBand(fileTotPower,256,36.e9,44.e9);
132 // On utilise une carte deja integree sur la reponse spectrale
133 // Dont les bornes en frequence sont specifiees dans le constructeur
134 // Pour la forme en coherence avec ChannelC
135 // On utilse un lobe (large bande)
136
137
138 // Voila l'outil
139 SigCalcTool ChannelCtool(&LightMapInCBand,&LobeGaussChanC,&FlatFilter);
140
141 // La carte
142 SphereHEALPix<r_4> ChannelCConvolue(nlat);
143
144 // On convolue
145 addToSkyMap(ChannelCConvolue, ChannelCtool);
146
147 // On sauvegarde
148 sprintf(filename,"ChannelConvolue%04i.fits",nlat);
149 FitsServer.save(ChannelCConvolue, filename);
150
151 MinMaxSigMap(ChannelCConvolue,Min,Max,Moy,sigma);
152 cout <<"ChannelCConvolue: "<<"Min="<<Min<<" Max="<<Max;
153 cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
154
155 // On regarde la carte obtenue
156 sprintf(filename,"SinusCielCconvol%04i.fits",nlat);
157 FitsServer.sinus_picture_projection(ChannelCConvolue,filename);
158 }
159 */
160
161 /*
162 { // On Compare
163#ifndef __MWERKS__
164 char* PATHSTangoRes=getenv("PATHSTangoRes");
165#endif
166 SphereHEALPix<r_4> ObsConvolSTango(nlat);
167 sprintf(filename, "%schannelC_convolved.fits",PATHSTangoRes);
168 FitsServer.load(ObsConvolSTango,filename);
169 // Unite nW/m2/St
170
171 // On met a la meme echelle
172 LightSrcMapPowerInband LightMapInCBand(fileTotPower,256,36.e9,44.e9);
173 SigCalcTool ChannelCtool(&LightMapInCBand,&LobeGaussChanC,&FlatFilter);
174 double lobeSize=ChannelCtool.CalcLobeSize(); //Steradian
175 double scaleFactor= 1.e-9 * lobeSize; // Pour se mettre W/m2
176 scaleMap(scaleFactor,ObsConvolSTango);
177
178 MinMaxSigMap(ObsConvolSTango,Min,Max,Moy,sigma);
179 cout <<"ObsConvolSTango W/m2/sr: "<<"Min="<<Min<<" Max="<<Max;
180 cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
181
182 SphereHEALPix<r_4> ChannelCConvolue(nlat);
183 sprintf(filename,"ChannelConvolue%04i.fits",nlat);
184 FitsServer.load(ChannelCConvolue,filename);
185
186 MinMaxSigMap(ChannelCConvolue,Min,Max,Moy,sigma);
187 cout <<"ChannelCConvolue W/m2/sr: "<<"Min="<<Min<<" Max="<<Max;
188 cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
189
190 // On divise
191 divMap1WithMap2(ChannelCConvolue,ObsConvolSTango);
192 MinMaxSigMap(ChannelCConvolue,Min,Max,Moy,sigma);
193 cout <<"Test: "<<"Min="<<Min<<" Max="<<Max;
194 cout<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
195
196 // On sauvegarde la carte
197 sprintf(filename,"SinusConvolDiv%04i.fits",nlat);
198 FitsServer.sinus_picture_projection(ChannelCConvolue,filename);
199 }
200*/
201 return 0;
202}
203
204
205
206 /*
207 // On lit les templates pour voir
208 sprintf(filename, "%sdust_res0256.fits",PATHDataLScr);
209 FitsServer.load(DiffDustMap,filename);
210 MinMaxSigMap(DiffDustMap,Min,Max,Moy,sigma);
211 flog<<"DiffDustMap: "<<"Min="<<Min<<" Max="<<Max;
212 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
213
214 */
215 /*
216 sprintf(filename, "%ssync_res0256.fits",PATHDataLScr);
217 FitsServer.load(SynchroMap,filename);
218 MinMaxSigMap(SynchroMap,Min,Max,Moy,sigma);
219 flog<<"SynchroMap: "<<"Min="<<Min<<" Max="<<Max;
220 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
221 */
222
223/*
224 sprintf(filename, "%scmb_res0256.fits",PATHDataLScr);
225 FitsServer.load(CMBPrimMap,filename);
226 MinMaxSigMap(CMBPrimMap,Min,Max,Moy,sigma);
227 flog<<"CMBPrimMap: "<<"Min="<<Min<<" Max="<<Max;
228 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
229 */
230
231 /*
232 // Test Module synchrotron
233 SphereHEALPix<r_4> ChanCSynchro(256);
234 SphereHEALPix<r_4> MYCSynchro(256);
235
236
237 // Sans lobes d'abord
238 // Lecture resultats de supertango
239 sprintf(filename, "%schannelC_sync_res0256.fits",PATHSTangoRes);
240 FitsServer.load(ChanCSynchro,filename);
241 MinMaxSigMap(ChanCSynchro,Min,Max,Moy,sigma);
242 flog<<"ChanCSynchro: "<<"Min="<<Min<<" Max="<<Max;
243 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
244
245 // Calcul SigPred
246 LightSynchro Synchro(256);
247 SigCalcTool ChanCTool(&Synchro,&ConeLobe,&FlatFilter);
248 //ChanCTool.SetLightScr(&Synchro);
249
250 addInBandPower(MYCSynchro,ChanCTool);
251 scaleMap(1.e9,MYCSynchro); // Pour se mettre a l'unite nW/m2/St
252 MinMaxSigMap(MYCSynchro,Min,Max,Moy,sigma);
253 flog<<"MYCSynchro: "<<"Min="<<Min<<" Max="<<Max;
254 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
255*/
256 // Validation
257/* //substractMap(MYCSynchro,ChanCSynchro);
258 divMap1WithMap2(MYCSynchro,ChanCSynchro);
259 MinMaxSigMap(MYCSynchro,Min,Max,Moy,sigma);
260 flog<<"Division des cartes: "<<"Min="<<Min<<" Max="<<Max;
261 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
262*/
263/* // test
264 divMap1WithMap2(MYCSynchro,SynchroMap);
265 MinMaxSigMap(MYCSynchro,Min,Max,Moy,sigma);
266 flog<<"MycoeffIntegrale MYCSynchro: "<<"Min="<<Min<<" Max="<<Max;
267 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
268
269 divMap1WithMap2(ChanCSynchro,SynchroMap);
270 MinMaxSigMap(ChanCSynchro,Min,Max,Moy,sigma);
271 flog<<"STangoCoeffIntegrale ChanCSynchro: "<<"Min="<<Min<<" Max="<<Max;
272 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
273
274*/
275
276/*
277 // Test Module CMB
278 SphereHEALPix<r_4> ChanCCMB(256);
279 SphereHEALPix<r_4> MYCCMB(256);
280
281 // Sans lobes d'abord: Integrale sur les dependances spectrales
282 //degrŽs, freqmin, Freqmax
283 LightCMBPrim CMB(256);
284 SigCalcTool ChanCTool(&CMB,&ConeLobe,&FlatFilter);
285
286 // Lecture resultats de supertango
287 sprintf(filename,"%schannelC_cmb_res0256.fits",PATHSTangoRes);
288 FitsServer.load(ChanCCMB,filename);
289 MinMaxSigMap(ChanCCMB,Min,Max,Moy,sigma);
290 flog<<"ChanCCMB: "<<"Min="<<Min<<" Max="<<Max;
291 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
292*/
293/* // Calcul SigPred
294 addInBandPower(MYCCMB,ChanCTool);
295 scaleMap(1.e9,MYCCMB); // Pour se mettre a l'unite nW/m2/St
296 MinMaxSigMap(MYCCMB,Min,Max,Moy,sigma);
297 flog<<"MYCCMB: "<<"Min="<<Min<<" Max="<<Max;
298 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
299*/
300 // Validation
301
302/* // test
303 divMap1WithMap2(MYCCMB,CMBPrimMap);
304 MinMaxSigMap(MYCCMB,Min,Max,Moy,sigma);
305 flog<<" MycoeffIntegrale MYCCMB: "<<"Min="<<Min<<" Max="<<Max;
306 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
307
308 divMap1WithMap2(ChanCCMB,CMBPrimMap);
309 MinMaxSigMap(ChanCCMB,Min,Max,Moy,sigma);
310 flog<<"STangoCoeffIntegrale ChanCCMB: "<<"Min="<<Min<<" Max="<<Max;
311 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
312*/
313
314/*
315 // Sans lobes d'abord: Integrale sur les dependances spectrales
316 //degrŽs, freqmin, Freqmax
317 LightDiffDust DiffDust(256);
318 SigCalcTool ChanCTool(&DiffDust,&ConeLobe,&FlatFilter);
319
320 // Test Module poussire diffuse
321 SphereHEALPix<r_4> ChanCDiffDust(256);
322 SphereHEALPix<r_4> MYCDiffDust(256);
323
324
325 // Lecture resultats de supertango
326 sprintf(filename,"%schannelC_dust_res0256.fits",PATHSTangoRes);
327 FitsServer.load(ChanCDiffDust,filename);
328 MinMaxSigMap(ChanCDiffDust,Min,Max,Moy,sigma);
329 flog<<"ChanCDiffDust: "<<"Min="<<Min<<" Max="<<Max;
330 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
331*/
332
333/* // Calcul SigPred
334 addInBandPower(MYCDiffDust,ChanCTool);
335 scaleMap(1.e9,MYCDiffDust); // Pour se mettre a l'unite nW/m2/St
336 MinMaxSigMap(MYCDiffDust,Min,Max,Moy,sigma);
337 flog<<"MYCDiffDust: "<<"Min="<<Min<<" Max="<<Max;
338 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
339*/
340 // Validation
341
342/* // test
343 divMap1WithMap2(MYCDiffDust,DiffDustMap);
344 MinMaxSigMap(MYCDiffDust,Min,Max,Moy,sigma);
345 flog<<"MycoeffIntegrale MYCDiffDust: "<<"Min="<<Min<<" Max="<<Max;
346 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl;
347
348 divMap1WithMap2(ChanCDiffDust,DiffDustMap);
349 MinMaxSigMap(ChanCDiffDust,Min,Max,Moy,sigma);
350 flog<<"STangoCoeffIntegrale ChanCDiffDust: "<<"Min="<<Min<<" Max="<<Max;
351 flog<<" Moy="<<Moy<<" sigma="<<sigma<<endl<<endl;
352*/
353
354
355int Basictests()
356{
357// MakeTrap99Timeline();
358// cerr.setf(ios::scientific);
359 cerr<< "C'est parti"<<endl;
360
361
362
363 SphereHEALPix<r_4>* pTestMap;
364 pTestMap= new SphereHEALPix<r_4>(64);
365 SphereHEALPix<r_4>& TestMap=(*pTestMap);
366
367#ifndef __MWERKS__
368 char* PATHResults=getenv("PATHResults");
369 char* PATHDataLScr=getenv("PATHDataLScr");
370 char* PATHCarteLobe=getenv("PATHCarteLobe");
371#endif
372
373//_________________ Test implementation des lobes ___________________
374 LobeGaussien GaussLobe(0.3,75.e9,125.e9); //degrŽs, freqmin, Freqmax
375 LobeConique ConeLobe(1.,75.e9,125.e9);
376 Lobe4PiGaussien Lobe4Pi1(1.,75.e9,125.e9);
377
378 // Test de la classe Cartelobe
379 char lobeFileNameRoot[150];
380 sprintf(lobeFileNameRoot, "%sdetector100_1_cf",PATHCarteLobe);
381 char LobeFitsFile[150];
382 sprintf(LobeFitsFile, "%sCarteLobe100GHz_1_cf.fits",PATHResults);
383
384
385// Test LobeCartoMoyen ___________________________________________
386// CarteLobe CarteLobe100GHzCF(lobeFileNameRoot,100.e9); // Hz
387// LobeCartoMoyen LobeCartoTest(&CarteLobe100GHzCF,75.e9,125.e9);
388
389 // Un filtre bete........
390 SquareFilter FlatFilter(85.e9, 115.e9); // Hz
391
392 cout<<"param specifiques a BoloFilterFlat"<<endl;
393 cout<<"mean: "<< FlatFilter.meanFreq()<<"/tPeakTransmission: "<<
394 FlatFilter.transmission(FlatFilter.meanFreq())<<flush;
395 cout<<"Width: "<<FlatFilter.maxFreq()-FlatFilter.minFreq()<<endl;
396/*
397 cout<<" Lobe conique:"<<endl;
398 cout<<" Nlat gorsky: "<<LobeCartoTest.gorskyNlatFromRes();
399 cout<<" Resolution du lobe: "<<LobeCartoTest.lobeResol();
400 cout<<" Radian" <<endl;
401 */
402 char FileSinus[150];
403
404 // Du plus simple au plus compliquŽ
405/*
406 // Dipole
407 LightDipole Dipole;
408 SigCalcTool DipoleTool(&Dipole,&ConeLobe,&FlatFilter);
409// SigCalcTool DipoleTool(&Dipole,&GaussLobe,&FlatFilter);
410 compSkyMap(TestMap,DipoleTool);
411 sprintf(FileSinus,"%sMapSinus_Dipole.fits",PATHResults);
412 FitsServer.sinus_picture_projection(TestMap,FileSinus);
413
414 double MaxDipole=0;
415 for (long i=0; i<TestMap.NbPixels(); i++)
416 if(TestMap(i)>MaxDipole) MaxDipole=TestMap(i);
417 cout<<"Maximum du dipole: "<<MaxDipole<<endl;
418*/
419 /*
420 // Quadrupole
421 LightQuadrupole Quadrupole;
422 SigCalcTool QuadTool(&Quadrupole,&ConeLobe,&FlatFilter);
423 compSkyMap(TestMap,QuadTool);
424 sprintf(FileSinus,"%sMapSinus_quadrupole.fits",PATHResults);
425 FitsServer.sinus_picture_projection(TestMap,FileSinus);
426
427 // Total des contributions du dipole
428 addToSkyMap(TestMap,DipoleTool);
429 sprintf(FileSinus, "%sMapSinus_Dip_Quad.fits",PATHResults);
430 FitsServer.sinus_picture_projection(TestMap,FileSinus);
431*/
432
433/* // Effet du soleil
434 double SunAngsize = 4.833e-3; // Radian;
435 // double SunAngsize = 4.833e-2; // Debug
436 QuasiPtSources Sun(5800.,SunAngsize,M_PI/2.,M_PI);
437 SigCalcTool SunTool(&Sun,&GaussLobe,&FlatFilter);
438 //SigCalcTool SunTool(&Sun,&ConeLobe,&FlatFilter);
439 compSkyMap(TestMap,SunTool);
440
441// sprintf(FileSinus,"%sSphereGorsky_Soleil.fits",PATHResults);
442// FitsServer.save(TestMap, FileSinus);
443
444 sprintf(FileSinus,"%sMapSinus_Soleil.fits",PATHResults);
445 FitsServer.sinus_picture_projection(TestMap,FileSinus);
446
447 double MaxSun=0.;
448 for (long i=0; i<TestMap.NbPixels(); i++)
449 if(TestMap(i)>MaxSun) MaxSun=TestMap(i);
450 cout<<"Maximum du soleil: "<<MaxSun<<endl;
451*/
452
453/*
454 // Fluctuation primordiales
455 LightCMBPrim CMBFlucLight(256,1.e-5);
456 SigCalcTool CMBGraalTool(&CMBFlucLight,&ConeLobe,&FlatFilter);
457 compSkyMap(TestMap,CMBGraalTool);
458 sprintf(FileSinus,"%sMapSinus_CMBPrimordial.fits",PATHResults);
459 FitsServer.sinus_picture_projection(TestMap,FileSinus);
460*/
461
462/*
463 // Carte synchrotron
464 LightSynchro Synchro(256);
465 SigCalcTool SynchroTool(&Synchro,&ConeLobe,&FlatFilter);
466 compSkyMap(TestMap,SynchroTool);
467
468 sprintf(FileSinus,"%sMapSinus_Synchro.fits",PATHResults);
469 FitsServer.sinus_picture_projection(TestMap,FileSinus);
470*/
471
472/*
473 // Diffuse Dust
474 LightDiffDust DiffDust(256);
475 SigCalcTool DiffDustTool(&DiffDust,&ConeLobe,&FlatFilter);
476 compSkyMap(TestMap,DiffDustTool);
477 sprintf(FileSinus,"%sMapSinus_DiffDust.fits",PATHResults);
478 FitsServer.sinus_picture_projection(TestMap,FileSinus);
479
480*/
481
482/*
483 LightGalaxResol Galaxresol(128);
484 SigCalcTool GalacResolTool(&Galaxresol,&ConeLobe,&FlatFilter);
485 compSkyMap(TestMap,GalacResolTool);
486 sprintf(FileSinus,"%sMapSinus_GalaxResol.fits",PATHResults);
487 FitsServer.sinus_picture_projection(TestMap,FileSinus);
488*/
489
490/*
491 // Fond galactique non resolues
492 LightGalaxNoResol GalaxNoresol(128);
493 SigCalcTool GalacNoResolTool(&GalaxNoresol,&ConeLobe,&FlatFilter);
494 compSkyMap(TestMap,GalacNoResolTool);
495 sprintf(FileSinus,"%sMapSinus_GalaxNOResol.fits",PATHResults);
496 FitsServer.sinus_picture_projection(TestMap,FileSinus);
497
498*/
499
500return 0;
501}
Note: See TracBrowser for help on using the repository browser.