48 :fZoffset(zOffset),invertX(false),invertY(false),invertZ(false)
51 double lenUnit=
meter;
52 double fieldUnit=
tesla;
53 G4cout <<
"\n-----------------------------------------------------------"
54 <<
"\n Magnetic field"
55 <<
"\n-----------------------------------------------------------";
57 G4cout <<
"\n ---> " "Reading the field grid from " << filename <<
" ... " << endl;
62 G4AutoLock lock(&myPurgMagTabulatedField3DLock);
64 ifstream file( filename );
69 ed <<
"Could not open input file " << filename <<
G4endl;
70 G4Exception(
"PurgMagTabulatedField3D::PurgMagTabulatedField3D",
76 file.getline(buffer,256);
81 G4cout <<
" [ Number of values x,y,z: "
82 <<
nx <<
" " <<
ny <<
" " << nz <<
" ] "
90 for (ix=0; ix<
nx; ix++) {
94 for (iy=0; iy<
ny; iy++) {
105 file.getline(buffer,256);
106 }
while ( buffer[1]!=
'0');
109 double xval,yval,zval,bx,by,bz;
111 for (ix=0; ix<
nx; ix++) {
112 for (iy=0; iy<
ny; iy++) {
113 for (iz=0; iz<
nz; iz++) {
114 file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
115 if ( ix==0 && iy==0 && iz==0 ) {
116 minx = xval * lenUnit;
117 miny = yval * lenUnit;
118 minz = zval * lenUnit;
120 xField[ix][iy][
iz] = bx * fieldUnit;
121 yField[ix][iy][
iz] = by * fieldUnit;
122 zField[ix][iy][
iz] = bz * fieldUnit;
130 maxx = xval * lenUnit;
131 maxy = yval * lenUnit;
132 maxz = zval * lenUnit;
134 G4cout <<
"\n ---> ... done reading " << endl;
137 G4cout <<
" ---> assumed the order: x, y, z, Bx, By, Bz "
138 <<
"\n ---> Min values x,y,z: "
140 <<
"\n ---> Max values x,y,z: "
142 <<
"\n ---> The field will be offset by " << zOffset/
cm <<
" cm " << endl;
148 G4cout <<
"\nAfter reordering if neccesary"
149 <<
"\n ---> Min values x,y,z: "
151 <<
" \n ---> Max values x,y,z: "
157 G4cout <<
"\n ---> Dif values x,y,z (range): "
159 <<
"\n-----------------------------------------------------------" << endl;
163 double *Bfield )
const
177 double xfraction = (x -
minx) /
dx;
178 double yfraction = (y -
miny) /
dy;
179 double zfraction = (z -
minz) /
dz;
181 if (
invertX) { xfraction = 1 - xfraction;}
182 if (
invertY) { yfraction = 1 - yfraction;}
183 if (
invertZ) { zfraction = 1 - zfraction;}
187 double xdindex, ydindex, zdindex;
191 double xlocal = ( std::modf(xfraction*(
nx-1), &xdindex));
192 double ylocal = ( std::modf(yfraction*(
ny-1), &ydindex));
193 double zlocal = ( std::modf(zfraction*(
nz-1), &zdindex));
197 int xindex =
static_cast<int>(xdindex);
198 int yindex =
static_cast<int>(ydindex);
199 int zindex =
static_cast<int>(zdindex);
202 #ifdef DEBUG_INTERPOLATING_FIELD
203 G4cout <<
"Local x,y,z: " << xlocal <<
" " << ylocal <<
" " << zlocal << endl;
204 G4cout <<
"Index x,y,z: " << xindex <<
" " << yindex <<
" " << zindex << endl;
205 double valx0z0, mulx0z0, valx1z0, mulx1z0;
206 double valx0z1, mulx0z1, valx1z1, mulx1z1;
207 valx0z0= table[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
208 valx1z0= table[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
209 valx0z1= table[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
210 valx1z1= table[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
215 xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
216 xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
217 xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
218 xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
219 xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
220 xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
221 xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
222 xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
224 yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
225 yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
226 yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
227 yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
228 yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
229 yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
230 yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
231 yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
233 zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
234 zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
235 zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
236 zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
237 zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
238 zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
239 zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
240 zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
std::ostringstream G4ExceptionDescription
vector< vector< vector< double > > > zField
#define G4MUTEX_INITIALIZER
void GetFieldValue(const double Point[4], double *Bfield) const
G4GLOB_DLL std::ostream G4cout
static const double meter
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
vector< vector< vector< double > > > xField
vector< vector< vector< double > > > yField
static const double tesla
PurgMagTabulatedField3D(const char *filename, double zOffset)