Commit in GeomConverter/src/org/lcsim/geometry/field on MAIN | |||
FieldMap3D.java | +94 | -79 | 1.1 -> 1.2 |
implementation of 3D field map; needs to be merged with Norman's class
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;
} };
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1