source: Sophya/trunk/SigPredictor/sigcalchackingtools.cc@ 2222

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

mise a jour
File size: 7.5 KB
Line 
1 // Dominique YVON, CEA/DAPNIA/SPP 02/2000
2
3#include <math.h>
4#include <iostream>
5#include <fstream>
6#ifdef __MWERKS__
7
8 #include "unixmac.h"
9#endif
10
11#include "sigcalctools.h"
12//#include "integ.h"
13
14 int kmg_eulerRad(double ai, double bi, int select, double *ao, double *bo) {
15 // All coordinates are in Radian.
16 /* kmg_euler.c
17 *
18 * Converts between different coordinate systems.
19 *
20 * SELECT From To | SELECT From To
21 * 1 RA-Dec (2000) Galactic | 4 Ecliptic RA-Dec
22 * 2 Galactic RA-DEC | 5 Ecliptic
23Galactic
24 * 3 RA-Dec Ecliptic | 6 Galactic
25Ecliptic
26 *
27 * Date Programmer Remarks
28 * ----------- ---------- -------
29 * 08-Aug-1999 K. Ganga First version. Copied and modified EULER from
30 * the IDL Astrolib.
31 * May 2000, D. Yvon Change coordinates units to Radians
32 */
33
34 /* Local Declarations */
35 double a, b, sb, cb, cbsa;
36 long i;
37 const double TWOPI = 2.0*M_PI;
38 const double FOURPI = 4.0*M_PI;
39 const double DEG2RAD = 180.0/M_PI;
40
41 /* J2000 coordinate conversions are based on the following constants
42 * eps = 23.4392911111 Obliquity of the ecliptic
43 * alphaG = 192.85948 Right Ascension of Galactic North Pole
44 * deltaG = 27.12825 Declination of Galactic North Pole
45 * lomega = 32.93192 Galactic longitude of celestial equator
46 * alphaE = 180.02322 Ecliptic longitude of Galactic North Pole
47 * deltaE = 29.811438523 Ecliptic latitude of Galactic North Pole
48 * Eomega = 6.3839743 Galactic longitude of ecliptic equator
49 */
50
51 const double psi[6] = {0.57477043300, 4.9368292465 ,
52 0.00000000000, 0.0000000000 ,
53 0.11142137093, 4.71279419371};
54 const double stheta[6] = {0.88998808748,-0.88998808748,
55 0.39777715593,-0.39777715593,
56 0.86766622025,-0.86766622025};
57 const double ctheta[6] = {0.45598377618, 0.45598377618,
58 0.91748206207, 0.91748206207,
59 0.49714719172, 0.49714719172};
60 const double phi[6] = {4.9368292465 , 0.57477043300,
61 0.0000000000 , 0.00000000000,
62 4.71279419371, 0.11142137093};
63
64 i = select - 1;
65 a = ai - phi[i];
66 b = bi;
67 sb = sin(b);
68 cb = cos(b);
69 cbsa = cb*sin(a);
70 b = -stheta[i]*cbsa + ctheta[i]*sb;
71 b = ( (b > 1.0) ? 1.0 : b );
72 b = ( (b < -1.0) ? -1.0 : b );
73 *bo = asin(b);
74
75 a = atan2( ctheta[i] * cbsa + stheta[i] * sb, cb * cos(a) );
76 *ao = fmod(a + psi[i] + FOURPI, TWOPI);
77
78 /* Later */
79 return 0;
80}
81
82
83template <class T> void addToSkyMap(PixelMap<T>& Map, SigCalcTool& Tool)
84{
85 double theta, phi;
86 long PixelNumber= Map.NbPixels();
87 cout<<"Nbre de pIxel a calculer dans addToSkyMap: "<<PixelNumber<<endl;
88 T temp;
89 for(long k=0; k<PixelNumber; k++)
90 {
91 Map.PixThetaPhi(k,theta,phi);
92 temp=(T) Tool.compPixel(theta,phi);
93// if (temp<1.e-35) temp=0.; // BUGG XXXXXX
94 // pour vivre avec les pb de fichiers fits D.Y.
95 Map(k)+=temp;
96 if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
97 }
98 return;
99}
100
101template <class T> void compSkyMap(PixelMap<T>& Map, SigCalcTool& Tool)
102{
103 double theta, phi;
104 long PixelNumber= Map.NbPixels();
105 cout<<"Nbre de pIxel a calculer dans compSkyMap: "<<PixelNumber<<endl;
106 T temp;
107 for(long k=0; k<PixelNumber; k++)
108 {
109 Map.PixThetaPhi(k,theta,phi);
110 temp= (T) Tool.compPixel(theta,phi);
111// if (temp<1.e-35) temp=0.;
112 // BUGG XXXXXX
113 // pour vivre avec les pb de fichiers fits D.Y.
114 //if (out!=0.) cout<< out <<'\t'<< temp<<endl;
115 Map(k)= temp;
116 if((k%500)==0) cout<<"500 points calculŽs "<<"NbPoint Total= "<<k<<endl;
117 }
118 return;
119}
120
121template <class T1, class T2> void addMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
122{
123 double theta, phi;
124 long PixelNumber= Map.NbPixels();
125 //T1 temp;
126 for(long k=0; k<PixelNumber; k++)
127 {
128 Map.PixThetaPhi(k,theta,phi);
129 // temp=(T1) Map2.PixValSph(theta,phi);
130 if(Map2.ContainsSph(theta,phi)) Map(k) += (T1) Map2.PixValSph(theta,phi);
131 }
132 return;
133}
134
135template <class T1, class T2> void substractMap(PixelMap<T1>& Map, PixelMap<T2>& Map2)
136{
137 double theta, phi;
138 long PixelNumber= Map.NbPixels();
139 //T1 temp;
140 for(long k=0; k<PixelNumber; k++)
141 {
142 Map.PixThetaPhi(k,theta,phi);
143 // temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
144 if(Map2.ContainsSph(theta,phi)) Map(k) = Map(k)-(T1) Map2.PixValSph(theta,phi);
145 }
146 return;
147}
148
149template <class T1, class T2> void divMap1WithMap2(PixelMap<T1>& Map, PixelMap<T2>& Map2)
150{
151 double theta, phi;
152 long PixelNumber= Map.NbPixels();
153 //T1 temp;
154 for(long k=0; k<PixelNumber; k++)
155 {
156 Map.PixThetaPhi(k,theta,phi);
157 // temp=(T1) Map2.PixValSph(theta,phi); cout<<" Map2.PixValSph(theta,phi): "<<temp<<endl;
158 if(Map2.ContainsSph(theta,phi))
159 if (Map2.PixValSph(theta,phi)!=0.)
160 Map(k) = Map(k)/(T1) Map2.PixValSph(theta,phi);
161 }
162 return;
163}
164
165template <class T> void scaleMap(double scalefactor, PixelMap<T>& Map)
166{
167 long PixelNumber= Map.NbPixels();
168 //T temp;
169 for(long k=0; k<PixelNumber; k++)
170 {
171// temp= (T) ( scalefactor*Map2(k) );
172 Map(k) = (T) ( scalefactor*Map(k) );
173 }
174}
175
176
177template <class T> int MinMaxSigMap(PixelMap<T>& Map, double& Min,
178 double& Max, double& Moy, double& sigma)
179{
180 long PixelNumber= Map.NbPixels();
181 double val=0.;
182 double variance=0;
183 Moy=0.;
184 Min=1.e36;
185 Max=-1.e-36;
186 sigma=0.;
187
188 for(long k=0; k<PixelNumber; k++)
189 { val=Map(k); // cout<<k<<'\t'<<val<<'\n';
190 if(val>Max) Max=val;
191 if(val<Min) Min=val;
192 Moy+=val;
193 }
194
195 Moy/=PixelNumber;
196 for(long k=0; k<PixelNumber; k++) variance+=(val-Moy)*(val-Moy);
197 variance/=PixelNumber;
198 sigma=sqrt(variance);
199
200 if(Min==Max)
201 { cerr<< "bizarre, une carte vide ou uniforme"<<endl;
202 return 1;
203 }
204 return 0;
205}
206
207// On instancie les plus courants
208template void addToSkyMap(PixelMap<double>& Map, SigCalcTool& Tool);
209template void addToSkyMap(PixelMap<float>& Map, SigCalcTool& Tool);
210template void compSkyMap(PixelMap<double>& Map, SigCalcTool& Tool);
211template void compSkyMap(PixelMap<float>& Map, SigCalcTool& Tool);
212
213// Idem pour les outils de cartes
214
215template void scaleMap(double scalefactor, PixelMap<double>& Map2);
216template void scaleMap(double scalefactor, PixelMap<float>& Map2);
217
218template int MinMaxSigMap(PixelMap<float>& Map, double& Min,
219 double& Max, double& Moy, double& sigma);
220template int MinMaxSigMap(PixelMap<double>& Map, double& Min,
221 double& Max, double& Moy, double& sigma);
222
223
224template void substractMap(PixelMap<float>& Map, PixelMap<float>& Map2);
225template void addMap(PixelMap<float>& Map, PixelMap<float>& Map2);
226template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<float>& Map2);
227
228template void substractMap(PixelMap<float>& Map, PixelMap<double>& Map2);
229template void addMap(PixelMap<float>& Map, PixelMap<double>& Map2);
230template void divMap1WithMap2(PixelMap<float>& Map, PixelMap<double>& Map2);
231
232template void substractMap(PixelMap<double>& Map, PixelMap<float>& Map2);
233template void addMap(PixelMap<double>& Map, PixelMap<float>& Map2);
234template void divMap1WithMap2(PixelMap<double>& Map, PixelMap<float>& Map2);
235
236template void substractMap(PixelMap<double>& Map, PixelMap<double>& Map2);
237template void addMap(PixelMap<double>& Map, PixelMap<double>& Map2);
238template void divMap1WithMap2(PixelMap<double>& Map, PixelMap<double>& Map2);
239
Note: See TracBrowser for help on using the repository browser.