Hi Mike,
I am implementing a 3D interpolation function which I have copied from the purging magnet example. In addition, I have implemented a messenger function /field/getField so I can determine the field at any given point.
At the bottom of this message is the code of the interpolation function (if your´re interested !). I wrote some code so that I could have a coarser field further away from the tip, but this adds a complication, so I have returned to the simpler code using just one scale for the interpolation.
I have also attached a screen shot of the trajectories. The electrons are launched just above the surface with an energy of 1 eV. Some electrons feel the field and bend back to the surface as they should, but others just fly away as if they do not feel the field at all.
Best Wishes
Chris
void PurgMagTabulatedEField3D::GetFieldValue(const double point[4],
double *BEfield ) const
{
BEfield[0] = 0.0;
BEfield[1] = 0.0;
BEfield[2] = 0.0;
double x = point[0];
double y = point[1];
double z = point[2] + fZoffset;
#ifdef DEBUG_INTERPOLATING_FIELD
G4cout << "x, y, z: " << x << " " << y << " " << z << endl;
G4cout << "minx, miny, minz: " << minx << " " << miny << " " << minz << endl;
G4cout << "maxx, maxy, maxz: " << maxx << " " << maxy << " " << maxz << endl;
#endif
// Check that the point is within the defined region
if ( x>=minx && x<=maxx &&
y>=miny && y<=maxy &&
z>=minz && z<=maxz ) {
// Position of given point within region, normalized to the range
// [0,1]
double xfraction = (x - minx) / dx;
double yfraction = (y - miny) / dy;
double zfraction = (z - minz) / dz;
if (invertX) { xfraction = 1 - xfraction;}
if (invertY) { yfraction = 1 - yfraction;}
if (invertZ) { zfraction = 1 - zfraction;}
// Need addresses of these to pass to modf below.
// modf uses its second argument as an OUTPUT argument.
double xdindex, ydindex, zdindex;
// Position of the point within the cuboid defined by the
// nearest surrounding tabulated points
double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));
// The indices of the nearest tabulated point whose coordinates
// are all less than those of the given point
int xindex = static_cast<int>(xdindex);
int yindex = static_cast<int>(ydindex);
int zindex = static_cast<int>(zdindex);
#ifdef DEBUG_INTERPOLATING_FIELD
G4cout << "Fraction x,y,z: " << xfraction << " " << yfraction << " " << zfraction << endl;
G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << endl;
G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << endl;
double valx0z0, mulx0z0, valx1z0, mulx1z0;
double valx0z1, mulx0z1, valx1z1, mulx1z1;
valx0z0= xField[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
valx1z0= xField[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
valx0z1= xField[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
valx1z1= xField[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
#endif
// Full 3-dimensional version
BEfield[3] =
xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
BEfield[4] =
yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
BEfield[5] =
zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
#ifdef DEBUG_INTERPOLATING_FIELD
G4cout << "minx: " << minx << " maxx: " << maxx << " x: " << x << " BEfield[3]: " << BEfield[3] << G4endl;
G4cout << "miny: " << miny << " maxy: " << maxy << " y: " << y << " BEfield[4]: " << BEfield[4] << G4endl;
G4cout << "minz: " << minz << " maxz: " << maxz << " z: " << z << " BEfield[5]: " << BEfield[5] << G4endl;
#endif
} else {
BEfield[3] = 0.0;
BEfield[4] = 0.0;
BEfield[5] = 0.0;
#ifdef DEBUG_INTERPOLATING_FIELD
// G4cout << "minx: " << minx << " maxx: " << maxx << " x: " << x << G4endl;
// G4cout << "miny: " << miny << " maxy: " << maxy << " y: " << y << G4endl;
// G4cout << "minz: " << minz << " maxz: " << maxz << " z: " << z << G4endl;
#endif
}
}