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

Last change on this file since 1184 was 1154, checked in by ansari, 25 years ago

portage Cxx
D.Y.
  • Property svn:executable set to *
File size: 9.1 KB
RevLine 
[1149]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
[1154]61/* // Pour compilation sur cxx boeuf!
[1149]62double AbsCalcTool::max(double a, double b) const{
63 if(a>b) return a;
64 else return b;
65}
66
67double AbsCalcTool::min(double a, double b) const{
68 if(a<b) return a;
69 else return b;
70}
[1154]71*/
[1149]72int kmg_eulerRad(double ai, double bi, int select, double *ao, double *bo) {
73 // All coordinates are in Radian.
74 /* kmg_euler.c
75 *
76 * Converts between different coordinate systems.
77 *
78 * SELECT From To | SELECT From To
79 * 1 RA-Dec (2000) Galactic | 4 Ecliptic RA-Dec
80 * 2 Galactic RA-DEC | 5 Ecliptic
81Galactic
82 * 3 RA-Dec Ecliptic | 6 Galactic
83Ecliptic
84 *
85 * Date Programmer Remarks
86 * ----------- ---------- -------
87 * 08-Aug-1999 K. Ganga First version. Copied and modified EULER from
88 * the IDL Astrolib.
89 * May 2000, D. Yvon Change coordinates units to Radians
90 */
91
92 /* Local Declarations */
93 double a, b, sb, cb, cbsa;
94 long i;
95 const double TWOPI = 2.0*M_PI;
96 const double FOURPI = 4.0*M_PI;
97 const double DEG2RAD = 180.0/M_PI;
98
99 /* J2000 coordinate conversions are based on the following constants
100 * eps = 23.4392911111 Obliquity of the ecliptic
101 * alphaG = 192.85948 Right Ascension of Galactic North Pole
102 * deltaG = 27.12825 Declination of Galactic North Pole
103 * lomega = 32.93192 Galactic longitude of celestial equator
104 * alphaE = 180.02322 Ecliptic longitude of Galactic North Pole
105 * deltaE = 29.811438523 Ecliptic latitude of Galactic North Pole
106 * Eomega = 6.3839743 Galactic longitude of ecliptic equator
107 */
108
109 const double psi[6] = {0.57477043300, 4.9368292465 ,
110 0.00000000000, 0.0000000000 ,
111 0.11142137093, 4.71279419371};
112 const double stheta[6] = {0.88998808748,-0.88998808748,
113 0.39777715593,-0.39777715593,
114 0.86766622025,-0.86766622025};
115 const double ctheta[6] = {0.45598377618, 0.45598377618,
116 0.91748206207, 0.91748206207,
117 0.49714719172, 0.49714719172};
118 const double phi[6] = {4.9368292465 , 0.57477043300,
119 0.0000000000 , 0.00000000000,
120 4.71279419371, 0.11142137093};
121
122 i = select - 1;
123 a = ai - phi[i];
124 b = bi;
125 sb = sin(b);
126 cb = cos(b);
127 cbsa = cb*sin(a);
128 b = -stheta[i]*cbsa + ctheta[i]*sb;
129 b = ( (b > 1.0) ? 1.0 : b );
130 b = ( (b < -1.0) ? -1.0 : b );
131 *bo = asin(b);
132
133 a = atan2( ctheta[i] * cbsa + stheta[i] * sb, cb * cos(a) );
134 *ao = fmod(a + psi[i] + FOURPI, TWOPI);
135
136 /* Later */
137 return 0;
138}
139
140template <class T> void addToSkyMap(PixelMap<T>& Map, AbsCalcTool& Tool)
141{
142 double theta, phi;
143 long PixelNumber= Map.NbPixels();
144 cout<<"Nbre de pIxel a calculer dans addToSkyMap: "<<PixelNumber<<endl;
145 T temp;
146 for(long k=0; k<PixelNumber; k++)
147 {
148 Map.PixThetaPhi(k,theta,phi);
149 temp=(T) Tool.compPixelQD(theta,phi);
150// if (temp<1.e-35) temp=0.; // BUGG XXXXXX
151 // pour vivre avec les pb de fichiers fits D.Y.
152 Map(k)+=temp;
153 if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
154 }
155 return;
156}
157
158template <class T> void compSkyMap(PixelMap<T>& Map, AbsCalcTool& Tool)
159{
160 double theta, phi;
161 long PixelNumber= Map.NbPixels();
162 cout<<"Nbre de pIxel a calculer dans compSkyMap: "<<PixelNumber<<endl;
163 T temp;
164 for(long k=0; k<PixelNumber; k++)
165 {
166 Map.PixThetaPhi(k,theta,phi);
167 temp= (T) Tool.compPixelQD(theta,phi);
168// if (temp<1.e-35) temp=0.;
169 // BUGG XXXXXX
170 // pour vivre avec les pb de fichiers fits D.Y.
171 //if (out!=0.) cout<< out <<'\t'<< temp<<endl;
172 Map(k)= temp;
173 if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
174 }
175 return;
176}
177
178template <class T1, class T2> void addMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
179{
180 double theta, phi;
181 long PixelNumber= Map.NbPixels();
182 //T1 temp;
183 for(long k=0; k<PixelNumber; k++)
184 {
185 Map.PixThetaPhi(k,theta,phi);
186 // temp=(T1) Map2.PixValSph(theta,phi);
187 if(Map2.ContainsSph(theta,phi)) Map(k) += (T1) Map2.PixValSph(theta,phi);
188 }
189 return;
190}
191
192template <class T1, class T2> void substractMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
193{
194 double theta, phi;
195 long PixelNumber= Map.NbPixels();
196 //T1 temp;
197 for(long k=0; k<PixelNumber; k++)
198 {
199 Map.PixThetaPhi(k,theta,phi);
200 // temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
201 if(Map2.ContainsSph(theta,phi)) Map(k) = Map(k)-(T1) Map2.PixValSph(theta,phi);
202 }
203 return;
204}
205
206template <class T1, class T2> void divMap1WithMap2(PixelMap<T1>& Map, PixelMap<T2>& Map2)
207{
208 double theta, phi;
209 long PixelNumber= Map.NbPixels();
210 //T1 temp;
211 for(long k=0; k<PixelNumber; k++)
212 {
213 Map.PixThetaPhi(k,theta,phi);
214 // temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
215 if(Map2.ContainsSph(theta,phi))
216 if (Map2.PixValSph(theta,phi)!=0.)
217 Map(k) = Map(k)/(T1) Map2.PixValSph(theta,phi);
218 }
219 return;
220}
221
222template <class T> void scaleMap(double scalefactor, PixelMap<T>& Map)
223{
224 long PixelNumber= Map.NbPixels();
225 //T temp;
226 for(long k=0; k<PixelNumber; k++)
227 {
228// temp= (T) ( scalefactor*Map2(k) );
229 Map(k) = (T) ( scalefactor*Map(k) );
230 }
231}
232
233
234template <class T> int MinMaxSigMap(PixelMap<T>& Map, double& Min,
235 double& Max, double& Moy, double& sigma)
236{
237 long PixelNumber= Map.NbPixels();
238 double val=0.;
239 double variance=0;
240 Moy=0.;
241 Min=1.e36;
242 Max=-1.e-36;
243 sigma=0.;
244
245 for(long k=0; k<PixelNumber; k++)
246 { val=Map(k); // cout<<k<<'\t'<<val<<'\n';
247 if(val>Max) Max=val;
248 if(val<Min) Min=val;
249 Moy+=val;
250 }
251
252 Moy/=PixelNumber;
253 for(long k=0; k<PixelNumber; k++) variance+=(val-Moy)*(val-Moy);
254 variance/=PixelNumber;
255 sigma=sqrt(variance);
256
257 if(Min==Max)
258 { cerr<< "bizarre, une carte vide ou uniforme"<<endl;
259 return 1;
260 }
261 return 0;
262}
263
264// On instancie les plus courants
265template void addToSkyMap(PixelMap<double>& Map, AbsCalcTool& Tool);
266template void addToSkyMap(PixelMap<float>& Map, AbsCalcTool& Tool);
267template void compSkyMap(PixelMap<double>& Map, AbsCalcTool& Tool);
268template void compSkyMap(PixelMap<float>& Map, AbsCalcTool& Tool);
269
270// Idem pour les outils de cartes
271
272template void scaleMap(double scalefactor, PixelMap<double>& Map2);
273template void scaleMap(double scalefactor, PixelMap<float>& Map2);
274
275template int MinMaxSigMap(PixelMap<float>& Map, double& Min,
276 double& Max, double& Moy, double& sigma);
277template int MinMaxSigMap(PixelMap<double>& Map, double& Min,
278 double& Max, double& Moy, double& sigma);
279
280
281template void substractMap(PixelMap<float>& Map, PixelMap<float>& Map2);
282template void addMap(PixelMap<float>& Map, PixelMap<float>& Map2);
283template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<float>& Map2);
284
285template void substractMap(PixelMap<float>& Map, PixelMap<double>& Map2);
286template void addMap(PixelMap<float>& Map, PixelMap<double>& Map2);
287template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<double>& Map2);
288
289template void substractMap(PixelMap<double>& Map, PixelMap<float>& Map2);
290template void addMap(PixelMap<double>& Map, PixelMap<float>& Map2);
291template void divMap1WithMap2(PixelMap<double>& Map, PixelMap<float>& Map2);
292
293template void substractMap(PixelMap<double>& Map, PixelMap<double>& Map2);
294template void addMap(PixelMap<double>& Map, PixelMap<double>& Map2);
295template void divMap1WithMap2(PixelMap<double>& Map, PixelMap<double>& Map2);
296
Note: See TracBrowser for help on using the repository browser.