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

Last change on this file since 3769 was 3769, checked in by ansari, 16 years ago

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