nanobeam/src/TabulatedField3D.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Please cite the following papers if you use this software
0027 // Nucl.Instrum.Meth.B260:20-27, 2007
0028 // IEEE TNS 51, 4:1395-1401, 2004
0029 //
0030 // Based on purging magnet advanced example
0031 //
0032
0033 #include "TabulatedField3D.hh"
0034 #include "G4SystemOfUnits.hh"
0035 #include "G4Exp.hh"
0036
0037 #include "G4AutoLock.hh"
0038
0039 namespace
0040 {
0041   G4Mutex myTabulatedField3DLock = G4MUTEX_INITIALIZER;
0042 }
0043
0044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0045
0046 TabulatedField3D::TabulatedField3D(G4float gr1, G4float gr2, G4float gr3, G4float gr4, G4int choiceModel)
0047 {
0048
0049   G4cout << " ********************** " << G4endl;
0050   G4cout << " **** CONFIGURATION *** " << G4endl;
0051   G4cout << " ********************** " << G4endl;
0052
0053   G4cout<< G4endl;
0054   G4cout << "=====> You have selected :" << G4endl;
0055   if (choiceModel==1) G4cout<< "-> Square quadrupole field"<<G4endl;
0056   if (choiceModel==2) G4cout<< "-> 3D quadrupole field"<<G4endl;
0057   if (choiceModel==3) G4cout<< "-> Enge quadrupole field"<<G4endl;
0058   G4cout << "   G1 (T/m) = "<< gr1 << G4endl;
0059   G4cout << "   G2 (T/m) = "<< gr2 << G4endl;
0060   G4cout << "   G3 (T/m) = "<< gr3 << G4endl;
0061   G4cout << "   G4 (T/m) = "<< gr4 << G4endl;
0062
0063   fGradient1 = gr1;
0064   fGradient2 = gr2;
0065   fGradient3 = gr3;
0066   fGradient4 = gr4;
0067   fModel = choiceModel;
0068
0069   if (fModel==2)
0070   {
0071   //
0072   //This is a thread-local class and we have to avoid that all workers open the
0073   //file at the same time
0074   G4AutoLock lock(&myTabulatedField3DLock);
0075   //
0076
0077   const char * filename ="OM50.grid";
0078
0079   double lenUnit= mm;
0080   G4cout << "\n-----------------------------------------------------------"
0081      << "\n      3D Magnetic field from OPERA software "
0082      << "\n-----------------------------------------------------------";
0083
0084   G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << G4endl;
0085   std::ifstream file( filename ); // Open the file for reading.
0086
0087   // Read table dimensions
0088   file >> fNx >> fNy >> fNz; // Note dodgy order
0089
0090   G4cout << "  [ Number of values x,y,z: "
0091      << fNx << " " << fNy << " " << fNz << " ] "
0092      << G4endl;
0093
0094   // Set up storage space for table
0095   fXField.resize( fNx );
0096   fYField.resize( fNx );
0097   fZField.resize( fNx );
0098   int ix, iy, iz;
0099   for (ix=0; ix<fNx; ix++)
0100   {
0101     fXField[ix].resize(fNy);
0102     fYField[ix].resize(fNy);
0103     fZField[ix].resize(fNy);
0104     for (iy=0; iy<fNy; iy++)
0105     {
0106       fXField[ix][iy].resize(fNz);
0107       fYField[ix][iy].resize(fNz);
0108       fZField[ix][iy].resize(fNz);
0109     }
0110   }
0111
0112   // Read in the data
0113   double xval,yval,zval,bx,by,bz;
0114   double permeability; // Not used in this example.
0115   for (ix=0; ix<fNx; ix++)
0116   {
0117     for (iy=0; iy<fNy; iy++)
0118     {
0119       for (iz=0; iz<fNz; iz++)
0120       {
0121         file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
0122         if ( ix==0 && iy==0 && iz==0 )
0123     {
0124           fMinix = xval * lenUnit;
0125           fMiniy = yval * lenUnit;
0126           fMiniz = zval * lenUnit;
0127         }
0128         fXField[ix][iy][iz] = bx ;
0129         fYField[ix][iy][iz] = by ;
0130         fZField[ix][iy][iz] = bz ;
0131       }
0132     }
0133   }
0134   file.close();
0135
0136   //
0137   lock.unlock();
0138   //
0139
0140   fMaxix = xval * lenUnit;
0141   fMaxiy = yval * lenUnit;
0142   fMaxiz = zval * lenUnit;
0143
0144   G4cout << "\n ---> ... done reading " << G4endl;
0145
0146   // G4cout << " Read values of field from file " << filename << G4endl;
0147
0148   G4cout << " ---> assumed the order:  x, y, z, Bx, By, Bz "
0149      << "\n ---> Min values x,y,z: "
0150      << fMinix/cm << " " << fMiniy/cm << " " << fMiniz/cm << " cm "
0151      << "\n ---> Max values x,y,z: "
0152      << fMaxix/cm << " " << fMaxiy/cm << " " << fMaxiz/cm << " cm " << G4endl;
0153
0154   fDx = fMaxix - fMinix;
0155   fDy = fMaxiy - fMiniy;
0156   fDz = fMaxiz - fMiniz;
0157
0158   G4cout << "\n ---> Dif values x,y,z (range): "
0159      << fDx/cm << " " << fDy/cm << " " << fDz/cm << " cm in z "
0160      << "\n-----------------------------------------------------------" << G4endl;
0161
0162
0163   // Table normalization
0164
0165   for (ix=0; ix<fNx; ix++)
0166   {
0167     for (iy=0; iy<fNy; iy++)
0168     {
0169       for (iz=0; iz<fNz; iz++)
0170       {
0171
0172     fXField[ix][iy][iz] = (fXField[ix][iy][iz]/197.736);
0173         fYField[ix][iy][iz] = (fYField[ix][iy][iz]/197.736);
0174         fZField[ix][iy][iz] = (fZField[ix][iy][iz]/197.736);
0175
0176       }
0177     }
0178   }
0179
0180   } // fModel==2
0181
0182 }
0183
0184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0185
0186
0187 void TabulatedField3D::GetFieldValue(const double point[4],
0188                       double *Bfield ) const
0189 {
0190   //G4cout << fGradient1 << G4endl;
0191   //G4cout << fGradient2 << G4endl;
0192   //G4cout << fGradient3 << G4endl;
0193   //G4cout << fGradient4 << G4endl;
0194   //G4cout << "---------" << G4endl;
0195
0196   G4double coef, G0;
0197   G0 = 0;
0198
0199   coef=1; // for protons
0200   //coef=2; // for alphas
0201
0202 //******************************************************************
0203
0204 // MAP
0205
0206 if (fModel==2)
0207 {
0208   Bfield[0] = 0.0;
0209   Bfield[1] = 0.0;
0210   Bfield[2] = 0.0;
0211   Bfield[3] = 0.0;
0212   Bfield[4] = 0.0;
0213   Bfield[5] = 0.0;
0214
0215   double x = point[0];
0216   double y = point[1];
0217   double z = point[2];
0218
0219   G4int quad;
0220   G4double gradient[5];
0221
0222   gradient[0]=fGradient1*(tesla/m)/coef;
0223   gradient[1]=fGradient2*(tesla/m)/coef;
0224   gradient[2]=fGradient3*(tesla/m)/coef;
0225   gradient[3]=fGradient4*(tesla/m)/coef;
0226   gradient[4]=-fGradient3*(tesla/m)/coef;
0227
0228   for (quad=0; quad<=4; quad++)
0229   {
0230   if ((quad+1)==1) {z = point[2] + 3720 * mm;}
0231   if ((quad+1)==2) {z = point[2] + 3580 * mm;}
0232   if ((quad+1)==3) {z = point[2] + 330  * mm;}
0233   if ((quad+1)==4) {z = point[2] + 190  * mm;}
0234   if ((quad+1)==5) {z = point[2] + 50   * mm;}
0235
0236   // Check that the point is within the defined region
0237
0238   if
0239   (
0240     x>=fMinix && x<=fMaxix &&
0241     y>=fMiniy && y<=fMaxiy &&
0242     z>=fMiniz && z<=fMaxiz
0243   )
0244   {
0245     // Position of given point within region, normalized to the range
0246     // [0,1]
0247     double xfraction = (x - fMinix) / fDx;
0248     double yfraction = (y - fMiniy) / fDy;
0249     double zfraction = (z - fMiniz) / fDz;
0250
0251     // Need addresses of these to pass to modf below.
0252     // modf uses its second argument as an OUTPUT argument.
0253     double xdindex, ydindex, zdindex;
0254
0255     // Position of the point within the cuboid defined by the
0256     // nearest surrounding tabulated points
0257     double xlocal = ( std::modf(xfraction*(fNx-1), &xdindex));
0258     double ylocal = ( std::modf(yfraction*(fNy-1), &ydindex));
0259     double zlocal = ( std::modf(zfraction*(fNz-1), &zdindex));
0260
0261     // The indices of the nearest tabulated point whose coordinates
0262     // are all less than those of the given point
0263
0264     //int xindex = static_cast<int>(xdindex);
0265     //int yindex = static_cast<int>(ydindex);
0266     //int zindex = static_cast<int>(zdindex);
0267
0268     // SI 15/12/2016: modified according to bugzilla 1879
0269     int xindex = static_cast<int>(std::floor(xdindex));
0270     int yindex = static_cast<int>(std::floor(ydindex));
0271     int zindex = static_cast<int>(std::floor(zdindex));
0272
0273     // Interpolated field
0274     Bfield[0] =
0275      (fXField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
0276       fXField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
0277       fXField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
0278       fXField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
0279       fXField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
0280       fXField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
0281       fXField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
0282       fXField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal)*gradient[quad]
0283       + Bfield[0];
0284
0285     Bfield[1] =
0286      (fYField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
0287       fYField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
0288       fYField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
0289       fYField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
0290       fYField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
0291       fYField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
0292       fYField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
0293       fYField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal)*gradient[quad]
0294       + Bfield[1];
0295
0296     Bfield[2] =
0297      (fZField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
0298       fZField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
0299       fZField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
0300       fZField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
0301       fZField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
0302       fZField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
0303       fZField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
0304       fZField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal)*gradient[quad]
0305       + Bfield[2];
0306
0307      }
0308
0309 } // loop on quads
0310
0311 } //end  MAP
0312
0313
0314 //******************************************************************
0315 // SQUARE
0316
0317 if (fModel==1)
0318 {
0319   Bfield[0] = 0.0;
0320   Bfield[1] = 0.0;
0321   Bfield[2] = 0.0;
0322   Bfield[3] = 0.0;
0323   Bfield[4] = 0.0;
0324   Bfield[5] = 0.0;
0325
0326   // Field components
0327   G4double Bx = 0;
0328   G4double By = 0;
0329   G4double Bz = 0;
0330
0331   G4double x = point[0];
0332   G4double y = point[1];
0333   G4double z = point[2];
0334
0335   if (z>=-3770*mm && z<=-3670*mm)  G0 = (fGradient1/coef)* tesla/m;
0336   if (z>=-3630*mm && z<=-3530*mm)  G0 = (fGradient2/coef)* tesla/m;
0337
0338   if (z>=-380*mm  && z<=-280*mm)   G0 = (fGradient3/coef)* tesla/m;
0339   if (z>=-240*mm  && z<=-140*mm)   G0 = (fGradient4/coef)* tesla/m;
0340   if (z>=-100*mm  && z<=0*mm)      G0 = (-fGradient3/coef)* tesla/m;
0341
0342   Bx = y*G0;
0343   By = x*G0;
0344   Bz = 0;
0345
0346   Bfield[0] = Bx;
0347   Bfield[1] = By;
0348   Bfield[2] = Bz;
0349
0350 }
0351
0352 // end SQUARE
0353
0354 //******************************************************************
0355 // ENGE
0356
0357 if (fModel==3)
0358 {
0359
0360   // X POSITION OF FIRST QUADRUPOLE
0361   // G4double lineX = 0*mm;
0362
0363   // Z POSITION OF FIRST QUADRUPOLE
0364   G4double lineZ = -3720*mm;
0365
0366   // QUADRUPOLE HALF LENGTH
0367   // G4double quadHalfLength = 50*mm;
0368
0369   // QUADRUPOLE CENTER COORDINATES
0370   G4double zoprime;
0371
0372   G4double Grad1, Grad2, Grad3, Grad4, Grad5;
0373   Grad1=fGradient1;
0374   Grad2=fGradient2;
0375   Grad3=fGradient3;
0376   Grad4=fGradient4;
0377   Grad5=-Grad3;
0378
0379   Bfield[0] = 0.0;
0380   Bfield[1] = 0.0;
0381   Bfield[2] = 0.0;
0382   Bfield[3] = 0.0;
0383   Bfield[4] = 0.0;
0384   Bfield[5] = 0.0;
0385
0386   double x = point[0];
0387   double y = point[1];
0388   double z = point[2];
0389
0390   if ( (z>=-3900*mm && z<-3470*mm)  || (z>=-490*mm && z<100*mm) )
0391   {
0392   G4double Bx=0;
0393   G4double By=0;
0394   G4double Bz=0;
0395
0396   // FRINGING FILED CONSTANTS
0397   G4double c0[5], c1[5], c2[5], z1[5], z2[5], a0[5], gradient[5];
0398
0399   // DOUBLET***************
0400
0401   // QUADRUPOLE 1
0402   c0[0] = -10.;         // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
0403   c1[0] = 3.08874;
0404   c2[0] = -0.00618654;
0405   z1[0] = 28.6834*mm;       // Fringing field lower limit
0406   z2[0] = z1[0]+50*mm;      // Fringing field upper limit
0407   a0[0] = 7.5*mm;               // Bore Radius
0408   gradient[0] =Grad1*(tesla/m)/coef;
0409
0410   // QUADRUPOLE 2
0411   c0[1] = -10.;         // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
0412   c1[1] = 3.08874;
0413   c2[1] = -0.00618654;
0414   z1[1] = 28.6834*mm;       // Fringing field lower limit
0415   z2[1] = z1[1]+50*mm;      // Fringing field upper limit
0416   a0[1] = 7.5*mm;               // Bore Radius
0417   gradient[1] =Grad2*(tesla/m)/coef;
0418
0419   // TRIPLET**********
0420
0421   // QUADRUPOLE 3
0422   c0[2] = -10.;         // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
0423   c1[2] = 3.08874;
0424   c2[2] = -0.00618654;
0425   z1[2] = 28.6834*mm;       // Fringing field lower limit
0426   z2[2] = z1[2]+50*mm;      // Fringing field upper limit
0427   a0[2] = 7.5*mm;               // Bore Radius
0428   gradient[2] = Grad3*(tesla/m)/coef;
0429
0430   // QUADRUPOLE 4
0431   c0[3] = -10.;         // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
0432   c1[3] = 3.08874;
0433   c2[3] = -0.00618654;
0434   z1[3] = 28.6834*mm;       // Fringing field lower limit
0435   z2[3] = z1[3]+50*mm;      // Fringing field upper limit
0436   a0[3] = 7.5*mm;               // Bore Radius
0437   gradient[3] = Grad4*(tesla/m)/coef;
0438
0439    // QUADRUPOLE 5
0440   c0[4] = -10.;         // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
0441   c1[4] = 3.08874;
0442   c2[4] = -0.00618654;
0443   z1[4] = 28.6834*mm;       // Fringing field lower limit
0444   z2[4] = z1[4]+50*mm;      // Fringing field upper limit
0445   a0[4] = 7.5*mm;               // Bore Radius
0446   gradient[4] = Grad5*(tesla/m)/coef;
0447
0448   // FIELD CREATED BY A QUADRUPOLE IN ITS LOCAL FRAME
0449   G4double Bx_local,By_local,Bz_local;
0450   Bx_local = 0; By_local = 0; Bz_local = 0;
0451
0452   // QUADRUPOLE FRAME
0453   G4double x_local,y_local,z_local;
0454   x_local= 0; y_local=0; z_local=0;
0455
0456   G4double myVars = 0;          // For Enge formula
0457   G4double G1, G2, G3;          // For Enge formula
0458   G4double K1, K2, K3;      // For Enge formula
0459   G4double P0, P1, P2, cte; // For Enge formula
0460
0461   K1=0;
0462   K2=0;
0463   K3=0;
0464
0465   P0=0;
0466   P1=0;
0467   P2=0;
0468
0469   G0=0;
0470   G1=0;
0471   G2=0;
0472   G3=0;
0473
0474   cte=0;
0475
0476   for (G4int i=0;i<5; i++) // LOOP ON MAGNETS
0477   {
0478
0479      if (i<2) // (if Doublet)
0480      {
0481        zoprime = lineZ + i*140*mm; // centre of magnet nbr i
0482        x_local = x;
0483        y_local = y;
0484        z_local = (z - zoprime);
0485      }
0486      else    // else the current magnet is in the triplet
0487      {
0488        zoprime = lineZ + i*140*mm +(3150-40)*mm;
0489
0490        x_local = x;
0491        y_local = y;
0492        z_local = (z - zoprime);
0493
0494      }
0495
0496      if ( z_local < -z2[i] || z_local > z2[i])  // Outside the fringing field
0497      {
0498       G0=0;
0499       G1=0;
0500       G2=0;
0501       G3=0;
0502      }
0503
0504      if ( (z_local>=-z1[i]) && (z_local<=z1[i]) ) // inside the quadrupole but outside the fringefield
0505      {
0506       G0=gradient[i];
0507       G1=0;
0508       G2=0;
0509       G3=0;
0510      }
0511
0512      if ( ((z_local>=-z2[i]) && (z_local<-z1[i])) ||  ((z_local>z1[i]) && (z_local<=z2[i])) ) // inside the fringefield
0513      {
0514
0515           myVars = ( z_local - z1[i]) / a0[i];     // se (8) p1397 TNS 51
0516           if (z_local<-z1[i])  myVars = ( - z_local - z1[i]) / a0[i];  // see (9) p1397 TNS 51
0517
0518
0519       P0 = c0[i]+c1[i]*myVars+c2[i]*myVars*myVars;
0520
0521       P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i]; // dP/fDz
0522       if (z_local<-z1[i])  P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i];
0523
0524       P2 = 2*c2[i]/a0[i]/a0[i];     // d2P/fDz2
0525
0526       cte = 1 + G4Exp(c0[i]);   // (1+e^c0)
0527
0528       K1 = -cte*P1*G4Exp(P0)/( (1+G4Exp(P0))*(1+G4Exp(P0)) );  // see (11) p1397 TNS 51
0529
0530       K2 = -cte*G4Exp(P0)*(                 // see (12) p1397 TNS 51
0531        P2/( (1+G4Exp(P0))*(1+G4Exp(P0)) )
0532        +2*P1*K1/(1+G4Exp(P0))/cte
0533        +P1*P1/(1+G4Exp(P0))/(1+G4Exp(P0))
0534        );
0535
0536       K3 = -cte*G4Exp(P0)*(             // see (13) p1397 TNS 51
0537        (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp(P0))
0538        +4*K1*(P1*P1+P2)/(1+G4Exp(P0))/cte
0539        +2*P1*(K1*K1/cte/cte+K2/(1+G4Exp(P0))/cte)
0540        );
0541
0542       G0 = gradient[i]*cte/(1+G4Exp(P0));       // G = G0*K(z) see (7) p1397 TNS 51
0543       G1 = gradient[i]*K1;              // dG/fDz
0544       G2 = gradient[i]*K2;              // d2G/fDz2
0545       G3 = gradient[i]*K3;              // d3G/fDz3
0546
0547      }
0548
0549      Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2);    // see (4) p1396 TNS 51
0550      By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2);    // see (5) p1396 TNS 51
0551      Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3);  // see (6) p1396 TNS 51
0552
0553      // TOTAL MAGNETIC FIELD
0554
0555      Bx = Bx + Bx_local ;
0556      By = By + By_local ;
0557      Bz = Bz + Bz_local ;
0558
0559
0560   } // LOOP ON QUADRUPOLES
0561
0562   Bfield[0] = Bx;
0563   Bfield[1] = By;
0564   Bfield[2] = Bz;
0565   }
0566
0567
0568 } // end ENGE
0569
0570 }