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

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

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