#include <BLFieldMap.hh>
Initially this class simply reads an input file to define the map. Eventually it will also generate a map automatically from a BLCoil, and from an arbitrary collection of magnets (to improve tracking efficiency for a large number of overlapping solenoids).
Field components are interpolated in the map; the values used are: Bx = (interpolated value)*normB*(element_current/current_param) Ex = (interpolated value)*normE*(element_gradient/gradient_param) [other components are similar] where element_current and element_gradient are from the element definition command in the input file, and normB, normE, curren_param, and gradient_param come from the input file. The presence of both norm and current/gradient in the input file is to accommodate diverse sources of input files and units.
Outside the map the fields are zero, so if the map is truncated you should ensure that no particles are tracked outside the map in that region. This happens, for instance, for a beam solenoid truncated in radius to the inside of the coil -- particles which enter the coil will see a zero field, so the coil should be set to kill them; particles outside the map along z but inside the bore will also see a zero field, but that is usually OK as long as the map extends far enough along Z to include the nonzero field region.
Input File format: Blank lines, and lines beginning with # or * are comments. Lines beginning with * are printed to stdout. Units are mm for coordinates, Tesla for B, and MegaVolts/meter for E; use normB and normE if the data points use different units.
The input file starts with a set of commands to define the parameters of the map, followed by blocks of lines containing the values of the field components. The field component names depend on the type of map (grid: Bx,By,Bz,Ex,Ey,Ez; cylinder: Br,Bz,Er,Ez). Each command has a specific list of arguments to define parameters of the map.
BEWARE: the parsing is not exhaustive. For instance, invalid arguments are silently ignored (which means you must verify the spelling and capitalization of argument names). Correct inputs will yield correct results, but invalid inputs may not be detected and may yield seemingly-correct but unintended results.
The first command is usually a param command, which has the following arguments: maxline The maximum number of characters per line (default=1023) current The current corresponding to this map (default=1.0) gradient The gradient corresponding to the map (default=1.0) normE A normalization factor for E components (default=1.0) normB A normalization factor for B components (default=1.0)
Two types of maps are implemented: grid and cylinder.
grid maps are a 3-D grid, with each block of data being a single X-Y plane; within a block the lines are Y and the columns of each line are values along X. The grid command has the following arguments: X0 The X value for the first value in each line Y0 The Y value for the first line in each block Z0 The Z value for the first block of each field component nX The number of columns per line nY The number of lines per block nZ The number of blocks per field component dX The X increment between values in each line dY The Y increment between lines dZ The Z increment between blocks tolerance The tolerance for pointwise data (default=0.01 mm) After the grid command, the following optional commands can be given: extendX flip=... extendY flip=... extendZ flip=... These commands permit a half-map to be extended to the full map around X=0, Y=0, or Z=0 respectively. The optional flip argument is a comma-separated list of field components whose signs will be inverted for negative values of the coordinate. For example, "extendZ flip=Bx,Ex" means the map from Z=0 to Z=(nZ-1)*dZ is extended symmetrically around Z=0 to negative Z values, flipping the signs of Bx and Ex when Z<0. This could be followed by "extendX flip=Bx,Ex", and the field flips will be the products of both commands.
cylinder maps are a 2-D map with rotational symmetry around the Z axis. Each field component has a single block with lines being Z and the columns being R. The cylinder command has the following arguments: Z0 The Z value for the first line in each block nR The number of columns per line nZ The number of lines per block dR The R increment between colums dZ The Z increment between lines tolerance The tolerance for pointwise data (default=0.01 mm) After the cylinder command, the following optional commands can be given: extendZ flip=... This command behaves the same as for the grid map.
After the commands, each block consists of a line containing the name of the field component, followed by the lines of the block. The values within a line can be separated by whitespace or a ',' followed by optional whitespace. Field components that are not given are set to 0.0 everywhere. Missing values will be considered to be 0.0. For grid maps the first block is for Z=Z0, and successive blocks increment Z by dZ; the first line in a block is for Y=Y0 and the first column in each line is for X=X0. For cylinder maps, the first line in each block is for Z=Z0 and the first column in each line is for R=0.
Instead of the blocked input format, a pointwise data format can be used. This is introduced by a line containing the command "data", followed by the individual points of the map, one per line. for a grid field, each line contains values for X,Y,Z,Bx,By,Bz,Ex,Ey,Ez separated by either a comma and optional whitespece or by whitespace. for a cylinder field each line contains values for R,Z,Br,Bz,Er,Ez. The order of the points does not matter; omitted grid points will be 0.0, and for duplicates the last entry wins. If there is no E field, the Ex,Ey,Ez or Er,Ez entries should be omitted on every line. NOTE: every line's X,Y,Z or R,Z must be on a grid point as specified by the arguments to the grid or cylinder commands, to within the tolerance specified; if not, an error message is printed and the input line is ignored.
For time=dependent fields, the "time" command is used: time [period=12] period, if given, is in nanoseconds, and causes the interval [0,period) to be extended forever (before and after the values given). Because of the interpolation used, at least two points beyond the interval boundaries should be provided; there need not be a point at either boundary (but usually there are). Following the time command are lines containing 2 or 3 doubles: t B E where t is the time (nanoseconds), and B and E are factors for the fields. If E is omitted, the value for B is used. These values will be interpolated in time with a cubic spline that can handle either uniform or non-uniform spacing of points along t. The time command can come either before or after the cylinder or grid commands, but not within either of their sequences. Note a cubic spline is used to interpolate between points, and that can cause over/under-shoot near an abrupt change. Combined with period= this gives an excellent representation of sinewave/cosinewave.
Note that time dependence can currently only be specified via the time command in an input file (i.e. not programmable method exists).
Example block input file: * this is an example BLFieldMap input file, suitable for a solenoid # grid interval is 1 cm. # The region of validity is -390<=Z<=390 and 0<=R<=90 param normB=1.0 current=1.0 cylinder Z0=0.0 nR=10 nZ=40 dR=10.0 dZ=10.0 extendZ flip=Br Bz ... 40 lines of 10 values, Z=0 thru Z=390 Br ... 40 lines of 10 values, Z=0 thru Z=390 --EOF--
Example pointwise input file: * this is an example BLFieldMap input file, suitable for a solenoid # grid interval is 1 cm. # The region of validity is -390<=Z<=390 and 0<=R<=90 param normB=1.0 current=1.0 cylinder Z0=0.0 nR=10 nZ=40 dR=10.0 dZ=10.0 extendZ flip=Br data ... 400 lines of 4 values, giving R,Z,Br,Bz --EOF--
Public Member Functions | |
BLFieldMap () | |
default constructor. | |
virtual | ~BLFieldMap () |
destructor. | |
bool | readFile (G4String filename) |
readFile() reads a file to initialize the map. Returns true if OK, false on error. | |
bool | writeFile (G4String filename, G4String comment="") |
writeFile() writes the map to a file. Returns true if OK, false on error. | |
void | getFieldValue (const G4double local[4], G4double field[6], G4double current=1.0, G4double gradient=1.0) |
getFieldValue() gets the map's field to the field[] array. local[] must be LOCAL coordinates; local[3] (time) is ignored unless there is a time dependence specified. (elements using BLFieldMap perform the global->local conversion) | |
void | getFieldValueNoTimeNoScale (const G4double local[4], G4double field[6]) |
getFieldValueNoTimeNoScale() as above, except no time dependence is applies, and no scaling is performed. | |
void | getBoundingPoint (int i, G4double point[4]) |
getBoundingPoint() returns the i-th bounding point of the map. | |
bool | hasB () |
hasB() returns true if this map has a nonzero B field. | |
bool | hasE () |
hasE() returns true if this map has a nonzero E field. | |
bool | createGridMap (G4double X0, G4double Y0, G4double Z0, G4double dX, G4double dY, G4double dZ, int nX, int nY, int nZ, class G4ElectroMagneticField *field) |
createGridMap() will create a 3-d map from a pre-defined field. returns false if error. | |
bool | createCylinderMap (G4double Z0, G4double dR, G4double dZ, int nR, int nZ, class G4ElectroMagneticField *field) |
createCylinderMap() will create a 2-d map from a pre-defined cylindrically-symmetric field. NOTE: cylindrical symmetry is ASSUMED; the field in the Y=0,X>=0 half-plane is used (Bx->Br...). returns false if error. | |
bool | createTimeDependence (int n, G4double t[], G4double b[], G4double e[]=0, G4double period=-1.0) |
createTimeDependence() will apply the time dependence given. n is the # elements in the arrays; t[] is the time value for each point, b[] is the factor for B, and e[] is the factor for E. returns false if error. | |
bool | getTimeFactor (G4double t, G4double *b, G4double *e) |
getTimeFactor() returns the time factors for B and E at time t. returns false if error. | |
Private Attributes | |
G4int | maxline |
G4double | current |
G4double | gradient |
G4double | normB |
G4double | normE |
class FieldMapImpl * | impl |
class TimeImpl * | time |
Friends | |
class | FieldMapPlacement |
BLFieldMap::BLFieldMap | ( | ) |
BLFieldMap::~BLFieldMap | ( | ) | [virtual] |
bool BLFieldMap::readFile | ( | G4String | filename | ) |
readFile() reads a file to initialize the map. Returns true if OK, false on error.
References argDouble(), argInt(), InputFile::close(), current, InputFile::filename(), InputFile::getline(), InputFile::good(), gradient, FieldMapImpl::handleCommand(), impl, InputFile::linenumber(), maxline, normB, normE, BLCommand::parseArgs(), BLCommand::printError(), TimeImpl::readTime(), InputFile::setMaxline(), and time.
Referenced by BLCMDrfdevice::command(), BLCMDpillbox::command(), and BLCMDfieldmap::command().
00313 { 00314 InputFile in(filename,maxline); 00315 if(!in.good()) { 00316 BLCommand::printError("BLFieldMap Cannot open file '%s'", 00317 in.filename()); 00318 return false; 00319 } 00320 printf("BLFieldMap: reading file '%s'\n",in.filename()); 00321 00322 bool retval = true; 00323 00324 char *line; 00325 while((line=in.getline()) != 0) { 00326 BLArgumentVector argv; 00327 BLArgumentMap namedArgs; 00328 if(BLCommand::parseArgs(line,argv,namedArgs) < 0) 00329 goto invalid; 00330 if(argv[0] == "") continue; 00331 if(argv[0] == "param") { 00332 argInt(maxline,"maxline",namedArgs); 00333 argDouble(current,"current",namedArgs); 00334 argDouble(gradient,"gradient",namedArgs); 00335 argDouble(normB,"normB",namedArgs); 00336 argDouble(normE,"normE",namedArgs); 00337 in.setMaxline(maxline); 00338 } else if(argv[0] == "grid") { 00339 if(impl) goto invalid; 00340 impl = new GridImpl(argv,namedArgs); 00341 } else if(argv[0] == "cylinder") { 00342 if(impl) goto invalid; 00343 impl = new CylinderImpl(argv,namedArgs); 00344 } else if(argv[0] == "time") { 00345 if(time) goto invalid; 00346 time = TimeImpl::readTime(in,argv,namedArgs); 00347 if(!time) goto invalid; 00348 } else if(impl) { 00349 if(!impl->handleCommand(in,argv,namedArgs)) 00350 goto invalid; 00351 } else { 00352 invalid: BLCommand::printError("BLFieldMap file '%s' line %d invalid command '%s'\n", 00353 in.filename(),in.linenumber(), 00354 argv[0].c_str()); 00355 retval = false; 00356 break; 00357 } 00358 } 00359 in.close(); 00360 00361 return retval; 00362 }
bool BLFieldMap::writeFile | ( | G4String | filename, | |
G4String | comment = "" | |||
) |
writeFile() writes the map to a file. Returns true if OK, false on error.
References current, gradient, impl, maxline, normB, normE, and FieldMapImpl::writeFile().
Referenced by BLCMDprintfield::do_cylinder(), and BLCMDprintfield::do_grid().
00527 { 00528 if(!impl) return false; 00529 00530 FILE *f = fopen(filename.c_str(),"r"); 00531 if(f) { 00532 fclose(f); 00533 G4Exception("BLFieldMap","Output File Exists",FatalException, 00534 filename.c_str()); 00535 } 00536 00537 f = fopen(filename,"w"); 00538 if(!f) { 00539 fprintf(stderr,"BLFieldMap::writeFile CANNOT WRITE file '%s'\n", 00540 filename.c_str()); 00541 return false; 00542 } 00543 00544 fprintf(f,"# %s\n",comment.c_str()); 00545 fprintf(f,"param maxline=%d current=%g gradient=%g normB=%g normE=%g\n", 00546 maxline,current,gradient,normB,normE); 00547 00548 bool retval = impl->writeFile(f); 00549 00550 fclose(f); 00551 return retval; 00552 }
void BLFieldMap::getFieldValue | ( | const G4double | local[4], | |
G4double | field[6], | |||
G4double | current = 1.0 , |
|||
G4double | gradient = 1.0 | |||
) |
getFieldValue() gets the map's field to the field[] array. local[] must be LOCAL coordinates; local[3] (time) is ignored unless there is a time dependence specified. (elements using BLFieldMap perform the global->local conversion)
References current, TimeImpl::factorB(), TimeImpl::factorE(), FieldMapImpl::getFieldValue(), gradient, impl, normB, normE, and time.
Referenced by PillboxField::addFieldValue(), FieldMapPlacement::addFieldValue(), FieldExprPlacement::addFieldValue(), BLCMDfieldexpr::maxError(), and RFdeviceField::UserSteppingAction().
00283 { 00284 if(!impl) 00285 throw "BLFieldMap::getFieldValue called, no implementation"; 00286 00287 G4double thisField[6]; 00288 impl->getFieldValue(local,thisField); 00289 G4double timeB=1.0, timeE=1.0; 00290 if(time) { 00291 timeB = time->factorB(local[3]); 00292 timeE = time->factorE(local[3]); 00293 } 00294 field[0] = thisField[0] * normB * timeB * _current/current; 00295 field[1] = thisField[1] * normB * timeB * _current/current; 00296 field[2] = thisField[2] * normB * timeB * _current/current; 00297 field[3] = thisField[3] * normE * timeE * _gradient/gradient; 00298 field[4] = thisField[4] * normE * timeE * _gradient/gradient; 00299 field[5] = thisField[5] * normE * timeE * _gradient/gradient; 00300 }
void BLFieldMap::getFieldValueNoTimeNoScale | ( | const G4double | local[4], | |
G4double | field[6] | |||
) |
getFieldValueNoTimeNoScale() as above, except no time dependence is applies, and no scaling is performed.
References FieldMapImpl::getFieldValue(), and impl.
Referenced by RFdeviceField::addFieldValue().
00305 { 00306 if(!impl) 00307 throw "BLFieldMap::getFieldValueNoTimeNoScale called, no implementation"; 00308 00309 impl->getFieldValue(local,field); 00310 }
void BLFieldMap::getBoundingPoint | ( | int | i, | |
G4double | point[4] | |||
) |
getBoundingPoint() returns the i-th bounding point of the map.
References FieldMapImpl::getBoundingPoint(), and impl.
Referenced by FieldExprPlacement::FieldExprPlacement(), and FieldMapPlacement::FieldMapPlacement().
00366 { 00367 impl->getBoundingPoint(i,point); 00368 }
bool BLFieldMap::hasB | ( | ) |
hasB() returns true if this map has a nonzero B field.
References FieldMapImpl::hasB(), and impl.
Referenced by FieldMapPlacement::addFieldValue(), and FieldExprPlacement::addFieldValue().
bool BLFieldMap::hasE | ( | ) |
hasE() returns true if this map has a nonzero E field.
References FieldMapImpl::hasE(), and impl.
Referenced by FieldMapPlacement::addFieldValue(), and FieldExprPlacement::addFieldValue().
bool BLFieldMap::createGridMap | ( | G4double | X0, | |
G4double | Y0, | |||
G4double | Z0, | |||
G4double | dX, | |||
G4double | dY, | |||
G4double | dZ, | |||
int | nX, | |||
int | nY, | |||
int | nZ, | |||
class G4ElectroMagneticField * | field | |||
) |
createGridMap() will create a 3-d map from a pre-defined field. returns false if error.
References current, d2string(), gradient, i2string(), impl, maxline, normB, normE, and GridImpl::setField().
Referenced by BLCMDfieldexpr::command(), and BLCMDprintfield::do_grid().
00436 { 00437 if(impl) { 00438 delete impl; 00439 impl = 0; 00440 } 00441 maxline = 128; 00442 current = 1.0; 00443 gradient = 1.0; 00444 normB = 1.0; 00445 normE = 1.0; 00446 00447 BLArgumentVector argv; 00448 BLArgumentMap args; 00449 args["X0"] = d2string(X0); 00450 args["Y0"] = d2string(Y0); 00451 args["Z0"] = d2string(Z0); 00452 args["dX"] = d2string(dX); 00453 args["dY"] = d2string(dY); 00454 args["dZ"] = d2string(dZ); 00455 args["nX"] = i2string(nX); 00456 args["nY"] = i2string(nY); 00457 args["nZ"] = i2string(nZ); 00458 00459 GridImpl *grid = new GridImpl(argv,args); 00460 impl = grid; 00461 00462 bool retval = true; 00463 G4double pos[4], field[6]; 00464 pos[3] = 0.0; 00465 for(int i=0; i<nX; ++i) { 00466 pos[0] = X0 + i*dX; 00467 for(int j=0; j<nY; ++j) { 00468 pos[1] = Y0 + j*dY; 00469 for(int k=0; k<nZ; ++k) { 00470 pos[2] = Z0 + k*dZ; 00471 emField->GetFieldValue(pos,field); 00472 if(!grid->setField(pos[0],pos[1],pos[2], 00473 field[0],field[1],field[2], 00474 field[3],field[4],field[5],0)) 00475 retval = false; 00476 } 00477 } 00478 } 00479 00480 return retval; 00481 }
bool BLFieldMap::createCylinderMap | ( | G4double | Z0, | |
G4double | dR, | |||
G4double | dZ, | |||
int | nR, | |||
int | nZ, | |||
class G4ElectroMagneticField * | field | |||
) |
createCylinderMap() will create a 2-d map from a pre-defined cylindrically-symmetric field. NOTE: cylindrical symmetry is ASSUMED; the field in the Y=0,X>=0 half-plane is used (Bx->Br...). returns false if error.
References current, d2string(), gradient, i2string(), impl, maxline, normB, normE, and CylinderImpl::setField().
Referenced by BLCMDfieldexpr::command(), and BLCMDprintfield::do_cylinder().
00485 { 00486 if(impl) { 00487 delete impl; 00488 impl = 0; 00489 } 00490 maxline = 128; 00491 current = 1.0; 00492 gradient = 1.0; 00493 normB = 1.0; 00494 normE = 1.0; 00495 00496 BLArgumentVector argv; 00497 BLArgumentMap args; 00498 args["Z0"] = d2string(Z0); 00499 args["dR"] = d2string(dR); 00500 args["dZ"] = d2string(dZ); 00501 args["nR"] = i2string(nR); 00502 args["nZ"] = i2string(nZ); 00503 00504 CylinderImpl *cyl = new CylinderImpl(argv,args); 00505 impl = cyl; 00506 00507 // use the Y=0 plane for the R,Z plane 00508 bool retval = true; 00509 G4double pos[4], field[6]; 00510 pos[3] = 0.0; 00511 for(int i=0; i<nR; ++i) { 00512 pos[0] = 0.0 + i*dR; 00513 pos[1] = 0.0; 00514 for(int k=0; k<nZ; ++k) { 00515 pos[2] = Z0 + k*dZ; 00516 emField->GetFieldValue(pos,field); 00517 if(!cyl->setField(pos[0],pos[2],field[0],field[2], 00518 field[3],field[5],0)) 00519 retval = false; 00520 } 00521 } 00522 00523 return false; 00524 }
bool BLFieldMap::createTimeDependence | ( | int | n, | |
G4double | t[], | |||
G4double | b[], | |||
G4double | e[] = 0 , |
|||
G4double | period = -1.0 | |||
) |
createTimeDependence() will apply the time dependence given. n is the # elements in the arrays; t[] is the time value for each point, b[] is the factor for B, and e[] is the factor for E. returns false if error.
References TimeImpl::setPeriod(), and time.
Referenced by BLCMDfieldexpr::handleTimeDependence().
00556 { 00557 if(time) return false; 00558 time = new TimeImpl(n,t,b,e); 00559 if(period > 0.0) 00560 time->setPeriod(period); 00561 return true; 00562 }
bool BLFieldMap::getTimeFactor | ( | G4double | t, | |
G4double * | b, | |||
G4double * | e | |||
) |
getTimeFactor() returns the time factors for B and E at time t. returns false if error.
References TimeImpl::factorB(), TimeImpl::factorE(), and time.
00565 { 00566 if(b) *b = (time ? time->factorB(t) : 1.0); 00567 if(e) *e = (time ? time->factorE(t) : 1.0); 00568 return true; 00569 }
friend class FieldMapPlacement [friend] |
G4int BLFieldMap::maxline [private] |
Referenced by BLFieldMap(), createCylinderMap(), createGridMap(), readFile(), and writeFile().
G4double BLFieldMap::current [private] |
Referenced by BLFieldMap(), createCylinderMap(), createGridMap(), getFieldValue(), readFile(), and writeFile().
G4double BLFieldMap::gradient [private] |
Referenced by BLFieldMap(), createCylinderMap(), createGridMap(), getFieldValue(), readFile(), and writeFile().
G4double BLFieldMap::normB [private] |
Referenced by BLFieldMap(), createCylinderMap(), createGridMap(), getFieldValue(), readFile(), and writeFile().
G4double BLFieldMap::normE [private] |
Referenced by BLFieldMap(), createCylinderMap(), createGridMap(), getFieldValue(), readFile(), and writeFile().
class FieldMapImpl* BLFieldMap::impl [private] |
Referenced by BLFieldMap(), createCylinderMap(), createGridMap(), getBoundingPoint(), getFieldValue(), getFieldValueNoTimeNoScale(), hasB(), hasE(), readFile(), writeFile(), and ~BLFieldMap().
class TimeImpl* BLFieldMap::time [private] |
Referenced by BLFieldMap(), createTimeDependence(), getFieldValue(), getTimeFactor(), readFile(), and ~BLFieldMap().