Commit in GeomConverter/src/org/lcsim/geometry/field on MAIN | |||
FieldMap3D.java | -299 | 1.2 removed |
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; - } -};
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