GeomConverter/src/org/lcsim/geometry/field
diff -u -r1.1 -r1.2
--- FieldMap3D.java 29 Aug 2013 22:37:57 -0000 1.1
+++ FieldMap3D.java 3 Sep 2013 20:20:03 -0000 1.2
@@ -197,88 +197,103 @@
void getField(double x, double y, double z, BasicHep3Vector field)
{
- /*
- double x = point[0];
- double y = point[1];
- double z = point[2] + _zOffset;
+ //double x = point[0];
+ //double y = point[1];
+ //double z = point[2] + _zOffset;
+ z = z + _zOffset;
// Check that the point is within the defined region
if ( x>=_minx && x<=_maxx &&
- y>=_miny && y<=_maxy &&
- z>=_minz && z<=_maxz ) {
+ 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
- cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << endl;
- cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << endl;
- double valx0z0, mulx0z0, valx1z0, mulx1z0;
- double valx0z1, mulx0z1, valx1z1, mulx1z1;
- valx0z0= table[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
- valx1z0= table[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
- valx0z1= table[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
- valx1z1= table[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
- #endif
-
- // Full 3-dimensional version
- Bfield[0] =
- _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 ;
- Bfield[1] =
- _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 ;
- Bfield[2] =
- _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 ;
-
- } else {
- Bfield[0] = 0.0;
- Bfield[1] = 0.0;
- Bfield[2] = 0.0;
- }
- */
+ // 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;
+ double xlocal, ylocal, zlocal;
+ xdindex = ydindex = zdindex = 0;
+ double results[] = new double[2];
+
+ // Position of the point within the cuboid defined by the
+ // nearest surrounding tabulated points
+ modf(xfraction*(_nx-1), results);
+ xdindex = results[0];
+ xlocal = results[1];
+
+ modf(yfraction*(_ny-1), results);
+ ydindex = results[0];
+ ylocal = results[1];
+
+ modf(zfraction*(_nz-1), results);
+ zdindex = results[0];
+ zlocal = results[1];
+
+ // The indices of the nearest tabulated point whose coordinates
+ // are all less than those of the given point
+ int xindex = (int)xdindex;
+ int yindex = (int)ydindex;
+ int zindex = (int)zdindex;
+
+ // cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << endl;
+ // cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << endl;
+ // double valx0z0, mulx0z0, valx1z0, mulx1z0;
+ // double valx0z1, mulx0z1, valx1z1, mulx1z1;
+ // valx0z0= table[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
+ // valx1z0= table[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
+ // valx0z1= table[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
+ // valx1z1= table[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
+
+ // Full 3-dimensional version
+ double bx, by, bz;
+ bx =
+ _xField.get(xindex).get(yindex).get(zindex) * (1-xlocal) * (1-ylocal) * (1-zlocal) +
+ _xField.get(xindex).get(yindex).get(zindex+1) * (1-xlocal) * (1-ylocal) * zlocal +
+ _xField.get(xindex).get(yindex+1).get(zindex) * (1-xlocal) * ylocal * (1-zlocal) +
+ _xField.get(xindex).get(yindex+1).get(zindex+1) * (1-xlocal) * ylocal * zlocal +
+ _xField.get(xindex+1).get(yindex).get(zindex) * xlocal * (1-ylocal) * (1-zlocal) +
+ _xField.get(xindex+1).get(yindex).get(zindex+1) * xlocal * (1-ylocal) * zlocal +
+ _xField.get(xindex+1).get(yindex+1).get(zindex) * xlocal * ylocal * (1-zlocal) +
+ _xField.get(xindex+1).get(yindex+1).get(zindex+1) * xlocal * ylocal * zlocal;
+ by =
+ _yField.get(xindex).get(yindex).get(zindex) * (1-xlocal) * (1-ylocal) * (1-zlocal) +
+ _yField.get(xindex).get(yindex).get(zindex+1) * (1-xlocal) * (1-ylocal) * zlocal +
+ _yField.get(xindex).get(yindex+1).get(zindex) * (1-xlocal) * ylocal * (1-zlocal) +
+ _yField.get(xindex).get(yindex+1).get(zindex+1)* (1-xlocal) * ylocal * zlocal +
+ _yField.get(xindex+1).get(yindex).get(zindex) * xlocal * (1-ylocal) * (1-zlocal) +
+ _yField.get(xindex+1).get(yindex).get(zindex+1) * xlocal * (1-ylocal) * zlocal +
+ _yField.get(xindex+1).get(yindex+1).get(zindex) * xlocal * ylocal * (1-zlocal) +
+ _yField.get(xindex+1).get(yindex+1).get(zindex+1) * xlocal * ylocal * zlocal ;
+ bz =
+ _zField.get(xindex).get(yindex).get(zindex) * (1-xlocal) * (1-ylocal) * (1-zlocal) +
+ _zField.get(xindex).get(yindex).get(zindex+1) * (1-xlocal) * (1-ylocal) * zlocal +
+ _zField.get(xindex).get(yindex+1).get(zindex) * (1-xlocal) * ylocal * (1-zlocal) +
+ _zField.get(xindex).get(yindex+1).get(zindex+1) * (1-xlocal) * ylocal * zlocal +
+ _zField.get(xindex+1).get(yindex).get(zindex) * xlocal * (1-ylocal) * (1-zlocal) +
+ _zField.get(xindex+1).get(yindex).get(zindex+1) * xlocal * (1-ylocal) * zlocal +
+ _zField.get(xindex+1).get(yindex+1).get(zindex) * xlocal * ylocal * (1-zlocal) +
+ _zField.get(xindex+1).get(yindex+1).get(zindex+1) * xlocal * ylocal * zlocal ;
+
+ field.setV(bx, by, bz);
+
+ } else {
+ field.setV(0., 0., 0.);
+ }
+ }
+
+ static private final void modf(double value, double results[])
+ {
+ int whole = (int)value;
+ double frac = value - ((float)whole);
+ results[0] = whole;
+ results[1] = frac;
}
};