36 :fXoffset(xOffset),invertX(false),invertY(false),invertZ(false)
40 double lenUnit=
meter;
41 double fieldUnit=
tesla;
42 G4cout <<
"\n-----------------------------------------------------------"
43 <<
"\n Magnetic field"
44 <<
"\n-----------------------------------------------------------";
47 G4cout <<
"\n ---> " "Reading the field grid from " << filename <<
" ... " << endl;
50 ifstream file( filename );
54 file.getline(buffer,256);
59 G4cout <<
" [ Number of values x,y,z: "
60 <<
nx <<
" " <<
ny <<
" " << nz <<
" ] "
68 for (ix=0; ix<
nx; ix++) {
72 for (iy=0; iy<
ny; iy++) {
86 for (ix=0; ix<
nx; ix++) {
87 for (iy=0; iy<
ny; iy++) {
88 for (iz=0; iz<
nz; iz++) {
89 file >> xval >> yval >> zval >> bx >> by >> bz ;
90 if ( ix==0 && iy==0 && iz==0 ) {
91 minx = xval * lenUnit;
92 miny = yval * lenUnit;
93 minz = zval * lenUnit;
105 maxx = xval * lenUnit;
106 maxy = yval * lenUnit;
107 maxz = zval * lenUnit;
109 G4cout <<
"\n ---> ... done reading " << endl;
112 G4cout <<
" ---> assumed the order: x, y, z, Bx, By, Bz "
113 <<
"\n ---> Min values x,y,z: "
115 <<
"\n ---> Max values x,y,z: "
117 <<
"\n ---> The field will be offset by " << xOffset/
cm <<
" cm " << endl;
123 G4cout <<
"\nAfter reordering if neccesary"
124 <<
"\n ---> Min values x,y,z: "
126 <<
" \n ---> Max values x,y,z: "
132 G4cout <<
"\n ---> Dif values x,y,z (range): "
134 <<
"\n-----------------------------------------------------------" << endl;
138 double *Bfield )
const
150 double xfraction = (x -
minx) /
dx;
151 double yfraction = (y -
miny) /
dy;
152 double zfraction = (z -
minz) /
dz;
154 if (
invertX) { xfraction = 1 - xfraction;}
155 if (
invertY) { yfraction = 1 - yfraction;}
156 if (
invertZ) { zfraction = 1 - zfraction;}
160 double xdindex, ydindex, zdindex;
164 double xlocal = ( std::modf(xfraction*(
nx-1), &xdindex));
165 double ylocal = ( std::modf(yfraction*(
ny-1), &ydindex));
166 double zlocal = ( std::modf(zfraction*(
nz-1), &zdindex));
170 int xindex =
static_cast<int>(xdindex);
171 int yindex =
static_cast<int>(ydindex);
172 int zindex =
static_cast<int>(zdindex);
175 #ifdef DEBUG_INTERPOLATING_FIELD
176 G4cout <<
"Local x,y,z: " << xlocal <<
" " << ylocal <<
" " << zlocal << endl;
177 G4cout <<
"Index x,y,z: " << xindex <<
" " << yindex <<
" " << zindex << endl;
178 double valx0z0, mulx0z0, valx1z0, mulx1z0;
179 double valx0z1, mulx0z1, valx1z1, mulx1z1;
180 valx0z0= table[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
181 valx1z0= table[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
182 valx0z1= table[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
183 valx1z1= table[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
187 xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
188 xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
189 xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
190 xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
191 xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
192 xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
193 xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
194 xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
197 yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
198 yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
199 yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
200 yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
201 yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
202 yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
203 yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
204 yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
207 zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
208 zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
209 zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
210 zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
211 zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
212 zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
213 zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
214 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
const G4double x[NPOINTSGL]
static const double tesla