Introducing Header instead of Log
[u/mrichter/AliRoot.git] / STEER / AliMagFDM.cxx
CommitLineData
aee8290b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
fb17acd4 16/* $Header$ */
aee8290b 17
8918e700 18#include <stdlib.h>
19
aee8290b 20#include "TSystem.h"
e2afb3b6 21
22#include "AliMagFDM.h"
aee8290b 23
24
25ClassImp(AliMagFDM)
26
e2afb3b6 27//_______________________________________________________________________
28AliMagFDM::AliMagFDM():
29 fSolenoid(0),
30 fInd(0),
31 fZmin(0),
32 fZmax(0),
33 fYmax(0),
34 fYmin(0),
35 fZpmx(0),
36 fZpmn(0),
37 fRmax(0),
38 fRmin(0),
39 fXdel(0),
40 fYdel(0),
41 fZdel(0),
42 fRdel(0),
43 fPhid(0),
44 fZpdl(0),
45 fCx1(0),
46 fCx2(0),
47 fAx1(0),
48 fAx2(0),
49 fXl(0),
50 fYl(0),
51 fZl(0),
52 fRn(0),
53 fPhin(0),
54 fZpl(0)
55{
56 //
57 // Default constructor for the Dipole field
58 //
59}
60
61//_______________________________________________________________________
aee8290b 62AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
e2afb3b6 63 const Float_t factor, const Float_t fmax):
64 AliMagF(name,title,integ,factor,fmax),
65 fSolenoid(0),
66 fInd(0),
67 fZmin(0),
68 fZmax(0),
69 fYmax(0),
70 fYmin(0),
71 fZpmx(0),
72 fZpmn(0),
73 fRmax(0),
74 fRmin(0),
75 fXdel(0),
76 fYdel(0),
77 fZdel(0),
78 fRdel(0),
79 fPhid(0),
80 fZpdl(0),
81 fCx1(0),
82 fCx2(0),
83 fAx1(0),
84 fAx2(0),
85 fXl(0),
86 fYl(0),
87 fZl(0),
88 fRn(0),
89 fPhin(0),
90 fZpl(0)
aee8290b 91{
ef42d733 92 //
93 // Standard constructor for the Dipole field
94 //
aee8290b 95 fType = kDipoMap;
d8408e76 96 fMap = 3;
4cc8933f 97 SetSolenoidField();
d8408e76 98
e2afb3b6 99 Info("ctor",
100 "Field Map for Muon Arm from IP till muon filter %s created: map= %d, integ= %d, factor= %f, file=%s\n",
101 fName.Data(), fMap ,integ,factor,fTitle.Data());
e678332c 102
aee8290b 103}
104
e2afb3b6 105//_______________________________________________________________________
aee8290b 106void AliMagFDM::Field(Float_t *xfi, Float_t *b)
107{
108 //
109 // Main routine to compute the field in a point
110 //
e2afb3b6 111 const Double_t keps=0.1E-06;
112 const Double_t PI2=2.*TMath::Pi();
113 const Double_t kone=1;
114
115 const Int_t kiip=33;
116 const Int_t kmiip=0;
117 const Int_t kliip=0;
118
119 const Int_t kiic=0;
120 const Int_t kmiic=0;
121 const Int_t kliic=0;
122
123 const Double_t kfZbg=502.92; // Start of Map using in z
124 const Double_t kfZL3=600; // Beginning of L3 door in z
aee8290b 125
126 Double_t x[3];
127 Double_t xL3[3];
128 Double_t bint[3];
129
130 Double_t r0;
e2afb3b6 131 Int_t iKvar,jb;
aee8290b 132
133 Double_t zp1, zp2,xp1,xp2,yp1,yp2;
134 Double_t zz1, zz2,yy1,yy2,x2,x1;
135
136// --- start the map fiel from z = 502.92 cm ---
137
138 x[0] = xfi[0];
139 x[1] = xfi[1];
140 x[2] = xfi[2];
e678332c 141 b[0]=b[1]=b[2]=0;
aee8290b 142 // printf("x[0] %f,x[1] %f,x[2] %f\n",x[0],x[1],x[2]);
143
144 Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
145 r0=rr/100;
146 Double_t rPmax;
803c9752 147 rPmax=fRmax;
148 if ( (-700<x[2] && x[2]<=kfZbg &&
aee8290b 149 (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
803c9752 150 || (kfZbg<x[2] && x[2]<=kfZL3 && (rr>rPmax*100 && rr< 560)) )
aee8290b 151 {
803c9752 152 b[0]=b[1]=0;
4cc8933f 153 b[2]=fSolenoid;
aee8290b 154 }
155
156 xL3[0]=x[0]/100;
157 xL3[1]=(x[1]+30)/100;
158 xL3[2]=x[2]/100;
803c9752 159
aee8290b 160
803c9752 161 if (TMath::Abs(xL3[0])<=0.1E-06) xL3[0]=TMath::Sign(0.1E-05,xL3[0]);
162 if (TMath::Abs(xL3[1])<=0.1E-06) xL3[1]=TMath::Sign(0.1E-05,xL3[1]);
163
164 Double_t xminn=xL3[2]*fAx1+fCx1;
165 Double_t xmaxx=xL3[2]*fAx2+fCx2;
aee8290b 166 Double_t zCmin,zCmax,yCmin,yCmax;
167
803c9752 168 zCmin=fZmin;
169 zCmax=fZmax;
170 yCmin=fYmin;
171 yCmax=fYmax;
803c9752 172if ((kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax) || ((zCmin<xL3[2] && xL3[2] <= zCmax ) && (yCmin<xL3[1] && xL3[1]<= yCmax) && (xminn<xL3[0] && xL3[0] <= xmaxx)))
aee8290b 173 {
174 if(fMap==3)
175 {
803c9752 176 if (kfZbg/100<xL3[2] && xL3[2]<=zCmin && r0<=rPmax)
aee8290b 177 {
178 //--------------------- Polar part ----------------------
179
180 Double_t yyp,ph0;
181 Int_t kp0, lp0, mp0;
182 Int_t kpi,lpi,mpi;
183 Double_t alp1,alp2,alp3;
184 Double_t zpz,rp,fip,cphi;
185
186 kpi=kiip;
187 lpi=kliip;
188 mpi=kmiip;
189
190 zpz=xL3[2];
e678332c 191 kp0=FZ(zpz, fZp ,fZpdl,kpi,fZpl) ;
192 zp2=(zpz-fZp[kp0])/(fZp[kp0+1]-fZp[kp0]);
193 zp1= kone-zp2;
aee8290b 194 yyp=xL3[1]- 0.3;
e678332c 195 cphi=TMath::Abs(yyp/r0);
196 Int_t kcphi=0;
197 if (cphi > kone) {
198 printf("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, cphi %e\n",xL3[0],xL3[1],xL3[2],yyp,r0,cphi);
199 cphi =kone;
200 kcphi=777;
201 }
aee8290b 202 ph0=TMath::ACos(cphi);
e678332c 203 if (xL3[0] < 0 && yyp > 0 ) {ph0=PI2/2 - ph0;}
204 if (xL3[0] < 0 && yyp < 0 ) {ph0=PI2/2 + ph0;}
205 if (xL3[0] > 0 && yyp < 0 ) {ph0=PI2 - ph0;}
206 if (ph0 > PI2) { ph0=ph0 - PI2;}
207 if (kcphi==777) {
208 printf("xL3[0] %e, xL3[1] %e, xL3[2] %e, yyp %e, r0 %e, ph0 %e\n",xL3[0],xL3[1],xL3[2],yyp,r0,ph0);
209 }
210 fip=ph0;
211 mp0=FZ(fip,fPhi,fPhid ,mpi,fPhin);
212 xp2=(fip-fPhi[mp0])/(fPhi[mp0+1]-fPhi[mp0]);
213 xp1=kone-xp2;
aee8290b 214
215 Double_t rDel;
803c9752 216 rDel=fRdel;
aee8290b 217
803c9752 218 if (r0<= fRdel)
aee8290b 219 {
220
221 if(r0< keps)
222 {
223
803c9752 224 bint[0]=(zp1*fBpx[kp0][0][0] + zp2*fBpx[kp0+1][0][0])*10;
225 bint[1]=(zp1*fBpy[kp0][0][0] + zp2*fBpy[kp0+1][0][0])*10;
226 bint[2]=(zp1*fBpz[kp0][0][0] + zp2*fBpz[kp0+1][0][0])*10;
aee8290b 227
803c9752 228 } else {
aee8290b 229
803c9752 230 alp2= fB[0][0][mp0]*yyp + fB[0][1][mp0]*xL3[0];
231 alp3= fB[1][0][mp0]*yyp + fB[1][1][mp0]*xL3[0];
aee8290b 232 alp1= kone - alp2 - alp3;
233
234 for (jb=0; jb<3 ; jb++)
803c9752 235 {
236 iKvar=jb;
e678332c 237 bint[jb] = Ba(iKvar,zp1,zp2,alp1,alp2,alp3, kp0,mp0)*10 ;
803c9752 238 }
239 } // end of keps <=r0
aee8290b 240 }
241 else
242 {
243 rp=r0;
244
e678332c 245 lp0=FZ(rp,fR ,fRdel,lpi,fRn);
246 yp2=(rp-fR[lp0])/(fR[lp0+1]-fR[lp0]);
247 yp1=kone-yp2;
aee8290b 248
249 for (jb=0; jb<3 ; jb++)
250 {
251 iKvar=jb;
e678332c 252 bint[jb] = Bb(zp1,zp2,yp1,yp2,xp1,xp2,iKvar,kp0,lp0,mp0)*10 ;
aee8290b 253 }
254 }
255
256 b[0]=bint[0];
257 b[1]=bint[1];
258 b[2]=bint[2];
259
aee8290b 260 }
261 else
262 {
263 //-------------- Cartensian part ------------------
264
265 Double_t zzc,yyc;
266 Int_t k0, l0,m0;
267 Double_t xx1, xx2,dx, xxx ,xXl;
268 Int_t kci,mci,lci;
269
270 kci=kiic;
271 lci=kliic;
272 mci=kmiic;
273
803c9752 274 xx1 = fAx1*xL3[2] + fCx1;
275 xx2 = fAx2*xL3[2] + fCx2;
aee8290b 276
277 zzc=xL3[2];
e678332c 278 k0=FZ(zzc, fZc ,fZdel, kci, fZl);
279 zz2=(zzc-fZc[k0])/(fZc[k0+1]-fZc[k0]);
280 zz1=kone - zz2;;
aee8290b 281
282 yyc=xL3[1];
e678332c 283 l0=FZ(yyc, fY , fYdel,lci,fYl);
284 yy2=(yyc-fY[l0])/(fY[l0+1]-fY[l0]);;
285 yy1=kone - yy2;
803c9752 286 xXl = fXl-kone;
aee8290b 287 dx = (xx2-xx1)/xXl;
288 xxx= xL3[0]-xx1;
aee8290b 289 m0 = int(xxx/dx);
290
291 if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
292 {
293 m0=m0+1;
294 printf(" m0 %d, m0+1 %d\n",m0,m0+1);
295 }
296
297 x2=(xL3[0]-( xx1+m0*dx))/dx;
298 x1=kone-x2;
299 m0=m0-1;
300 for (jb=3; jb<6; jb++)
301 {
302 iKvar=jb;
e678332c 303 bint[jb-3] = Bb(zz1,zz2,yy1,yy2,x1,x2,iKvar,k0, l0, m0)*10 ;
aee8290b 304 }
305
306 b[0]=bint[0];
307 b[1]=bint[1];
308 b[2]=bint[2];
309
aee8290b 310 }
e678332c 311
aee8290b 312
313 } else {
314 printf("Unknown map of Dipole region %d\n",fMap);
315 }
316
317} else {
318
319//This is the ZDC part
d95ea23f 320 Float_t rad2=x[0]*x[0]+x[1]*x[1];
006cb30c 321 if(x[2]>kCORBEG2 && x[2]<kCOREND2){
d95ea23f 322 if(rad2<kCOR2RA2){
323 b[0] = kFCORN2;
324 }
325 }
326 else if(x[2]>kZ1BEG && x[2]<kZ1END){
327 if(rad2<kZ1RA2){
328 b[0] = -kG1*x[1];
329 b[1] = -kG1*x[0];
330 }
331 }
332 else if(x[2]>kZ2BEG && x[2]<kZ2END){
333 if(rad2<kZ2RA2){
334 b[0] = kG1*x[1];
335 b[1] = kG1*x[0];
aee8290b 336 }
337 }
d95ea23f 338 else if(x[2]>kZ3BEG && x[2]<kZ3END){
339 if(rad2<kZ3RA2){
340 b[0] = kG1*x[1];
341 b[1] = kG1*x[0];
342 }
343 }
344 else if(x[2]>kZ4BEG && x[2]<kZ4END){
345 if(rad2<kZ4RA2){
346 b[0] = -kG1*x[1];
347 b[1] = -kG1*x[0];
348 }
349 }
350 else if(x[2]>kD1BEG && x[2]<kD1END){
351 if(rad2<kD1RA2){
352 b[1] = -kFDIP;
353 }
354 }
355 else if(x[2]>kD2BEG && x[2]<kD2END){
356 if(((x[0]-kXCEN1D2)*(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)*(x[1]-kYCEN1D2))<kD2RA2
357 || ((x[0]-kXCEN2D2)*(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)*(x[1]-kYCEN2D2))<kD2RA2){
358 b[1] = kFDIP;
359 }
360 }
361 }
362
aee8290b 363
364 if(fFactor!=1) {
365 b[0]*=fFactor;
366 b[1]*=fFactor;
367 b[2]*=fFactor;
368 }
369}
370
e2afb3b6 371//_______________________________________________________________________
372Int_t AliMagFDM::FZ(Double_t temp, Float_t *Ar,
373 Float_t delu, Int_t ik,Int_t nk)
aee8290b 374{
375 //
e678332c 376 // Quest of a point position at x,y,z (Cartensian) and R,Phi,z (Polar) axises
aee8290b 377 //
e678332c 378 Int_t l,ikj;
379 Double_t ddu,ar;
aee8290b 380
e678332c 381 Int_t ku,kf;
382 kf=0;
aee8290b 383 ar=Ar[ik];
384 ddu=temp-ar;
385
386 ku=int(ddu/delu);
e678332c 387 ikj=ik+ku+1;
aee8290b 388 if (ddu<=0) ikj=0;
389
390 for(l=ikj; l<nk; l++)
391 {
e678332c 392
803c9752 393 if(temp <= Ar[l])
aee8290b 394 {
803c9752 395
e678332c 396 kf=l-1;
397 break;
aee8290b 398 }
399 }
e678332c 400 return kf;
aee8290b 401 }
402
e2afb3b6 403//_______________________________________________________________________
404Double_t AliMagFDM::Ba(Int_t kaai,Double_t zaa1, Double_t zaa2,
405 Double_t alf1, Double_t alf2, Double_t alf3,
406 Int_t kaa, Int_t maa)
aee8290b 407{
408 //
e678332c 409 // Calculation of field componet for case (keps <r0<= fRdel) at a given axis
aee8290b 410 //
411 Double_t fa11,fa12,fa13;
412 Double_t fa21,fa22,fa23;
413 Double_t faY1,faY2;
414 Double_t bba;
e678332c 415
8918e700 416 switch (kaai) {
417 case 0:
803c9752 418 fa11 = fBpx[kaa][0][0];
419 fa12 = fBpx[kaa][0][maa];
420 fa13 = fBpx[kaa][0][maa+1];
421 fa21 = fBpx[kaa+1][0][0];
422 fa22 = fBpx[kaa+1][0][maa];
423 fa23 = fBpx[kaa+1][0][maa+1];
8918e700 424 break;
425 case 1:
803c9752 426 fa11 = fBpy[kaa][0][0];
427 fa12 = fBpy[kaa][0][maa];
428 fa13 = fBpy[kaa][0][maa+1];
429 fa21 = fBpy[kaa+1][0][0];
430 fa22 = fBpy[kaa+1][0][maa];
431 fa23 = fBpy[kaa+1][0][maa+1];
8918e700 432 break;
433 case 2:
803c9752 434 fa11 = fBpz[kaa][0][0];
435 fa12 = fBpz[kaa][0][maa];
436 fa13 = fBpz[kaa][0][maa+1];
437 fa21 = fBpz[kaa+1][0][0];
438 fa22 = fBpz[kaa+1][0][maa];
439 fa23 = fBpz[kaa+1][0][maa+1];
8918e700 440 break;
441 default:
e678332c 442 Fatal("Ba","Invalid value of kaai %d\n",kaai);
8918e700 443 exit(1);
aee8290b 444 }
445 faY1=alf1*fa11+alf2*fa12+alf3*fa13;
446 faY2=alf1*fa21+alf2*fa22+alf3*fa23;
447 bba = zaa1*faY1+zaa2*faY2;
e678332c 448 return bba;
aee8290b 449}
450
451
e2afb3b6 452//_______________________________________________________________________
453Double_t AliMagFDM::Bb(Double_t z1,Double_t z2, Double_t y1,Double_t y2,
454 Double_t x1,Double_t x2, Int_t kv, Int_t k, Int_t l, Int_t m)
aee8290b 455{
456 //
e678332c 457 // Calculation of field componet at a given axis (general case)
aee8290b 458 //
459 Double_t fy1, fy2, ffy;
460 Double_t gy1,gy2,ggy;
aee8290b 461 Double_t bbi;
aee8290b 462 Double_t bf11,bf12,bf21,bf22;
463 Double_t bg11,bg12,bg21,bg22;
e678332c 464
aee8290b 465
466 /*-----------------Polar part ------------------*/
467
8918e700 468 switch (kv) {
469 case 0:
803c9752 470 bf11=fBpx[k][l][m];
471 bf12=fBpx[k+1][l][m];
472 bf21=fBpx[k+1][l+1][m];
473 bf22=fBpx[k][l+1][m];
aee8290b 474
803c9752 475 bg11=fBpx[k][l][m+1];
476 bg12=fBpx[k+1][l][m+1];
477 bg21=fBpx[k+1][l+1][m+1];
478 bg22=fBpx[k][l+1][m+1];
8918e700 479 break;
480
481 case 1:
803c9752 482 bf11=fBpy[k][l][m];
483 bf12=fBpy[k+1][l][m];
484 bf21=fBpy[k+1][l+1][m];
485 bf22=fBpy[k][l+1][m];
aee8290b 486
803c9752 487 bg11=fBpy[k][l][m+1];
488 bg12=fBpy[k+1][l][m+1];
489 bg21=fBpy[k+1][l+1][m+1];
490 bg22=fBpy[k][l+1][m+1];
8918e700 491 break;
492
493 case 2:
803c9752 494 bf11=fBpz[k][l][m];
495 bf12=fBpz[k+1][l][m];
496 bf21=fBpz[k+1][l+1][m];
497 bf22=fBpz[k][l+1][m];
aee8290b 498
803c9752 499 bg11=fBpz[k][l][m+1];
500 bg12=fBpz[k+1][l][m+1];
501 bg21=fBpz[k+1][l+1][m+1];
502 bg22=fBpz[k][l+1][m+1];
8918e700 503 break;
aee8290b 504 /*-----------------Cartensian part ---------------*/
505
8918e700 506 case 3:
803c9752 507 bf11=fBcx[k][l][m];
508 bf12=fBcx[k+1][l][m];
509 bf21=fBcx[k+1][l+1][m];
510 bf22=fBcx[k][l+1][m];
aee8290b 511
803c9752 512 bg11=fBcx[k][l][m+1];
513 bg12=fBcx[k+1][l][m+1];
514 bg21=fBcx[k+1][l+1][m+1];
515 bg22=fBcx[k][l+1][m+1];
8918e700 516 break;
517
518 case 4:
803c9752 519 bf11=fBcy[k][l][m];
520 bf12=fBcy[k+1][l][m];
521 bf21=fBcy[k+1][l+1][m];
522 bf22=fBcy[k][l+1][m];
aee8290b 523
803c9752 524 bg11=fBcy[k][l][m+1];
525 bg12=fBcy[k+1][l][m+1];
526 bg21=fBcy[k+1][l+1][m+1];
527 bg22=fBcy[k][l+1][m+1];
8918e700 528 break;
529
530 case 5:
803c9752 531 bf11=fBcz[k][l][m];
532 bf12=fBcz[k+1][l][m];
533 bf21=fBcz[k+1][l+1][m];
534 bf22=fBcz[k][l+1][m];
aee8290b 535
803c9752 536 bg11=fBcz[k][l][m+1];
537 bg12=fBcz[k+1][l][m+1];
538 bg21=fBcz[k+1][l+1][m+1];
539 bg22=fBcz[k][l+1][m+1];
8918e700 540 break;
541
542 default:
e678332c 543 Fatal("Bb","Invalid value of kv %d\n",kv);
8918e700 544 exit(1);
aee8290b 545 }
546
547
548
549 fy1=z1*bf11+z2*bf12;
550 fy2=z2*bf21+z1* bf22;
551 ffy=y1*fy1+ y2*fy2;
552
553
554 gy1 = z1*bg11+z2*bg12;
555 gy2 = z2*bg21+z1*bg22;
556 ggy= y1*gy1 + y2*gy2;
557
558 bbi = x1*ffy+x2*ggy;
e678332c 559
560 return bbi;
aee8290b 561
562}
aee8290b 563
e2afb3b6 564//_______________________________________________________________________
aee8290b 565void AliMagFDM::ReadField()
566{
567 //
568 // Method to read the magnetic field from file
569 //
570 FILE *magfile;
571
572 Int_t ik, il, im;
573 Float_t zzp, rr,phii;
574 Float_t zz, yy, bx,by,bz,bb;
575
576 char *fname;
577 printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
578 fname = gSystem->ExpandPathName(fTitle.Data());
579 magfile=fopen(fname,"r");
580 delete [] fname;
581
582 printf("Cartensian part\n");
583
584 if (magfile) {
585
586// Cartensian part
587
803c9752 588 fscanf(magfile,"%d %d %d ",&fYl, &fXl, &fZl);
aee8290b 589
803c9752 590 printf("fYl %d, fXl %d, fZl %d\n",fYl, fXl, fZl);
aee8290b 591
803c9752 592 for (ik=0; ik<fZl; ik++)
aee8290b 593 {
594
595 fscanf(magfile, " %e ", &zz);
803c9752 596 fZc[ik]=zz;
aee8290b 597
598 }
599
803c9752 600 for (ik=0; ik<fYl; ik++)
aee8290b 601 {
602 fscanf(magfile, " %e ", &yy);
803c9752 603 fY[ik]=yy;
aee8290b 604
605 }
606 for (ik=0; ik<81; ik++)
607 {
803c9752 608 printf("fZc %e,fY %e\n", fZc[ik],fY[ik]);
aee8290b 609 }
610
803c9752 611 fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fYdel,&fXdel,&fZdel,&fZmax,&fZmin,&fYmax,&fYmin,&fAx1,&fCx1,&fAx2,&fCx2);
aee8290b 612
803c9752 613printf("fYdel %e, fXdel %e, fZdel %e\n",fYdel,fXdel,fZdel);
614printf("fZmax %e, fZmin %e, fYmax %e,fYmin %e\n",fZmax,fZmin,fYmax,fYmin);
615printf("fAx1 %e, fCx1 %e, fAx2 %e, fCx %e\n",fAx1,fCx1,fAx2,fCx2);
aee8290b 616
617 for (il=0; il<44; il++) {
618 for (im=0; im<81; im++) {
619 for (ik=0; ik<81; ik++) {
620
621 fscanf(magfile, " %e ", &by);
803c9752 622 fBcy[ik][im][il]=by;
aee8290b 623 }
624 }
625 }
626
627 for (il=0; il<44; il++) {
628 for (im=0; im<81; im++) {
629 for (ik=0; ik<81; ik++) {
630
631 fscanf(magfile, " %e ", &bx);
803c9752 632 fBcx[ik][im][il]=bx;
aee8290b 633 }
634 }
635 }
636
637 for (il=0; il<44; il++) {
638 for (im=0; im<81; im++) {
639 for (ik=0; ik<81; ik++) {
640
641 fscanf(magfile, " %e ", &bz);
803c9752 642 fBcz[ik][im][il]=bz;
aee8290b 643 }
644 }
645 }
646//---------------------- Polar part ---------------------------------
647
648 printf("Polar part\n");
803c9752 649 fscanf(magfile,"%d %d %d ", &fZpl, &fRn, &fPhin);
650 printf("fZpl %d, fRn %d, fPhin %d\n",fZpl,fRn,fPhin);
aee8290b 651
803c9752 652 printf(" fZp array\n");
aee8290b 653
654 for (ik=0; ik<51; ik++)
655 {
656 fscanf(magfile, " %e ", &zzp);
803c9752 657 fZp[ik]=zzp;
658 printf(" %e\n",fZp[ik]);
aee8290b 659 }
660
803c9752 661 printf(" fR array\n");
aee8290b 662
663 for (ik=0; ik<10; ik++)
664 {
665 fscanf(magfile, " %e ", &rr);
803c9752 666 fR[ik]=rr;
667 printf(" %e\n",fR[ik]);
aee8290b 668 }
669
803c9752 670// printf("fPhi array\n");
aee8290b 671
672 for (il=0; il<33; il++)
673 {
674 fscanf(magfile, " %e ", &phii);
803c9752 675 fPhi[il]=phii;
676// printf(" %e\n",fPhi[il]);
aee8290b 677 }
678
803c9752 679 fscanf(magfile," %e %e %e %e %e %e %e ",&fZpdl,&fPhid,&fRdel,&fZpmx,&fZpmn,&fRmax, &fRmin);
aee8290b 680
803c9752 681printf("fZpdl %e, fPhid %e, fRdel %e, fZpmx %e, fZpmn %e,fRmax %e,fRmin %e \n", fZpdl,fPhid, fRdel,fZpmx, fZpmn,fRmax, fRmin);
aee8290b 682
683
684 for (il=0; il<33; il++) {
685 for (im=0; im<10; im++) {
686 for (ik=0; ik<51; ik++) {
687 fscanf(magfile, " %e ", &by);
803c9752 688 fBpy[ik][im][il]=by;
aee8290b 689 }
690 }
691 }
692
693 for (il=0; il<33; il++) {
694 for (im=0; im<10; im++) {
695 for (ik=0; ik<51; ik++) {
696 fscanf(magfile, " %e ", &bx);
803c9752 697 fBpx[ik][im][il]=bx;
aee8290b 698 }
699 }
700 }
701
702
703 for (il=0; il<33; il++) {
704 for (im=0; im<10; im++) {
705 for (ik=0; ik<51; ik++) {
706 fscanf(magfile, " %e ", &bz);
803c9752 707 fBpz[ik][im][il]=bz;
aee8290b 708 }
709 }
710 }
711
712
713 for (il=0; il<32; il++) {
714 for (im=0; im<2; im++) {
715 for (ik=0; ik<2; ik++) {
716 fscanf(magfile, " %e ", &bb);
803c9752 717 fB[ik][im][il]=bb;
aee8290b 718 }
719 }
720 }
721//
722 } else {
723 printf("File %s not found !\n",fTitle.Data());
724 exit(1);
725 }
726}