purging_magnet/src/PurgMagTabulatedField3D.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 // Code developed by:
0027 //  S.Larsson and J. Generowicz.
0028 //
0029 //    *************************************
0030 //    *                                   *
0031 //    *    PurgMagTabulatedField3D.cc     *
0032 //    *                                   *
0033 //    *************************************
0034 //
0035 //
0036
0037 #include "PurgMagTabulatedField3D.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4AutoLock.hh"
0040
0041 namespace{
0042   G4Mutex myPurgMagTabulatedField3DLock = G4MUTEX_INITIALIZER;
0043 }
0044
0045 using namespace std;
0046
0047 PurgMagTabulatedField3D::PurgMagTabulatedField3D(const char* filename,
0048                          double zOffset )
0049   :fZoffset(zOffset),invertX(false),invertY(false),invertZ(false)
0050 {
0051
0052   double lenUnit= meter;
0053   double fieldUnit= tesla;
0054   G4cout << "\n-----------------------------------------------------------"
0055      << "\n      Magnetic field"
0056      << "\n-----------------------------------------------------------";
0057
0058   G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << G4endl;
0059
0060   //
0061   //This is a thread-local class and we have to avoid that all workers open the
0062   //file at the same time
0063   G4AutoLock lock(&myPurgMagTabulatedField3DLock);
0064
0065   ifstream file( filename ); // Open the file for reading.
0066
0067   if (!file.is_open())
0068     {
0069       G4ExceptionDescription ed;
0070       ed << "Could not open input file " << filename << std::endl;
0071       G4Exception("PurgMagTabulatedField3D::PurgMagTabulatedField3D",
0072           "pugmag001",FatalException,ed);
0073     }
0074
0075   // Ignore first blank line
0076   char buffer[256];
0077   file.getline(buffer,256);
0078
0079   // Read table dimensions
0080   file >> nx >> ny >> nz; // Note dodgy order
0081
0082   G4cout << "  [ Number of values x,y,z: "
0083      << nx << " " << ny << " " << nz << " ] "
0084      << G4endl;
0085
0086   // Set up storage space for table
0087   xField.resize( nx );
0088   yField.resize( nx );
0089   zField.resize( nx );
0090   int ix, iy, iz;
0091   for (ix=0; ix<nx; ix++) {
0092     xField[ix].resize(ny);
0093     yField[ix].resize(ny);
0094     zField[ix].resize(ny);
0095     for (iy=0; iy<ny; iy++) {
0096       xField[ix][iy].resize(nz);
0097       yField[ix][iy].resize(nz);
0098       zField[ix][iy].resize(nz);
0099     }
0100   }
0101
0102   // Ignore other header information
0103   // The first line whose second character is '0' is considered to
0104   // be the last line of the header.
0105   do {
0106     file.getline(buffer,256);
0107   } while ( buffer[1]!='0');
0108
0109   // Read in the data
0110   double xval,yval,zval,bx,by,bz;
0111   double permeability; // Not used in this example.
0112   for (ix=0; ix<nx; ix++) {
0113     for (iy=0; iy<ny; iy++) {
0114       for (iz=0; iz<nz; iz++) {
0115         file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
0116         if ( ix==0 && iy==0 && iz==0 ) {
0117           minx = xval * lenUnit;
0118           miny = yval * lenUnit;
0119           minz = zval * lenUnit;
0120         }
0121         xField[ix][iy][iz] = bx * fieldUnit;
0122         yField[ix][iy][iz] = by * fieldUnit;
0123         zField[ix][iy][iz] = bz * fieldUnit;
0124       }
0125     }
0126   }
0127   file.close();
0128
0129   lock.unlock();
0130
0131   maxx = xval * lenUnit;
0132   maxy = yval * lenUnit;
0133   maxz = zval * lenUnit;
0134
0135   G4cout << "\n ---> ... done reading " << G4endl;
0136
0137   // G4cout << " Read values of field from file " << filename << G4endl;
0138   G4cout << " ---> assumed the order:  x, y, z, Bx, By, Bz "
0139      << "\n ---> Min values x,y,z: "
0140      << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
0141      << "\n ---> Max values x,y,z: "
0142      << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm "
0143      << "\n ---> The field will be offset by " << zOffset/cm << " cm " << G4endl;
0144
0145   // Should really check that the limits are not the wrong way around.
0146   if (maxx < minx) {swap(maxx,minx); invertX = true;}
0147   if (maxy < miny) {swap(maxy,miny); invertY = true;}
0148   if (maxz < minz) {swap(maxz,minz); invertZ = true;}
0149   G4cout << "\nAfter reordering if neccesary"
0150      << "\n ---> Min values x,y,z: "
0151      << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
0152      << " \n ---> Max values x,y,z: "
0153      << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm ";
0154
0155   dx = maxx - minx;
0156   dy = maxy - miny;
0157   dz = maxz - minz;
0158   G4cout << "\n ---> Dif values x,y,z (range): "
0159      << dx/cm << " " << dy/cm << " " << dz/cm << " cm in z "
0160      << "\n-----------------------------------------------------------" << G4endl;
0161 }
0162
0163 void PurgMagTabulatedField3D::GetFieldValue(const double point[4],
0164                       double *Bfield ) const
0165 {
0166
0167   double x = point[0];
0168   double y = point[1];
0169   double z = point[2] + fZoffset;
0170
0171   // Check that the point is within the defined region
0172   if ( x>=minx && x<=maxx &&
0173        y>=miny && y<=maxy &&
0174        z>=minz && z<=maxz ) {
0175
0176     // Position of given point within region, normalized to the range
0177     // [0,1]
0178     double xfraction = (x - minx) / dx;
0179     double yfraction = (y - miny) / dy;
0180     double zfraction = (z - minz) / dz;
0181
0182     if (invertX) { xfraction = 1 - xfraction;}
0183     if (invertY) { yfraction = 1 - yfraction;}
0184     if (invertZ) { zfraction = 1 - zfraction;}
0185
0186     // Need addresses of these to pass to modf below.
0187     // modf uses its second argument as an OUTPUT argument.
0188     double xdindex, ydindex, zdindex;
0189
0190     // Position of the point within the cuboid defined by the
0191     // nearest surrounding tabulated points
0192     double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
0193     double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
0194     double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));
0195
0196     // The indices of the nearest tabulated point whose coordinates
0197     // are all less than those of the given point
0198     int xindex = static_cast<int>(xdindex);
0199     int yindex = static_cast<int>(ydindex);
0200     int zindex = static_cast<int>(zdindex);
0201
0202
0203 #ifdef DEBUG_INTERPOLATING_FIELD
0204     G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << G4endl;
0205     G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << G4endl;
0206     double valx0z0, mulx0z0, valx1z0, mulx1z0;
0207     double valx0z1, mulx0z1, valx1z1, mulx1z1;
0208     valx0z0= table[xindex  ][0][zindex];  mulx0z0=  (1-xlocal) * (1-zlocal);
0209     valx1z0= table[xindex+1][0][zindex];  mulx1z0=   xlocal    * (1-zlocal);
0210     valx0z1= table[xindex  ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
0211     valx1z1= table[xindex+1][0][zindex+1]; mulx1z1=  xlocal    * zlocal;
0212 #endif
0213
0214         // Full 3-dimensional version
0215     Bfield[0] =
0216       xField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
0217       xField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
0218       xField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
0219       xField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
0220       xField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
0221       xField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
0222       xField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
0223       xField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
0224     Bfield[1] =
0225       yField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
0226       yField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
0227       yField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
0228       yField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
0229       yField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
0230       yField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
0231       yField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
0232       yField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
0233     Bfield[2] =
0234       zField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
0235       zField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
0236       zField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
0237       zField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
0238       zField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
0239       zField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
0240       zField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
0241       zField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
0242
0243   } else {
0244     Bfield[0] = 0.0;
0245     Bfield[1] = 0.0;
0246     Bfield[2] = 0.0;
0247   }
0248 }
0249