41 :fXoffset(xOffset),invertX(false),invertY(false),invertZ(false)
45 double lenUnit=
meter;
46 double fieldUnit=
tesla;
47 G4cout <<
"\n-----------------------------------------------------------"
48 <<
"\n Magnetic field"
49 <<
"\n-----------------------------------------------------------";
52 G4cout <<
"\n ---> " "Reading the field grid from " << filename <<
" ... " << endl;
55 ifstream file( filename );
59 file.getline(buffer,256);
64 G4cout <<
" [ Number of values x,y,z: "
65 <<
nx <<
" " <<
ny <<
" " << nz <<
" ] "
73 for (ix=0; ix<
nx; ix++) {
77 for (iy=0; iy<
ny; iy++) {
91 for (ix=0; ix<
nx; ix++) {
92 for (iy=0; iy<
ny; iy++) {
93 for (iz=0; iz<
nz; iz++) {
94 file >> xval >> yval >> zval >> bx >> by >> bz ;
95 if ( ix==0 && iy==0 && iz==0 ) {
96 minx = xval * lenUnit;
97 miny = yval * lenUnit;
98 minz = zval * lenUnit;
100 xField[ix][iy][
iz] = bx * fieldUnit;
101 yField[ix][iy][
iz] = by * fieldUnit;
102 zField[ix][iy][
iz] = bz * fieldUnit;
110 maxx = xval * lenUnit;
111 maxy = yval * lenUnit;
112 maxz = zval * lenUnit;
114 G4cout <<
"\n ---> ... done reading " << endl;
117 G4cout <<
" ---> assumed the order: x, y, z, Bx, By, Bz "
118 <<
"\n ---> Min values x,y,z: "
120 <<
"\n ---> Max values x,y,z: "
122 <<
"\n ---> The field will be offset by " << xOffset/
cm <<
" cm " << endl;
128 G4cout <<
"\nAfter reordering if neccesary"
129 <<
"\n ---> Min values x,y,z: "
131 <<
" \n ---> Max values x,y,z: "
137 G4cout <<
"\n ---> Dif values x,y,z (range): "
139 <<
"\n-----------------------------------------------------------" << endl;
143 double *Bfield )
const
155 double xfraction = (x -
minx) /
dx;
156 double yfraction = (y -
miny) /
dy;
157 double zfraction = (z -
minz) /
dz;
159 if (
invertX) { xfraction = 1 - xfraction;}
160 if (
invertY) { yfraction = 1 - yfraction;}
161 if (
invertZ) { zfraction = 1 - zfraction;}
165 double xdindex, ydindex, zdindex;
169 double xlocal = ( std::modf(xfraction*(
nx-1), &xdindex));
170 double ylocal = ( std::modf(yfraction*(
ny-1), &ydindex));
171 double zlocal = ( std::modf(zfraction*(
nz-1), &zdindex));
175 int xindex =
static_cast<int>(xdindex);
176 int yindex =
static_cast<int>(ydindex);
177 int zindex =
static_cast<int>(zdindex);
180 #ifdef DEBUG_INTERPOLATING_FIELD
181 G4cout <<
"Local x,y,z: " << xlocal <<
" " << ylocal <<
" " << zlocal << endl;
182 G4cout <<
"Index x,y,z: " << xindex <<
" " << yindex <<
" " << zindex << endl;
183 double valx0z0, mulx0z0, valx1z0, mulx1z0;
184 double valx0z1, mulx0z1, valx1z1, mulx1z1;
185 valx0z0= table[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
186 valx1z0= table[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
187 valx0z1= table[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
188 valx1z1= table[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
192 xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
193 xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
194 xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
195 xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
196 xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
197 xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
198 xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
199 xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
202 yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
203 yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
204 yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
205 yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
206 yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
207 yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
208 yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
209 yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
212 zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
213 zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
214 zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
215 zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
216 zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
217 zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
218 zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
219 zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
vector< vector< vector< double > > > xField
HadrontherapyMagneticField3D(const char *filename, double xOffset)
vector< vector< vector< double > > > zField
#define G4MUTEX_INITIALIZER
void GetFieldValue(const double Point[4], double *Bfield) const
G4GLOB_DLL std::ostream G4cout
vector< vector< vector< double > > > yField
static const double meter
static const double tesla