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

Last change on this file since 3935 was 1191, checked in by ansari, 25 years ago

cleaned up for maching LevelS upgrade ongoing

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