GeomConverter/src/org/lcsim/geometry/field
diff -N FieldMap3D.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FieldMap3D.java 29 Aug 2013 22:37:57 -0000 1.1
@@ -0,0 +1,284 @@
+package org.lcsim.geometry.field;
+
+import hep.physics.vec.BasicHep3Vector;
+
+import java.io.BufferedReader;
+import java.io.FileInputStream;
+import java.io.IOException;
+import java.io.InputStream;
+import java.io.InputStreamReader;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.StringTokenizer;
+
+import org.jdom.Element;
+
+class FieldMap3D extends AbstractFieldMap
+{
+ // Storage space for the table
+ List<List<List<Double>>> _xField;
+ List<List<List<Double>>> _yField;
+ List<List<List<Double>>> _zField;
+
+ // The dimensions of the table
+ int _nx, _ny, _nz;
+
+ // The physical limits of the defined region
+ double _minx, _maxx, _miny, _maxy, _minz, _maxz;
+
+ // The physical extent of the defined region
+ double _dx, _dy, _dz;
+
+ // Offsets if field map is not in global coordinates
+ double _xOffset, _yOffset, _zOffset;
+
+ boolean _invertX, _invertY, _invertZ;
+
+ String _filename;
+
+ FieldMap3D(Element node)
+ {
+ super(node);
+ load(node.getAttributeValue("filename"));
+ }
+
+ private final void load(String filename)
+ {
+ System.out.println("-----------------------------------------------------------");
+ System.out.println("Magnetic field");
+ System.out.println("-----------------------------------------------------------");
+
+ System.out.println("Reading the field grid from " + filename + " ... ");
+
+ InputStream fis;
+ BufferedReader br;
+ String line;
+
+ try {
+ fis = new FileInputStream(filename);
+ br = new BufferedReader(new InputStreamReader(fis));
+ line = br.readLine();
+
+ // Read table dimensions.
+ readTableDimensions(line);
+
+ // Reserve space in table for field map data.
+ createTableBuffer();
+
+ // Ignore 7 lines of input.
+ for (int i=0; i<=7; i++)
+ {
+ line = br.readLine();
+ }
+
+ // Read the field data.
+ readFieldData(br);
+
+ br.close();
+
+ } catch (IOException x) {
+ throw new RuntimeException(x);
+ }
+ }
+
+ private final void readFieldData(BufferedReader br) {
+ System.out.println("readFieldData");
+
+ double xval, yval, zval, bx, by, bz;
+ xval = yval = zval = bx = by = bz = 0;
+ int ix, iy, iz;
+ ix = iy = iz = 0;
+ for (ix=0; ix<_nx; ix++) {
+ for (iy=0; iy<_ny; iy++) {
+ for (iz=0; iz<_nz; iz++) {
+ String line = null;
+ try {
+ line = br.readLine();
+ } catch (IOException x) {
+ throw new RuntimeException(x);
+ }
+ double fieldValues[] = new double[6];
+ StringTokenizer st = new StringTokenizer(" ", line);
+ int i = 0;
+ while (st.hasMoreElements()) {
+ fieldValues[i] = Double.parseDouble(st.nextToken());
+ i++;
+ }
+ xval = fieldValues[0];
+ yval = fieldValues[1];
+ zval = fieldValues[2];
+ bx = fieldValues[3];
+ by = fieldValues[4];
+ bz = fieldValues[5];
+
+ if ( ix==0 && iy==0 && iz==0 ) {
+ _minx = xval;
+ _miny = yval;
+ _minz = zval;
+ }
+ _xField.get(ix).get(iy).add(iz, bx);
+ _yField.get(ix).get(iy).add(iz, by);
+ _zField.get(ix).get(iy).add(iz, bz);
+ }
+ }
+ }
+ _maxx = xval;
+ _maxy = yval;
+ _maxz = zval;
+
+ // Should really check that the limits are not the wrong way around.
+ double temp = 0;
+ if (_maxx < _minx) {
+ temp = _maxx;
+ _maxx = _minx;
+ _minx = temp;
+ _invertX = true;
+ }
+
+ if (_maxy < _miny) {
+ temp = _maxy;
+ _maxy = _miny;
+ _miny = temp;
+ _invertY = true;
+ }
+
+ if (_maxz < _minz) {
+ temp = _maxz;
+ _maxz = _minz;
+ _minz = temp;
+ _invertZ = true;
+ }
+
+ // Set deltas from field values.
+ _dx = _maxx - _minx;
+ _dy = _maxy - _miny;
+ _dz = _maxz - _minz;
+
+ System.out.println("min values x, y, z: " + _minx + " " + _miny + " " + _minz);
+ System.out.println("max values x, y, z: " + _maxx + " " + _maxy + " " + _maxz);
+ System.out.println("offset x, y, z: " + _xOffset + " " + _yOffset + " " + _zOffset);
+ System.out.println("delta x, y, z: " + _dx + " " + _dy + " " + _dz);
+
+ System.out.println("done reading field data");
+ }
+
+ private final void readTableDimensions(String line) {
+ int gridSizes[] = new int[3];
+ StringTokenizer st = new StringTokenizer(" ", line);
+ int i = 0;
+ while (st.hasMoreElements()) {
+ gridSizes[i] = Integer.parseInt(st.nextToken());
+ ++i;
+ }
+ _nx = gridSizes[0];
+ _ny = gridSizes[1];
+ _nz = gridSizes[2];
+
+ System.out.println(" [ Number of values x,y,z: " + _nx + " " + _ny + " " + _nz + " ] ");
+ }
+
+ private final void createTableBuffer() {
+ _xField = new ArrayList<List<List<Double>>>(_nx);
+ _yField = new ArrayList<List<List<Double>>>(_nx);
+ _zField = new ArrayList<List<List<Double>>>(_nx);
+ int ix, iy, iz;
+ ix = iy = iz = 0;
+ for (ix=0; ix<_nx; ix++) {
+ _xField.add(ix, new ArrayList<List<Double>>(_ny));
+ _xField.add(iy, new ArrayList<List<Double>>(_ny));
+ _xField.add(iz, new ArrayList<List<Double>>(_ny));
+ for (iy=0; iy<_ny; iy++) {
+ _xField.get(ix).add(iy, new ArrayList<Double>(_nz));
+ _yField.get(ix).add(iy, new ArrayList<Double>(_nz));
+ _zField.get(ix).add(iy, new ArrayList<Double>(_nz));
+ }
+ }
+ }
+
+ void getField(double x, double y, double z, BasicHep3Vector field)
+ {
+ /*
+ double x = point[0];
+ double y = point[1];
+ double z = point[2] + _zOffset;
+
+ // 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
+ 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;
+ }
+ */
+ }
+};