source: Sophya/trunk/Cosmo/RadioBeam/mdish.cc@ 3789

Last change on this file since 3789 was 3789, checked in by ansari, 15 years ago

Corrections diverses, Reza 27/06/2010
File size: 8.9 KB
Line 
1// Classes to compute 2D
2// R. Ansari - Nov 2008, May 2010
3
4#include "mdish.h"
5
6
7//--------------------------------------------------
8// -- Four2DResponse class
9//--------------------------------------------------
10// Constructor
11Four2DResponse::Four2DResponse(int typ, double dx, double dy, double lambda)
12 : typ_(typ), dx_((dx>1.e-3)?dx:1.), dy_((dy>1.e-3)?dy:1.)
13{
14 setLambdaRef(lambda);
15 setLambda(lambda);
16}
17
18// Return the response for the wave vecteor (kx,ky)
19double Four2DResponse::Value(double kx, double ky)
20{
21 kx *= lambda_ratio_;
22 ky *= lambda_ratio_;
23 double wk,wkx,wky;
24 switch (typ_)
25 {
26 case 1: // Reponse gaussienne parabole diametre D exp[ - 0.5 (lambda k_g / D )^2 ]
27 wk = sqrt(kx*kx+ky*ky)/dx_;
28 wk = 0.5*wk*wk;
29 return exp(-wk);
30 break;
31 case 2: // Reponse parabole diametre D Triangle <= kmax= 2 pi D / lambda
32 wk = sqrt(kx*kx+ky*ky)/dx_/2./M_PI;
33 return ( (wk<1.)?(1.-wk):0.);
34 break;
35 case 3: // Reponse rectangle Dx x Dy Triangle (|kx|,|k_y|) <= (2 pi Dx / lambda, 2 pi Dx / lambda)
36 wkx = kx/2./M_PI/dx_;
37 wky = ky/2./M_PI/dy_;
38 return ( ((wkx<1.)&&(wky<1.))?((1.-wkx)*(1-wky)):0.);
39 break;
40 default:
41 return 1.;
42 }
43}
44// Return a vector representing the power spectrum (for checking)
45Histo2D Four2DResponse::GetResponse(int nx, int ny)
46{
47 double kxmx = 1.2*DeuxPI*dx_;
48 double kymx = 1.2*DeuxPI*dy_;
49 if (typ_<3) kymx=kxmx;
50 Histo2D h2(0.,kxmx,nx,0.,kymx,ny);
51
52 for(int j=0; j<h2.NBinY(); j++)
53 for(int i=0; i<h2.NBinX(); i++)
54 h2(i,j) = Value((i+0.5)*h2.WBinX(), (j+0.5)*h2.WBinY());
55 return h2;
56}
57
58//---------------------------------------------------------------
59// -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi)
60//---------------------------------------------------------------
61Four2DRespTable::Four2DRespTable(Histo2D const & hrep, double d)
62 : Four2DResponse(0,d,d) , hrep_(hrep)
63{
64}
65
66double Four2DRespTable::Value(double kx, double ky)
67{
68 kx *= lambda_ratio_;
69 ky *= lambda_ratio_;
70 int_4 i,j;
71 if ( (kx<=hrep_.XMin())||(kx>=hrep_.XMax()) ||
72 (ky<=hrep_.YMin())||(ky>=hrep_.YMax()) ) return 0.;
73 hrep_.FindBin(kx, ky, i, j);
74 return hrep_(i, j);
75}
76
77//---------------------------------------------------------------
78// -- Four2DRespRatio : rapport de la reponse entre deux objets Four2DResponse
79//---------------------------------------------------------------
80Four2DRespRatio::Four2DRespRatio(Four2DResponse& a, Four2DResponse& b, double divzthr)
81 : Four2DResponse(0, a.D(), a.D()), a_(a), b_(b), divzthr_(divzthr)
82{
83}
84
85double Four2DRespRatio::Value(double kx, double ky)
86{
87 double ra = a_.Value(kx,ky);
88 double rb = b_.Value(kx,ky);
89 if (rb<divzthr_) {
90 if (ra<rb) return 0.;
91 else rb=divzthr_;
92 }
93 return (ra/rb);
94}
95
96//---------------------------------------------------------------
97//--- Classe simple pour le calcul de rotation
98class Rotation {
99public:
100 Rotation(double tet, double phi)
101 {
102// (Teta,Phi) = Direction de visee
103// Les angles d'Euler correspondants sont Teta, Phi+Pi/2
104// Le Pi/2 vient que les rotations d'euler se font dans l'ordre
105// Autour de oZ d'angle Phi, autour de oN (nouvel axe X) d'angle Teta
106// Autour du nouvel axe Z (x3) d'angle Psi (Psi=0 -> cp=1, sp=0.;
107 double ct = cos(tet);
108 double st = sin(tet);
109 // Le Pi/2 echange les axes X<>Y pour theta=phi=0 !
110 // double cf = cos(phi+M_PI/2);
111 // double sf = sin(phi+M_PI/2);
112 double cf = cos(phi);
113 double sf = sin(phi);
114 double cp = 1.; // cos((double)pO);
115 double sp = 0.; // sin((double)pO);
116 RE[0][0] = cf*cp-sf*ct*sp; RE[0][1] = sf*cp+cf*ct*sp; RE[0][2] = st*sp;
117 RE[1][0] = -cf*sp-sf*ct*cp; RE[1][1] = -sf*sp+cf*ct*cp; RE[1][2] = st*cp;
118 RE[2][0] = sf*st; RE[2][1] = -cf*st; RE[2][2] = ct;
119 }
120 inline void Do(double& x, double& y)
121 {
122 double xx=x;
123 double yy=y;
124 x = RE[0][0]*xx+RE[0][1]*yy;
125 y = RE[1][0]*xx+RE[1][1]*yy;
126 }
127 double RE[3][3];
128};
129
130
131//----------------------------------------------------------------------
132// -- Pour calculer la reponse ds le plan kx,ky d'un system MultiDish
133//----------------------------------------------------------------------
134MultiDish::MultiDish(double lambda, double dmax, vector<Dish>& dishes, bool fgnoauto)
135 : lambda_(lambda), dmax_(dmax), dishes_(dishes), fgnoauto_(fgnoauto)
136{
137 SetThetaPhiRange();
138 SetRespHisNBins();
139 mcnt_=0;
140}
141
142Histo2D MultiDish::GetResponse()
143{
144 cout << " MultiDish::GetResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
145 double kmx = 1.2*DeuxPI*dmax_/lambda_;
146 double dkmx = kmx/(double)nx_;
147 double dkmy = kmx/(double)ny_;
148 double kmxx = ((double)nx_+0.5)*dkmx;
149 double kmxy = ((double)ny_+0.5)*dkmy;
150 h2w_.Define(-kmxx,kmxx,2*nx_+1,-kmxy,kmxy,2*ny_+1);
151 h2w_.SetZeroBin(0.,0.);
152
153 double dold = dishes_[0].D/lambda_;
154 double dolx = dishes_[0].Dx/lambda_;
155 double doly = dishes_[0].Dy/lambda_;
156
157 Four2DResponse rd(2, dold, dold);
158 Four2DResponse rdr(3, dolx, doly);
159
160 if (!dishes_[0].isCircular()) rd = rdr;
161
162 double dtet = thetamax_/(double)ntet_;
163 double dphi = phimax_/(double)ntet_;
164
165 double sumw = 0.;
166 for(int kt=0; kt<ntet_; kt++)
167 for(int jp=0; jp<nphi_; jp++)
168 sumw += CumulResp(rd, (double)kt*dtet, (double)jp*dphi);
169
170 double kx1 = DeuxPI*(dishes_[0].DiameterX())/lambda_;
171 double ky1 = DeuxPI*(dishes_[0].DiameterY())/lambda_;
172 int_4 ib,jb;
173 // h2w_ /= h2cnt_;
174 Histo2D h2 = h2w_.Convert();
175 h2.FindBin(kx1, ky1, ib, jb);
176 if ((kx1<0)||(ky1<0)||(kx1>=h2.NBinX())||(ky1>=h2.NBinY())) {
177 cout << " MultiDish::GetResponse[1]/ERROR kx1,ky1=" << kx1 <<","<< ky1 << " --> ib,jb=" << ib <<","<< jb << endl;
178 ib=jb=0;
179 }
180 double vmax=h2.VMax();
181 cout << " MultiDish::GetResponse[1] VMin=" << h2.VMin() << " VMax= " << vmax
182 << " h(0,0)=" << h2(0,0) << " kx1,ky1->h(" << ib <<"," << jb << ")=" << h2(ib,jb) <<endl;
183 // double fnorm=sqrt((double)dishes_.size())/h2.VMax();
184 double fnorm=1.;
185 if (vmax > sumw) {
186 fnorm=(double)dishes_.size()/h2.VMax();
187 cout << " MultiDish::GetResponse[2]/Warning h2.VMax()=" << vmax << " > sumw=" << sumw << endl;
188 cout << " ... NDishes=" << dishes_.size() << " sumw=" << sumw
189 << " Renormalizing x NDishes/VMax " << fnorm << endl;
190 }
191 else {
192 fnorm=(double)dishes_.size()/sumw;
193 cout << " MultiDish::GetResponse[3] NDishes=" << dishes_.size() << " sumw=" << sumw
194 << " Renormalizing x NDishes/sumw " << fnorm << endl;
195 }
196 h2 *= fnorm;
197 cout << " ---- MultiDish::GetResponse/[4] APRES VMin=" << h2.VMin() << " VMax= " << h2.VMax() << " h(0,0)="
198 << h2(0,0) << endl;
199 return h2;
200}
201
202Histo2D MultiDish::PosDist(int nx, int ny, double dmax)
203{
204 if (dmax<1e-3) dmax=nx*dishes_[0].Diameter();
205 double dd = dmax/nx/2.;
206 Histo2D hpos(-dd,dmax+dd,nx+1,-dd,dmax+dd,ny+1);
207 for(size_t i=0; i<NbDishes(); i++) {
208 hpos.Add(dishes_[i].X, dishes_[i].Y);
209 }
210 return hpos;
211}
212
213double MultiDish::AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh)
214{
215 double xxp = kx0+x;
216 double yyp = ky0+y;
217 double sumw=0.;
218 sumw += h2w_.Add(xxp, yyp, w, fgfh);
219 double xxm=kx0-x;
220 double yym=ky0-y;
221 // if (xxm>0.) {
222 sumw += h2w_.Add(xxm, yyp, w, fgfh);
223 // if (yym>0.)
224 sumw += h2w_.Add(xxm, yym, w, fgfh);
225 // }
226 // if (yym>0.)
227 sumw += h2w_.Add(xxp, yym, w, fgfh);
228 return sumw;
229}
230
231double MultiDish::CumulResp(Four2DResponse& rd, double theta, double phi)
232{
233 // cout << " MultiDish::CumulResp() theta=" << theta << " phi=" << phi << endl;
234
235 double dx = h2w_.WBinX()/5;
236 double dy = h2w_.WBinY()/5;
237 int nbx = DeuxPI*rd.Dx()/dx+1;
238 int nby = DeuxPI*rd.Dy()/dy+1;
239 dx = DeuxPI*rd.Dx()/(double)nbx;
240 dy = DeuxPI*rd.Dy()/(double)nby;
241 if (mcnt_==0)
242 cout << " CumulResp() nbx=" << nbx << " nby=" << nby << " dx=" << dx << " dy=" << dy << endl;
243 mcnt_++;
244
245 double sumw = 0.;
246 Rotation rot(theta, phi);
247
248 for(size_t i=0; i<dishes_.size(); i++) {
249 for(size_t j=0; j<dishes_.size(); j++) {
250 double kx0 = DeuxPI*(dishes_[i].X-dishes_[j].X)/lambda_;
251 double ky0 = DeuxPI*(dishes_[i].Y-dishes_[j].Y)/lambda_;
252 rot.Do(kx0, ky0);
253 // if (kx0<0) kx0=-kx0;
254 // if (ky0<0) ky0=-ky0;
255 bool fgfh= (!fgnoauto_ || (dishes_[i].ReflectorId()!=dishes_[j].ReflectorId()));
256 for(int ix=0; ix<nbx; ix++)
257 for(int jy=0; jy<nby; jy++) {
258 double x = ix*dx;
259 double y = jy*dy;
260 if ((ix>0)&&(jy>0)) {
261 sumw += AddToHisto(kx0, ky0, x, y, rd(x,y), fgfh);
262 }
263 else {
264 if ((ix==0)&&(jy==0))
265 sumw += h2w_.Add(kx0, ky0, rd(0.,0.), fgfh);
266 else {
267 double w = rd(x,y);
268 if (ix==0) {
269 sumw += h2w_.Add(kx0, ky0+y, w, fgfh);
270 sumw += h2w_.Add(kx0, ky0-y, w, fgfh);
271 }
272 else {
273 sumw += h2w_.Add(kx0+x, ky0, w, fgfh);
274 sumw += h2w_.Add(kx0-x, ky0, w, fgfh);
275 }
276 }
277 //
278 }
279 }
280 // if (i%10==0)
281 // cout << " MultiDish::CumulResp() done i=" << i << " / imax=" << dishes_.size()
282 // << " theta=" << theta << " phi=" << phi << endl;
283 }
284 }
285 return sumw;
286}
287
Note: See TracBrowser for help on using the repository browser.