42 double lenUnit=
meter;
43 double fieldUnit=
tesla;
44 G4cout <<
"\n-----------------------------------------------------------"
45 <<
"\n Magnetic field"
46 <<
"\n-----------------------------------------------------------";
49 G4cout <<
"\n ---> " "Reading the field grid from " << filename <<
" ... " <<
G4endl;
52 ifstream
file( filename );
56 file.getline(buffer,256);
61 G4cout <<
" [ Number of values x,y,z: "
62 <<
nx <<
" " <<
ny <<
" " << nz <<
" ] "
70 for (ix=0; ix<
nx; ix++) {
74 for (iy=0; iy<
ny; iy++) {
88 for (ix=0; ix<
nx; ix++) {
89 for (iy=0; iy<
ny; iy++) {
90 for (iz=0; iz<
nz; iz++) {
91 file >> xval >> yval >> zval >> bx >> by >> bz ;
92 if ( ix==0 && iy==0 && iz==0 ) {
93 minx = xval * lenUnit;
94 miny = yval * lenUnit;
95 minz = zval * lenUnit;
97 xField[ix][iy][iz] = bx * fieldUnit;
98 yField[ix][iy][iz] = by * fieldUnit;
99 zField[ix][iy][iz] = bz * fieldUnit;
107 maxx = xval * lenUnit;
108 maxy = yval * lenUnit;
109 maxz = zval * lenUnit;
114 G4cout <<
" ---> assumed the order: x, y, z, Bx, By, Bz "
115 <<
"\n ---> Min values x,y,z: "
117 <<
"\n ---> Max values x,y,z: "
119 <<
"\n ---> The field will be offset by " << xOffset/
cm <<
" cm " <<
G4endl;
125 G4cout <<
"\nAfter reordering if neccesary"
126 <<
"\n ---> Min values x,y,z: "
128 <<
" \n ---> Max values x,y,z: "
134 G4cout <<
"\n ---> Dif values x,y,z (range): "
136 <<
"\n-----------------------------------------------------------" <<
G4endl;
140 double *Bfield )
const
148 double xfraction = (x -
minx) /
dx;
149 double yfraction = (y -
miny) /
dy;
150 double zfraction = (z -
minz) /
dz;
152 if (
invertX) { xfraction = 1 - xfraction;}
153 if (
invertY) { yfraction = 1 - yfraction;}
154 if (
invertZ) { zfraction = 1 - zfraction;}
158 double xdindex, ydindex, zdindex;
162 double xlocal = ( std::modf(xfraction*(
nx-1), &xdindex));
163 double ylocal = ( std::modf(yfraction*(
ny-1), &ydindex));
164 double zlocal = ( std::modf(zfraction*(
nz-1), &zdindex));
168 int xindex =
static_cast<int>(std::floor(xdindex));
169 int yindex =
static_cast<int>(std::floor(ydindex));
170 int zindex =
static_cast<int>(std::floor(zdindex));
173 if ((xindex < 0) || (xindex >=
nx - 1) ||
174 (yindex < 0) || (yindex >=
ny - 1) ||
175 (zindex < 0) || (zindex >=
nz - 1))
184 #ifdef DEBUG_INTERPOLATING_FIELD
185 G4cout <<
"Local x,y,z: " << xlocal <<
" " << ylocal <<
" " << zlocal <<
G4endl;
186 G4cout <<
"Index x,y,z: " << xindex <<
" " << yindex <<
" " << zindex <<
G4endl;
187 double valx0z0, mulx0z0, valx1z0, mulx1z0;
188 double valx0z1, mulx0z1, valx1z1, mulx1z1;
189 valx0z0= table[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
190 valx1z0= table[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
191 valx0z1= table[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
192 valx1z1= table[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
197 xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
198 xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
199 xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
200 xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
201 xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
202 xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
203 xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
204 xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
207 yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
208 yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
209 yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
210 yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
211 yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
212 yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
213 yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
214 yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
217 zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
218 zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
219 zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
220 zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
221 zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
222 zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
223 zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
224 zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;