source: Sophya/trunk/SigPredictor/sigcalctools.cc@ 1147

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

Fichiers au format unix

dominique
File size: 16.2 KB
Line 
1 // Dominique YVON, CEA/DAPNIA/SPP 02/2000
2
3#include <math.h>
4#include <iostream>
5#include <iostream>
6#include <fstream>
7#ifdef __MWERKS__
8 #include "mwerksmath.h"
9 #include "unixmac.h"
10#endif
11#include "sigcalctools.h"
12#include "lightdipole.h"
13
14//_______________ ici toutes les frequences sont en Hz ___________________________
15
16static SigCalcTool* pSigToolcur;
17
18double SigCalGLFreqFunc1(double freq) {
19 double temp1=(pSigToolcur->pLSrc)->spectre(freq);
20 double temp2=(pSigToolcur->pLobe)->spectre(freq);
21 double temp3=(pSigToolcur->pFilter)->transmission(freq);
22
23 return temp1*temp2*temp3;
24}
25
26double SigCalGLFreqFunc2(double freq)
27{
28 // Integration function for GLInteg
29 double temp1=
30 (pSigToolcur->pLSrc)->powSpecDens((pSigToolcur->VPointe).Theta(),(pSigToolcur->VPointe).Phi(),freq);
31 double temp2=(pSigToolcur->pLobe)->weigth(pSigToolcur->VCur,pSigToolcur->VPointe,pSigToolcur->VY,freq);
32 double temp3=(pSigToolcur->pFilter)->transmission(freq);
33 return temp1*temp2*temp3;
34}
35
36SigCalcTool::SigCalcTool(AbsLightSource* pLightSrc, AbsLobeNoPolar* pLobeNoPolar,
37 SpectralResponse* pFilt):pLSrc(pLightSrc),pLobe(pLobeNoPolar),pFilter(pFilt)
38{
39 SigCalcToolInit();
40}
41
42void SigCalcTool::SigCalcToolInit()
43{ emptySignal=false;
44// Compute frequency integration boundaries
45 cout<< "Initialisation Calctool"<<endl;
46 FreqMin=max(pLobe->minFreq(), pFilter->minFreq());
47 FreqMax=min(pLobe->maxFreq(), pFilter->maxFreq());
48 if(FreqMax<FreqMin) {
49 emptySignal=true;
50 cerr<< "Frequency max is lower than Frequency Min in SigCalcTool"<<endl;
51 cerr<< "check consistency of lobes and Filters"<<endl;
52 }
53// Computation Options
54 if(pLSrc->IsMappedPowerSrc())
55 { if(!pLobe->IsFreqSep())
56 { cerr<<" Sigcalctool error: using a LightMapPowerInband with a lobe non freq separable"<<endl;
57 cerr<<" Did you change lobe between constructing the map and running sigcalctool?"<<endl;
58 cerr<<" Program exited"<<endl;
59 exit(-1.);
60 }
61 Option=IsLightMapPowerInband;
62 pIntegrale= new GLInteg();
63 // Pour eviter un plantage dans ~SigCalcTool
64 }
65
66 else if(pLSrc->IsFreqSep()&&pLobe->IsFreqSep()) {
67 Option=AllSeparable;
68 pIntegrale= new GLInteg(SigCalGLFreqFunc1,FreqMin,FreqMax); //en Hz.
69 pSigToolcur=this;
70 pIntegrale->NStep(200); // Integration tres sŽrieuse
71 IntegSpectOverFreq=pIntegrale->Value();
72 }
73
74 else
75 { Option=NonSeparable;
76 pIntegrale= new GLInteg(SigCalGLFreqFunc2,FreqMin,FreqMax);
77 pIntegrale->NStep(10); // Pour aller plus vite. Serieux si le filtre est "compact"
78 }
79// Computation Resolution
80 RAngComp=pLSrc->LSrcResol(); // On integre sur la resolution de la carte
81 if(RAngComp==0.)
82 { RAngComp=pLobe->lobeResol();
83 if(RAngComp==0.)
84 { cerr<<" Bizarre un lobe de resolution nulle?"<<endl;
85 RAngComp= 5.e-4; // Radians
86 // On prend la resolution nominale de Planck
87 }
88 }
89 if(RAngComp<pLobe->lobeResol())
90 { cerr<<" SigCalcTool: LightSource resolution lower than expected lobe resolution"<<endl;
91 cerr<<" Not healthy: Ckeck consistency"<<endl;
92 }
93 cout<<"Resolution de calcul: "<<RAngComp<<" Radian"<<endl<<endl;
94}
95
96double SigCalcTool::compPixel(double theta,double phi) {
97 UnitVector VP(theta,phi);
98 UnitVector VYbidon=VP.VperpPhi();
99// Compute unit vector perpendicular to Vpoin at same theta
100 return compPixel(VP,VYbidon);
101}
102
103double SigCalcTool::compPixel(UnitVector& VP, UnitVector& VdirectY){
104 double returnRes=0.;
105 VPointe=VP;
106 VY=VdirectY;
107 VX=VY^VP;
108 if(!emptySignal) returnRes=powerInteg(); // On integre sur la sphere
109 return returnRes;
110}
111
112
113double SigCalcTool::calcPowerDens() const{
114// Compute the power integrated on frequency dependance, (Lobe and LightSource)
115 pSigToolcur=(SigCalcTool*) this;
116 double returnRes;
117 double poidlobe;
118 double Puiss;
119 switch (Option)
120 {
121 case AllSeparable:
122 {
123 poidlobe=(pSigToolcur->pLobe)->weigthAmpl(VCur,VPointe,VY); // ss dimensions
124/*
125 if (poidlobe>.1)
126 { cout<<poidlobe<<endl;
127 }
128*/
129
130 Puiss=(pSigToolcur->pLSrc)->powerDensAmpli(VCur.Theta(),VCur.Phi());
131 // W m-2 st-1 Hz-1
132 returnRes=Puiss * poidlobe * IntegSpectOverFreq; // W / m2 / st
133 return returnRes;
134 }
135 case IsLightMapPowerInband:
136 {
137 // cout<<"VCur.Theta: "<<VCur.Theta()<<"VCur.Phi(): "<<VCur.Phi()<<endl;
138 poidlobe= (pSigToolcur->pLobe)->weigthAmpl(VCur,VPointe,VY);
139 Puiss= (pSigToolcur->pLSrc)->powerDensAmpli(VCur.Theta(),VCur.Phi());
140 returnRes=Puiss * poidlobe;
141 return returnRes;
142 }
143
144 default:
145 { // Cas NonSeparable
146 // Integration over at coordinates
147 returnRes=pIntegrale->Value();
148 return returnRes;
149 }
150
151 }
152}
153
154
155#define NBStepCircleMin (12)
156
157double SigCalcTool::powerInteg() {
158 // compute power on detector
159
160 double powerInteg=0.;
161 // Sum of the incominig power on detector.
162 UnitVector VPoin;
163 // VPointe Boresigth du telescope microonde
164 // VPoin direction priviliegiee du lobe, autour de laquelle on calcule
165 // VCur, vecteur courant du calcul.
166// double thetaCur, phiCur; // Coordinates of VCur
167 // Units is radian
168
169
170
171 //------Initialisation of Lobe integration------------------------------------------
172 double angShift=0.; // Angular distance from VPoin
173 double angShiftLimit; // On calcule jusqu'a angShiftLimit de VPoin
174
175 if(pLSrc->IsQPtSrc())
176 { double ang1=pLSrc->getAngSize()+pLobe->AngleMax();
177 VPoin=pLobe->VecShift(VPointe, VY);
178 if (ang1>=M_PI) { } //rien
179 else
180 { double cosinus=VPoin*pLSrc->GetVSrcCenter();
181 if (cosinus<cos(ang1)) return 0.;
182 //C'est le cas ou la source est trop loin de la direction pointŽe
183 }
184 // Maintenant on intgre
185 angShiftLimit=ang1;
186 }
187 else
188 {
189 VPoin=pLobe->VecShift(VPointe, VY);
190 angShiftLimit=pLobe->AngleMax();
191 }
192
193 // On va tourner autour de VPoin
194 // Compute unit vector perpendicular to Vpoin at same theta
195 UnitVector VPerp;
196 VPerp=VPoin.VperpPhi();
197
198 double dAngShift=AngResComp(0.)*1.1;
199 // AngleSteps are not necessarily constant.
200 // factor 1.1 to raise ambiguities in nearby pixel integration
201 double lastAngShiftMax;
202 // Needed to compute accurately solid angle values
203
204 VCur=VPoin;
205
206 powerInteg+=calcPowerDens()*diffSolidAng(0.,dAngShift/2.);
207 lastAngShiftMax= dAngShift/2.;
208
209 long NbPasOneCircle;
210 long CircleNumber=0; // no du cercle en cour:
211 // Gestion des dŽcalages pour un echantillonnage en quinconce
212 double solidAngStepCircle;
213 float stepAngCircle;
214
215 ///---------- Lobe integration-----------------------------------------
216 // generate vectors around VPoin at angular distance angShift.
217 // Compute power flux from foreground in this direction
218 // Weigth with weigth function and solid angle
219 dAngShift=AngResComp(lastAngShiftMax);
220
221 while((lastAngShiftMax+dAngShift)<angShiftLimit){
222 CircleNumber++;
223 angShift=lastAngShiftMax+dAngShift/2.;
224
225 VCur=VPoin.Rotate(VPerp,angShift);
226
227 // Compute number of step and associates on a circle
228 NbPasOneCircle=(long) (2*M_PI*sin(angShift)/sin(dAngShift));
229 if(NbPasOneCircle<NBStepCircleMin) NbPasOneCircle=NBStepCircleMin;
230 stepAngCircle=2*M_PI/NbPasOneCircle;
231 solidAngStepCircle= diffSolidAng(lastAngShiftMax,angShift+dAngShift/2.)/NbPasOneCircle;
232 // MRotAround=RotVec(VPoin,stepAngCircle);
233
234 //----------- integrate on a circle -------------------
235 if((CircleNumber%2)==0) VCur=VCur.Rotate(VPoin,stepAngCircle/2.);
236 // Pour un echantillonnage en quinconce
237
238 for(long i=0;i<NbPasOneCircle;i++)
239 {
240 //cout<< "rotation numb: "<< i<<endl;
241 powerInteg+=calcPowerDens()*solidAngStepCircle;
242 VCur=VCur.Rotate(VPoin,stepAngCircle);
243 } // end of circle
244
245 lastAngShiftMax+=dAngShift;
246 dAngShift=AngResComp(lastAngShiftMax);
247 }
248
249 // On s'occupe des effets de bord: un dernier tour!
250 // On change le code pour eviter les instabilites dues a dAngShift tres petit
251 CircleNumber++;
252 angShift=(angShiftLimit+lastAngShiftMax)/2.;
253
254 VCur=VPoin.Rotate(VPerp,angShift);
255 // Compute number of step and associates on a circle
256 NbPasOneCircle=(long) 2*M_PI*sin(angShift)/sin(AngResComp(angShift));
257 if(NbPasOneCircle<NBStepCircleMin) NbPasOneCircle=NBStepCircleMin;
258 stepAngCircle=2*M_PI/NbPasOneCircle;
259 solidAngStepCircle= diffSolidAng(lastAngShiftMax,angShiftLimit)/NbPasOneCircle;
260
261 //----------- integrate on last circle -------------------
262 for(long i=0;i<NbPasOneCircle;i++)
263 {
264 powerInteg+=calcPowerDens()*solidAngStepCircle;
265 VCur=VCur.Rotate(VPoin,stepAngCircle);
266 }
267 //end of last circle
268
269 //end of integration
270
271// cout<<"On a termine un point, OUFF"<< endl;
272 return powerInteg;
273}
274
275double SigCalcTool::NormKelvinRayleighJeans()
276{
277 double tempeCNoir=10000.;
278 // Kelvin
279 double CutFreq=1.380662e-23*tempeCNoir/6.626176e-34/5.;
280 if(FreqMax>1.380662e-23*tempeCNoir/6.626176e-34/5.)
281 { cerr<< "RaleighJeans approximation is not valid for this frequency"<<endl;
282 cerr<< "Frequency: "<< FreqMax<<" in SigCalcTool::NormRayleighJeans"<<endl;
283 }
284
285 LightBlackBody CorpsNoir(tempeCNoir, RAngComp);
286 SigCalcTool ToolRJ(&CorpsNoir,pLobe,pFilter);
287 double puissNorm = ToolRJ.compPixel(M_PI/2.,M_PI); // Un pixel au hasard
288 return tempeCNoir/puissNorm; // Kelvin RaleighJeans/(Watt/m2)
289
290}
291
292double SigCalcTool::NormKelvinCMB()
293{
294 double deltatempeCNoir=1.; // Kelvin
295 LightNormTCMB DeltaCorpsNoir(deltatempeCNoir, RAngComp);
296 SigCalcTool ToolDeltaCMB(&DeltaCorpsNoir,pLobe,pFilter);
297 double puissNorm = ToolDeltaCMB.compPixel(M_PI/2.,M_PI); // Un pixel au hasard
298 return deltatempeCNoir/puissNorm; // KelvinCMB/(Watt/m2)
299
300}
301
302/*
303double SigCalcTool::CalcInBandPower(double theta, double phi)
304{
305 double returnRes=0.;
306 UnitVector VP(theta,phi);
307 UnitVector VYbidon=VP.VperpPhi();
308// Compute unit vector perpendicular to Vpoin at same theta
309 VCur=VP;
310 VPointe=VP;
311 VY=VYbidon;
312 VX=VY^VP;
313 if(!emptySignal) returnRes=calcPowerDens(); // On integre sur la frequence
314 return returnRes;
315}
316*/
317
318double SigCalcTool::AngResComp(double angle) const
319{
320 double AngRes;
321 if(pLSrc->IsQPtSrc()) AngRes=RAngComp;
322 else AngRes=RAngComp*pLobe->ResolutionCurve(angle);
323 return AngRes;
324}
325
326
327
328double SigCalcTool::max(double a, double b) const{
329 if(a>b) return a;
330 else return b;
331}
332
333 double SigCalcTool::min(double a, double b) const{
334 if(a<b) return a;
335 else return b;
336}
337
338double SigCalcTool::CalcLobeSize(double frequency)
339{
340 // Compute lobe extension in steradians
341
342 if(frequency== -10.) frequency=(FreqMin+FreqMax)/2.;
343
344 double SizeInteg=0.;
345 // Sum of the incominig power on detector.
346 UnitVector VPoin;
347 // VPointe Boresigth du telescope microonde
348 // VPoin direction priviliegiee du lobe, autour de laquelle on calcule
349 // VCur, vecteur courant du calcul.
350
351 //------Initialisation of Lobe integration------------------------------------------
352 double angShift=0.; // Angular distance from VPoin
353 double angShiftLimit=pLobe->AngleMax(); // On calcule jusqu'a angShiftLimit de VPoin
354
355
356 // On va tourner autour de VPoin
357 // Compute unit vector perpendicular to Vpoin at same theta
358 UnitVector VPerp;
359 VPerp=VPoin.VperpPhi();
360
361 double dAngShift=AngResComp(0.)*1.1;
362 // AngleSteps are not necessarily constant.
363 // factor 1.1 to raise ambiguities in nearby pixel integration.
364 double lastAngShiftMax;
365 // Needed to compute accurately solid angle values
366 UnitVector VCur;
367 VCur=VPoin;
368
369 SizeInteg+= pLobe->weigth(VCur,VPoin,VPerp,frequency)*diffSolidAng(0.,dAngShift/2.);
370 lastAngShiftMax= dAngShift/2.;
371
372 long NbPasOneCircle;
373 long CircleNumber=0; // no du cercle en cour:
374 // Gestion des dŽcalages pour un echantillonnage en quinconce
375 double solidAngStepCircle;
376 float stepAngCircle;
377
378 ///---------- Lobe integration-----------------------------------------
379 // generate vectors around VPoin at angular distance angShift.
380 // Compute power flux from foreground in this direction
381 // Weigth with weigth function and solid angle
382 dAngShift=AngResComp(lastAngShiftMax);
383
384 while((lastAngShiftMax+dAngShift)<angShiftLimit)
385 {
386 CircleNumber++;
387 angShift=lastAngShiftMax+dAngShift/2.;
388
389 VCur=VPoin.Rotate(VPerp,angShift);
390
391 // Compute number of step and associates on a circle
392 NbPasOneCircle=(long) (2*M_PI*sin(angShift)/sin(dAngShift));
393 if(NbPasOneCircle<NBStepCircleMin) NbPasOneCircle=NBStepCircleMin;
394 stepAngCircle=2*M_PI/NbPasOneCircle;
395 solidAngStepCircle= diffSolidAng(lastAngShiftMax,angShift+dAngShift/2.)/NbPasOneCircle;
396
397 //----------- integrate on a circle -------------------
398 if((CircleNumber%2)==0) VCur=VCur.Rotate(VPoin,stepAngCircle/2.);
399 // Pour un echantillonnage en quinconce
400
401 for(long i=0;i<NbPasOneCircle;i++)
402 {
403 SizeInteg+= pLobe->weigth(VCur,VPoin,VPerp,frequency)*diffSolidAng(0.,dAngShift/2.);
404 VCur=VCur.Rotate(VPoin,stepAngCircle);
405 } // end of circle
406
407 lastAngShiftMax+=dAngShift;
408 dAngShift=AngResComp(lastAngShiftMax);
409 }
410
411 // On s'occupe des effets de bord: un dernier tour!
412 // On change le code pour eviter les instabilites dues a dAngShift tres petit
413 CircleNumber++;
414 angShift=(angShiftLimit+lastAngShiftMax)/2.;
415
416 VCur=VPoin.Rotate(VPerp,angShift);
417 // Compute number of step and associates on a circle
418 NbPasOneCircle=(long) 2*M_PI*sin(angShift)/sin(AngResComp(angShift));
419 if(NbPasOneCircle<NBStepCircleMin) NbPasOneCircle=NBStepCircleMin;
420 stepAngCircle=2*M_PI/NbPasOneCircle;
421 solidAngStepCircle= diffSolidAng(lastAngShiftMax,angShiftLimit)/NbPasOneCircle;
422
423 //----------- integrate on last circle -------------------
424 for(long i=0;i<NbPasOneCircle;i++)
425 {
426 SizeInteg+= pLobe->weigth(VCur,VPoin,VPerp,frequency)*diffSolidAng(0.,dAngShift/2.);
427 VCur=VCur.Rotate(VPoin,stepAngCircle);
428 }
429 //end of last circle
430
431 //end of integration
432
433 return SizeInteg;
434}
435
436double SigCalcTool::diffSolidAng(double ang1,double ang2) const
437{ double returnVal; // Steradians
438 if(!pLSrc->IsPtSourceS()) returnVal= fabs(2*M_PI*(cos(ang1)-cos(ang2)));
439 // Cas d'une source Žtendue.
440 else returnVal= 1.;
441 // Cas d'une source ponctuelle. Sont flux est exprime en W/m2/Hz
442 // Pas d'angle solide.
443 return returnVal;
444
445}
446
447// should be included as a class member, would template member function
448// work on all compilers
449
450static AbsLobeNoPolar* AddInBandPowerpLobe;
451static AbsLightSource* AddInBandPowerpLSrc;
452static SpectralResponse* AddInBandPowerpFilter;
453static double AIBtheta;
454static double AIBphi;
455
456static double AddInBandPowerFreqFunc1(double freq)
457{ // Integration function for GLInteg
458 double temp1= AddInBandPowerpLSrc->powSpecDens(AIBtheta,AIBphi,freq);
459 double temp2= AddInBandPowerpLobe->spectre(freq);
460 double temp3= AddInBandPowerpFilter->transmission(freq);
461 return temp1*temp2*temp3;
462}
463
464template <class T> void addInInBandPowerMap(PixelMap<T>& Map, SigCalcTool& Tool)
465{ // No spatial integration on the lobe
466 // Valid if lobe is separable in frequency
467 // Test
468 AddInBandPowerpLobe=Tool.getpLobe();
469 AddInBandPowerpLSrc=Tool.getpLSrc();
470 AddInBandPowerpFilter=Tool.getpFilter();
471 if(!AddInBandPowerpLobe->IsFreqSep())
472 { cerr<<" Adding power to a map using a lobe non separable in frequency is inconsistent"<<endl;
473 cerr<<" No power added, addInBandPower skipped"<<endl;
474 return;
475 }
476
477 long PixelNumber= Map.NbPixels();
478 double out;
479 T temp;
480 if(Tool.getOption()==AllSeparable)
481 { // Fast !
482 double FreqIntFactor=Tool.getIntegSpectOverFreq();
483 for(long k=0; k<PixelNumber; k++)
484 { Map.PixThetaPhi(k,AIBtheta,AIBphi);
485 out= AddInBandPowerpLSrc->powerDensAmpli(AIBtheta,AIBphi)*FreqIntFactor;
486 // Lobe weigth do no enters here
487 temp= (T) out;
488 Map(k)+= temp;
489 // if((k%200)==0) cout<<"200 points calculŽs "<<"NbPoint Total= "<<k<<endl;
490 }
491
492 }
493 else
494 {
495 if(AddInBandPowerpLSrc->IsFreqSep())
496 { double FreqMax=Tool.getFreqMax();
497 double FreqMin=Tool.getFreqMin();
498 double out;
499 GLInteg Integrale(AddInBandPowerFreqFunc1,FreqMin,FreqMax);
500 Integrale.NStep(10); // Serieux!
501 for(long k=0; k<PixelNumber; k++)
502 {
503 Map.PixThetaPhi(k,AIBtheta,AIBphi);
504 // Lobe weigth do no enters here
505 out=Integrale.Value();
506 // Lobe weigth do no enters here
507 temp= (T) out;
508 Map(k)+= temp;
509 }
510 }
511 }
512 return;
513}
514
515template void addInInBandPowerMap(PixelMap<float>& Map, SigCalcTool& tool);
516template void addInInBandPowerMap(PixelMap<double>& Map, SigCalcTool& tool);
Note: See TracBrowser for help on using the repository browser.