GeomConverter/src/org/lcsim/geometry/field
diff -N FieldMap3D.java
--- FieldMap3D.java 3 Sep 2013 20:20:03 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,299 +0,0 @@
-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;
- 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 ) {
-
- // 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;
- }
-};