Commit in GeomConverter/src/org/lcsim/geometry/field on MAIN | |||
FieldMap3D.java | +284 | added 1.1 |
work in progress on porting C++ 3D field class to Java
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; + } + */ + } +};
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