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

Last change on this file since 975 was 798, checked in by ansari, 26 years ago

Creation du module SigPredictor (Simulation de signal Archeops/Planck)

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