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