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

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

New Files
  • Property svn:executable set to *
File size: 9.1 KB
Line 
1
2
3 // Dominique YVON, CEA/DAPNIA/SPP 02/2000
4
5#include <math.h>
6#include <iostream>
7#include <fstream>
8#ifdef __MWERKS__
9
10 #include "unixmac.h"
11#endif
12
13//#include "integ.h"
14
15#include "abscalctool.h"
16
17#include "lightdipole.h"
18#include "sigcalctools.h"
19
20
21double AbsCalcTool::compPixelQD(double theta,double phi)
22{
23// cout<<"theta: "<<theta<<"\tphi: "<<phi<<endl;
24 UnitVector VP(theta,phi);
25 UnitVector VYbidon=VP.VperpPhi();
26// Compute unit vector perpendicular to Vpoin at same theta
27 return compPixel(VP,VYbidon);
28}
29
30double AbsCalcTool::NormKelvinRayleighJeans()
31{
32 double tempeCNoir=10000.;
33 // Kelvin
34 double CutFreq=1.380662e-23*tempeCNoir/6.626176e-34/5.;
35 if(FreqMax>1.380662e-23*tempeCNoir/6.626176e-34/5.)
36 { cerr<< "RaleighJeans approximation is not valid for this frequency"<<endl;
37 cerr<< "Frequency: "<< FreqMax<<" in SigCalcTool::NormRayleighJeans"<<endl;
38 }
39
40 LightBlackBody CorpsNoir(tempeCNoir, RAngComp);
41 SigCalcTool ToolRJ(&CorpsNoir,pLobe,pFilter);
42 double puissNorm = ToolRJ.compPixelQD(M_PI/2.,M_PI); // Un pixel au hasard
43 return tempeCNoir/puissNorm; // Kelvin RaleighJeans/(Watt/m2)
44
45}
46
47double AbsCalcTool::NormKelvinCMB()
48{
49 double deltatempeCNoir=1.; // Kelvin
50 LightNormTCMB DeltaCorpsNoir(deltatempeCNoir, RAngComp);
51 SigCalcTool ToolDeltaCMB(&DeltaCorpsNoir,pLobe,pFilter);
52 double puissNorm = ToolDeltaCMB.compPixelQD(M_PI/2.,M_PI); // Un pixel au hasard
53 return deltatempeCNoir/puissNorm; // KelvinCMB/(Watt/m2)
54}
55
56double AbsCalcTool::diffSolidAng(double ang1,double ang2) const
57{ return fabs(2*M_PI*(cos(ang1)-cos(ang2))); // Steradians
58 // Cas d'une source Žtendue.
59}
60
61double AbsCalcTool::max(double a, double b) const{
62 if(a>b) return a;
63 else return b;
64}
65
66double AbsCalcTool::min(double a, double b) const{
67 if(a<b) return a;
68 else return b;
69}
70
71int kmg_eulerRad(double ai, double bi, int select, double *ao, double *bo) {
72 // All coordinates are in Radian.
73 /* kmg_euler.c
74 *
75 * Converts between different coordinate systems.
76 *
77 * SELECT From To | SELECT From To
78 * 1 RA-Dec (2000) Galactic | 4 Ecliptic RA-Dec
79 * 2 Galactic RA-DEC | 5 Ecliptic
80Galactic
81 * 3 RA-Dec Ecliptic | 6 Galactic
82Ecliptic
83 *
84 * Date Programmer Remarks
85 * ----------- ---------- -------
86 * 08-Aug-1999 K. Ganga First version. Copied and modified EULER from
87 * the IDL Astrolib.
88 * May 2000, D. Yvon Change coordinates units to Radians
89 */
90
91 /* Local Declarations */
92 double a, b, sb, cb, cbsa;
93 long i;
94 const double TWOPI = 2.0*M_PI;
95 const double FOURPI = 4.0*M_PI;
96 const double DEG2RAD = 180.0/M_PI;
97
98 /* J2000 coordinate conversions are based on the following constants
99 * eps = 23.4392911111 Obliquity of the ecliptic
100 * alphaG = 192.85948 Right Ascension of Galactic North Pole
101 * deltaG = 27.12825 Declination of Galactic North Pole
102 * lomega = 32.93192 Galactic longitude of celestial equator
103 * alphaE = 180.02322 Ecliptic longitude of Galactic North Pole
104 * deltaE = 29.811438523 Ecliptic latitude of Galactic North Pole
105 * Eomega = 6.3839743 Galactic longitude of ecliptic equator
106 */
107
108 const double psi[6] = {0.57477043300, 4.9368292465 ,
109 0.00000000000, 0.0000000000 ,
110 0.11142137093, 4.71279419371};
111 const double stheta[6] = {0.88998808748,-0.88998808748,
112 0.39777715593,-0.39777715593,
113 0.86766622025,-0.86766622025};
114 const double ctheta[6] = {0.45598377618, 0.45598377618,
115 0.91748206207, 0.91748206207,
116 0.49714719172, 0.49714719172};
117 const double phi[6] = {4.9368292465 , 0.57477043300,
118 0.0000000000 , 0.00000000000,
119 4.71279419371, 0.11142137093};
120
121 i = select - 1;
122 a = ai - phi[i];
123 b = bi;
124 sb = sin(b);
125 cb = cos(b);
126 cbsa = cb*sin(a);
127 b = -stheta[i]*cbsa + ctheta[i]*sb;
128 b = ( (b > 1.0) ? 1.0 : b );
129 b = ( (b < -1.0) ? -1.0 : b );
130 *bo = asin(b);
131
132 a = atan2( ctheta[i] * cbsa + stheta[i] * sb, cb * cos(a) );
133 *ao = fmod(a + psi[i] + FOURPI, TWOPI);
134
135 /* Later */
136 return 0;
137}
138
139template <class T> void addToSkyMap(PixelMap<T>& Map, AbsCalcTool& Tool)
140{
141 double theta, phi;
142 long PixelNumber= Map.NbPixels();
143 cout<<"Nbre de pIxel a calculer dans addToSkyMap: "<<PixelNumber<<endl;
144 T temp;
145 for(long k=0; k<PixelNumber; k++)
146 {
147 Map.PixThetaPhi(k,theta,phi);
148 temp=(T) Tool.compPixelQD(theta,phi);
149// if (temp<1.e-35) temp=0.; // BUGG XXXXXX
150 // pour vivre avec les pb de fichiers fits D.Y.
151 Map(k)+=temp;
152 if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
153 }
154 return;
155}
156
157template <class T> void compSkyMap(PixelMap<T>& Map, AbsCalcTool& Tool)
158{
159 double theta, phi;
160 long PixelNumber= Map.NbPixels();
161 cout<<"Nbre de pIxel a calculer dans compSkyMap: "<<PixelNumber<<endl;
162 T temp;
163 for(long k=0; k<PixelNumber; k++)
164 {
165 Map.PixThetaPhi(k,theta,phi);
166 temp= (T) Tool.compPixelQD(theta,phi);
167// if (temp<1.e-35) temp=0.;
168 // BUGG XXXXXX
169 // pour vivre avec les pb de fichiers fits D.Y.
170 //if (out!=0.) cout<< out <<'\t'<< temp<<endl;
171 Map(k)= temp;
172 if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
173 }
174 return;
175}
176
177template <class T1, class T2> void addMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
178{
179 double theta, phi;
180 long PixelNumber= Map.NbPixels();
181 //T1 temp;
182 for(long k=0; k<PixelNumber; k++)
183 {
184 Map.PixThetaPhi(k,theta,phi);
185 // temp=(T1) Map2.PixValSph(theta,phi);
186 if(Map2.ContainsSph(theta,phi)) Map(k) += (T1) Map2.PixValSph(theta,phi);
187 }
188 return;
189}
190
191template <class T1, class T2> void substractMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
192{
193 double theta, phi;
194 long PixelNumber= Map.NbPixels();
195 //T1 temp;
196 for(long k=0; k<PixelNumber; k++)
197 {
198 Map.PixThetaPhi(k,theta,phi);
199 // temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
200 if(Map2.ContainsSph(theta,phi)) Map(k) = Map(k)-(T1) Map2.PixValSph(theta,phi);
201 }
202 return;
203}
204
205template <class T1, class T2> void divMap1WithMap2(PixelMap<T1>& Map, PixelMap<T2>& Map2)
206{
207 double theta, phi;
208 long PixelNumber= Map.NbPixels();
209 //T1 temp;
210 for(long k=0; k<PixelNumber; k++)
211 {
212 Map.PixThetaPhi(k,theta,phi);
213 // temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
214 if(Map2.ContainsSph(theta,phi))
215 if (Map2.PixValSph(theta,phi)!=0.)
216 Map(k) = Map(k)/(T1) Map2.PixValSph(theta,phi);
217 }
218 return;
219}
220
221template <class T> void scaleMap(double scalefactor, PixelMap<T>& Map)
222{
223 long PixelNumber= Map.NbPixels();
224 //T temp;
225 for(long k=0; k<PixelNumber; k++)
226 {
227// temp= (T) ( scalefactor*Map2(k) );
228 Map(k) = (T) ( scalefactor*Map(k) );
229 }
230}
231
232
233template <class T> int MinMaxSigMap(PixelMap<T>& Map, double& Min,
234 double& Max, double& Moy, double& sigma)
235{
236 long PixelNumber= Map.NbPixels();
237 double val=0.;
238 double variance=0;
239 Moy=0.;
240 Min=1.e36;
241 Max=-1.e-36;
242 sigma=0.;
243
244 for(long k=0; k<PixelNumber; k++)
245 { val=Map(k); // cout<<k<<'\t'<<val<<'\n';
246 if(val>Max) Max=val;
247 if(val<Min) Min=val;
248 Moy+=val;
249 }
250
251 Moy/=PixelNumber;
252 for(long k=0; k<PixelNumber; k++) variance+=(val-Moy)*(val-Moy);
253 variance/=PixelNumber;
254 sigma=sqrt(variance);
255
256 if(Min==Max)
257 { cerr<< "bizarre, une carte vide ou uniforme"<<endl;
258 return 1;
259 }
260 return 0;
261}
262
263// On instancie les plus courants
264template void addToSkyMap(PixelMap<double>& Map, AbsCalcTool& Tool);
265template void addToSkyMap(PixelMap<float>& Map, AbsCalcTool& Tool);
266template void compSkyMap(PixelMap<double>& Map, AbsCalcTool& Tool);
267template void compSkyMap(PixelMap<float>& Map, AbsCalcTool& Tool);
268
269// Idem pour les outils de cartes
270
271template void scaleMap(double scalefactor, PixelMap<double>& Map2);
272template void scaleMap(double scalefactor, PixelMap<float>& Map2);
273
274template int MinMaxSigMap(PixelMap<float>& Map, double& Min,
275 double& Max, double& Moy, double& sigma);
276template int MinMaxSigMap(PixelMap<double>& Map, double& Min,
277 double& Max, double& Moy, double& sigma);
278
279
280template void substractMap(PixelMap<float>& Map, PixelMap<float>& Map2);
281template void addMap(PixelMap<float>& Map, PixelMap<float>& Map2);
282template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<float>& Map2);
283
284template void substractMap(PixelMap<float>& Map, PixelMap<double>& Map2);
285template void addMap(PixelMap<float>& Map, PixelMap<double>& Map2);
286template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<double>& Map2);
287
288template void substractMap(PixelMap<double>& Map, PixelMap<float>& Map2);
289template void addMap(PixelMap<double>& Map, PixelMap<float>& Map2);
290template void divMap1WithMap2(PixelMap<double>& Map, PixelMap<float>& Map2);
291
292template void substractMap(PixelMap<double>& Map, PixelMap<double>& Map2);
293template void addMap(PixelMap<double>& Map, PixelMap<double>& Map2);
294template void divMap1WithMap2(PixelMap<double>& Map, PixelMap<double>& Map2);
295
Note: See TracBrowser for help on using the repository browser.