lcdd/src

G4RZFieldMap.cc 1.2 -> 1.3

diff -u -r1.2 -r1.3
--- G4RZFieldMap.cc	30 Aug 2005 21:22:07 -0000	1.2
+++ G4RZFieldMap.cc	31 Aug 2005 01:21:26 -0000	1.3
@@ -1,77 +1,237 @@

-// $Header: /cvs/lcd/lcdd/src/G4RZFieldMap.cc,v 1.2 2005/08/30 21:22:07 jeremy Exp $

+// $Header: /cvs/lcd/lcdd/src/G4RZFieldMap.cc,v 1.3 2005/08/31 01:21:26 jeremy Exp $

 #include "G4RZFieldMap.hh"

-G4RZFieldMap::G4RZFieldMap()
-  : arraySize(0),
-    array(0),
-    madeArray(false)
-{}

+#include "globals.hh"
+
+#include <cmath>
+
+G4RZFieldMap::G4RZFieldMap(double spanR, double gridSizeR, double spanZ, double gridSizeZ)
+  : m_numBinsR(0),
+    m_numBinsZ(0),
+    m_spanR(spanR),
+    m_gridSizeR(gridSizeR),
+    m_spanZ(spanZ),
+    m_gridSizeZ(gridSizeZ),
+    m_BrArray(0),
+    m_BzArray(0),
+    m_madeArray(false)
+{
+  computeNumBins();
+}

 
 G4RZFieldMap::~G4RZFieldMap()
 {
   /* delete the array */

-  if ( madeArray ) {
-    for (size_t i=0; i<arraySize; i++) {
-      delete[] array[i];

+  if ( m_madeArray ) {
+    for (size_t i=0; i<m_numBinsR; i++) {
+      delete[] m_BrArray[i];
+      delete[] m_BzArray[i];

-    delete[] array;

+    delete[] m_BrArray;
+    delete[] m_BzArray;
+  }
+}
+
+void G4RZFieldMap::computeNumBins()
+{
+  /* r bins */
+  double rem = fmod( (float)m_spanR, (float)m_gridSizeR);
+
+  if ( rem != 0 ) {
+    G4Exception("spanR is not evenly divisible by gridSizeR.");
+  }
+
+  m_numBinsR = (size_t)(m_spanR / m_gridSizeR);
+
+  if ( m_numBinsR == 0 ) {
+    G4Exception("numBinsR == 0");
+  }
+
+  /* z bins */
+  rem = fmod( (float)m_spanZ, (float)m_gridSizeZ);
+
+  if ( rem != 0 ) {
+    G4Exception("spanZ is not evenly disible by gridSizeZ.");
+  }
+
+  m_numBinsZ = (size_t)(m_spanZ / m_gridSizeZ);
+
+  if ( m_numBinsZ == 0 ) {
+    G4Exception("numBinsZ == 0");

   }
 }
 
 /*

- *  FIXME: Needs implementation!
- *  See Takashi's FORTRAN code at /afs/slac/www/accel/nlc/local/systems/beamdelivery/geant/SDNEW/solenoid.f

+ * Compute B field at given point using the arrays of Br and Bz field strengths.
+ *
+ * This function is based on Takashi's FORTRAN code at
+ *
+ * /afs/slac/www/accel/nlc/local/systems/beamdelivery/geant/SDNEW/solenoid.f
+ *

  */
 void G4RZFieldMap::GetFieldValue( const double Point[3], double *Bfield ) const
 {

-  /* make compiler happy */
-  if ( Point[0] );
-  if ( Bfield );

+  std::cout << std::endl << "G4RZFieldMap::GetFieldValue" << std::endl;
+  std::cout << "Point = " << Point[0] << ", " << Point[1] << ", " << Point[2] << std::endl;
+
+  double rcyl = sqrt( ( Point[0] * Point[0] ) + ( Point[1] * Point[1] ) );
+  double z = Point[2];
+
+  size_t binR = (int)floor( rcyl / m_gridSizeR );
+  binR += 1;
+
+  std::cout << "binR = " << binR << std::endl;
+
+  if ( binR >= m_numBinsR ) {
+    std::cerr << "OVERFLOW: binR > m_numBinsR - 1 ... RETURNING" << std::endl;
+    return;
+  }
+
+  size_t binZ = (size_t) floor( ( fabs(z) / m_gridSizeZ ) );
+  binZ += 1;
+
+  std::cout << "binZ = " << binZ << std::endl;
+
+  if ( binZ >= m_numBinsZ ) {
+    std::cerr << "OVERFLOW: binZ > m_numBinsZ - 1 ... SETTING binZ = m_numBinZ" << std::endl;
+    binZ = m_numBinsZ;
+  }
+
+  double brmin = m_BrArray[binR][binZ];
+  double bzmin = m_BzArray[binR][binZ];
+
+  double remR = (double) ( rcyl - ( (binR - 1) * m_gridSizeR ) );
+  double brdz = (m_BrArray[binZ+1][binR] - brmin) / m_gridSizeZ;
+  double brdr = (m_BrArray[binZ][binR+1] - brmin) / m_gridSizeR;
+
+  double remZ = (double) fabs(z) - ( (binZ - 1) * m_gridSizeZ );
+  double bzdz = (m_BzArray[binZ+1][binR] - bzmin) / m_gridSizeZ;
+  double bzdr = (m_BzArray[binZ][binR+1] - bzmin) / m_gridSizeR;
+
+  double bz = bzmin + bzdz * remZ + bzdr * remR;
+  double br = brmin + brdz * remZ + brdr * remR;
+
+  double theta = atan2(Point[1], Point[0]);
+  double bx = br * cos(theta);
+  double by = br * sin(theta);
+
+  Bfield[0] = bx;
+  Bfield[1] = by;
+  Bfield[2] = bz;
+
+  // DEBUG: prints
+  std::cout << "brmin = " << brmin << std::endl;
+  std::cout << "bzmin = " << bzmin << std::endl;
+  std::cout << "remR = " << remR << std::endl;
+  std::cout << "brdz = " << brdz << std::endl;
+  std::cout << "brdr = " << brdr << std::endl;
+  std::cout << "bzdz = " << bzdz << std::endl;
+  std::cout << "bzdr = " << bzdr << std::endl;
+  std::cout << "bz = " << bz << std::endl;
+  std::cout << "br = " << br << std::endl;
+  std::cout << "Bfield = " << Bfield[0] << ", " << Bfield[1] << ", " << Bfield[2] << std::endl;

 }
 
 void G4RZFieldMap::addRZB(double r, double z, double Br, double Bz)
 {
   //  std::cout << "adding --> r z Br Bz = " << r << " " << z << " " << Br << " " << Bz << std::endl;

-  if ( !madeArray ) {
-    rows.push_back(G4RZB(r, z, Br, Bz));

+  if ( !m_madeArray ) {
+    m_rows.push_back(G4RZB(r, z, Br, Bz));

   }
   else {

-    std::cerr << "ERROR: Additional rows not allowed once G4RZFieldMap::makeArray() has been called." << std::endl;

+    std::cerr << "ERROR: Additional m_rows not allowed once G4RZFieldMap::makeArray() has been called." << std::endl;

   }
 }

-void G4RZFieldMap::makeArray(bool clearRows)

+void G4RZFieldMap::makeArray(bool clearM_Rows)

-  if ( !madeArray ) {

+  if ( !m_madeArray ) {

-    /* allocate array of [arraySize][4] */
-    arraySize = rows.size();
-    array = new double*[arraySize];

+    /* Make two 2D arrays of numBinsR x numBinsZ */
+    m_BrArray = new double*[m_numBinsR];
+    m_BzArray = new double*[m_numBinsR];
+
+    for (size_t i=0; i<m_numBinsR; i++) {
+      m_BrArray[i] = new double[m_numBinsZ];
+      m_BzArray[i] = new double[m_numBinsZ];
+    }

-    for (size_t i=0; i<arraySize; i++) {
-      array[i] = new double[4];

+    /* Init arrays to all 0's */
+    for ( size_t i=0; i<m_numBinsR; i++) {
+      for ( size_t j=0; j<m_numBinsZ; j++) {
+	m_BrArray[i][j]=0;
+	m_BzArray[i][j]=0;
+      }

-    /* fill array with values from object vector */
-    int cntr = 0;
-    for ( RowsType::const_iterator iter = rows.begin();
-	  iter != rows.end();

+    /* Fill arrays with values from object vector. */
+    size_t binR = 0;
+    size_t binZ = 0;
+    for ( RowsType::const_iterator iter = m_rows.begin();
+	  iter != m_rows.end();

 	  iter++ ) {

-      array[cntr][0] = (*iter).r();
-      array[cntr][1] = (*iter).z();
-      array[cntr][2] = (*iter).Br();
-      array[cntr][3] = (*iter).Bz();
-      ++cntr;

+
+      /* Compute the r bin. */
+      binR = (size_t) floor( (*iter).r() / m_gridSizeR );
+
+      double r = (*iter).r();
+      if ( r == m_spanR ) {
+	binR -= 1;
+      }
+      else if ( r > m_spanR ) {
+	G4Exception("r value is > spanR");
+      }
+
+      std::cout << std::endl;
+      std::cout << "r = " << (*iter).r() << std::endl;
+      std::cout << "binR = " << binR << std::endl;
+
+      /* r overflow? */
+      if ( binR > ( m_numBinsR - 1 ) ) {
+	G4Exception("OVERFLOW: R bin > numBinsR - 1");
+      }
+
+      /*
+       * Compute z bin, normalized to positive value.
+       * Forces z symmetry w.r.t. the +/- z halves.
+      */
+      double z = fabs((*iter).z());
+      binZ = (size_t) floor( z / m_gridSizeZ );
+
+      /* If Z value is equal to spanZ, set to last z bin. */
+      if (z == m_spanZ) {
+	binZ -= 1;
+      }
+      else if ( z > m_spanZ ) {
+	G4Exception("z value is > spanZ");
+      }
+
+      std::cout << "z = " << z << std::endl;
+      std::cout << "binZ = " << binZ << std::endl;
+
+      /* z overflow? */
+      if ( binZ > ( m_numBinsZ - 1) ) {
+	G4Exception("OVERFLOW: Z bin > numBinsZ - 1");
+      }
+
+      /* Assign Bz and Br array values. */
+      std::cout << "setting r,z bins = " << binR << ", " << binZ << std::endl;
+
+      m_BrArray[binR][binZ] = (*iter).Bz();
+      m_BzArray[binR][binZ] = (*iter).Br();
+
+      std::cout << "br[r][z] = " << m_BrArray[binR][binZ] << std::endl;
+      std::cout << "br[r][z] = " << m_BzArray[binR][binZ] << std::endl;

-    /* clear object vector rows */
-    if ( clearRows ) {
-      rows.clear();

+    /* clear object vector m_rows */
+    if ( clearM_Rows ) {
+      m_rows.clear();

-    madeArray = true;

+    m_madeArray = true;

-  else {

+    else {

     std::cerr << "WARNING: Ignoring additional call to G4RZFieldMap::makeArray()." << std::endl;
   }
 }

Commit in `lcdd` on MAIN
`include/G4RZFieldMap.hh`	+21	-7	1.2 -> 1.3
`/RZFieldMapType.hh`	+50	-2	1.1 -> 1.2
`/rz_field_map.hh`	+3	-3	1.1 -> 1.2
`schemas/lcdd/1.0/lcdd_fields.xsd`	+4		1.2 -> 1.3
`src/G4RZFieldMap.cc`	+199	-39	1.2 -> 1.3
`/rz_field_mapProcess.cc`	+7	-1	1.1 -> 1.2
`/rz_field_mapSubscriber.cc`	+20	-2	1.2 -> 1.3
	+304	-54