GridImpl Class Reference

Inheritance diagram for GridImpl:

Detailed Description

class GridImpl -- class for a Grid FieldMap implementation


Public Member Functions
	GridImpl (BLArgumentVector &argv, BLArgumentMap &namedArgs)
	~GridImpl ()
void	getFieldValue (const G4double local[4], G4double field[6]) const
bool	handleCommand (InputFile &in, BLArgumentVector &argv, BLArgumentMap &namedArgs)
virtual void	getBoundingPoint (int i, G4double point[4])
virtual bool	hasB ()
virtual bool	hasE ()
int	bits (G4String s)
const char *	bits2str (int v)
bool	setField (G4double X, G4double Y, G4double Z, G4double Bx, G4double By, G4double Bz, G4double Ex, G4double Ey, G4double Ez, int linenumber)
virtual bool	writeFile (FILE *f)
Private Attributes
G4int	nX
G4int	nY
G4int	nZ
G4double	dX
G4double	dY
G4double	dZ
G4double	X0
G4double	Y0
G4double	Z0
G4double	tolerance
float *	mapBx
float *	mapBy
float *	mapBz
float *	mapEx
float *	mapEy
float *	mapEz
bool	extendX
bool	extendY
bool	extendZ
int	extendXbits
int	extendYbits
int	extendZbits

Constructor & Destructor Documentation

GridImpl::GridImpl	(	BLArgumentVector &	argv,
		BLArgumentMap &	namedArgs
	)

References argDouble(), argInt(), dX, dY, dZ, extendX, extendXbits, extendY, extendYbits, extendZ, extendZbits, mapBx, mapBy, mapBz, mapEx, mapEy, mapEz, nX, nY, nZ, tolerance, X0, Y0, and Z0.

00578                                                         : FieldMapImpl()
00579 {
00580         nX = 2;
00581         nY = 2;
00582         nZ = 2;
00583         dX = 10.0*mm;
00584         dY = 10.0*mm;
00585         dZ = 10.0*mm;
00586         X0 = 0.0;
00587         Y0 = 0.0;
00588         Z0 = 0.0;
00589         tolerance = 0.01*mm;
00590         mapBx = 0;
00591         mapBy = 0;
00592         mapBz = 0;
00593         mapEx = 0;
00594         mapEy = 0;
00595         mapEz = 0;
00596         extendX = false;
00597         extendY = false;
00598         extendZ = false;
00599         extendXbits = 0;
00600         extendYbits = 0;
00601         extendZbits = 0;
00602         argInt(nX,"nX",namedArgs);
00603         argInt(nY,"nY",namedArgs);
00604         argInt(nZ,"nZ",namedArgs);
00605         argDouble(dX,"dX",namedArgs);
00606         argDouble(dY,"dY",namedArgs);
00607         argDouble(dZ,"dZ",namedArgs);
00608         argDouble(X0,"X0",namedArgs);
00609         argDouble(Y0,"Y0",namedArgs);
00610         argDouble(Z0,"Z0",namedArgs);
00611         argDouble(tolerance,"tolerance",namedArgs);
00612 }

GridImpl::~GridImpl ( )

References mapBx, mapBy, mapBz, mapEx, mapEy, and mapEz.

00615 {
00616         if(mapBx) delete mapBx;
00617         if(mapBy) delete mapBy;
00618         if(mapBz) delete mapBz;
00619         if(mapEx) delete mapEx;
00620         if(mapEy) delete mapEy;
00621         if(mapEz) delete mapEz;
00622 }

Member Function Documentation

void GridImpl::getFieldValue	(	const G4double	local[4],
		G4double	field[6]
	)			const `[virtual]`

Implements FieldMapImpl.

References BLAssert, COMPONENT, dX, dY, dZ, extendX, extendXbits, extendY, extendYbits, extendZ, extendZbits, mapBx, mapBy, mapBz, mapEx, mapEy, mapEz, nX, nY, nZ, X0, Y0, and Z0.

00626 {
00627         G4double x = local[0];
00628         G4double y = local[1];
00629         G4double z = local[2];
00630 
00631         x -= X0;
00632         y -= Y0;
00633         z -= Z0;
00634 
00635         G4double factor[6];
00636         factor[0]=factor[1]=factor[2]=factor[3]=factor[4]=factor[5]=1.0;
00637         if(extendX && x < 0.0) {
00638                 x = -x;
00639                 for(int i=0; i<6; ++i) {
00640                         if(extendXbits & (1<<i)) factor[i] = -factor[i];
00641                 }
00642         }
00643         if(extendY && y < 0.0) {
00644                 y = -y;
00645                 for(int i=0; i<6; ++i) {
00646                         if(extendYbits & (1<<i)) factor[i] = -factor[i];
00647                 }
00648         }
00649         if(extendZ && z < 0.0) {
00650                 z = -z;
00651                 for(int i=0; i<6; ++i) {
00652                         if(extendZbits & (1<<i)) factor[i] = -factor[i];
00653                 }
00654         }
00655 
00656         // We compute a 3D linear average of the 8 surrounding points in the map
00657         // First, get the X,Y,Z indices into i,j,k
00658         int i = (int)floor(x/dX);
00659         int j = (int)floor(y/dY);
00660         int k = (int)floor(z/dZ);
00661         if(i < 0 || i >= nX-1 || j < 0 || j >= nY-1 || k < 0 || k >= nZ-1) {
00662                 field[0] = field[1] = field[2] = field[3] = field[4] = 
00663                                                                 field[5] = 0.0;
00664                 return;
00665         }
00666         // m is the initial index (corner of the cube with minimum X, Y, and Z)
00667         int m = k*nY*nX + j*nX + i;
00668         BLAssert(m+nY*nX+nX+1 < nX*nY*nZ);
00669 
00670         // now compute the fractional weighting factors for X, Y, and Z
00671         float fx = 1.0 - (x - i*dX) / dX;
00672         BLAssert(fx >= 0.0 && fx <= 1.0);
00673         float fy = 1.0 - (y - j*dY) / dY;
00674         BLAssert(fy >= 0.0 && fy <= 1.0);
00675         float fz = 1.0 - (z - k*dZ) / dZ;
00676         BLAssert(fz >= 0.0 && fz <= 1.0);
00677 
00678         // now compute the fractional weighting factors for the 8 corners
00679         float f0 = fx*fy*fz;
00680         float f1 = (1.0-fx)*fy*fz;
00681         float f2 = fx*(1.0-fy)*fz;
00682         float f3 = (1.0-fx)*(1.0-fy)*fz;
00683         float f4 = fx*fy*(1.0-fz);
00684         float f5 = (1.0-fx)*fy*(1.0-fz);
00685         float f6 = fx*(1.0-fy)*(1.0-fz);
00686         float f7 = (1.0-fx)*(1.0-fy)*(1.0-fz);
00687 
00688         // Finally, compute the components of the field
00689 #define COMPONENT(C)                                            \
00690         G4double C = 0.0;                                       \
00691         if(map##C) C =  map##C[m]*f0 + map##C[m+1]*f1 +         \
00692                 map##C[m+nX]*f2 + map##C[m+nX+1]*f3 +           \
00693                 map##C[m+nY*nX]*f4 + map##C[m+nY*nX+1]*f5 +     \
00694                 map##C[m+nY*nX+nX]*f6 + map##C[m+nY*nX+nX+1]*f7;
00695         COMPONENT(Bx);
00696         COMPONENT(By);
00697         COMPONENT(Bz);
00698         COMPONENT(Ex);
00699         COMPONENT(Ey);
00700         COMPONENT(Ez);
00701 
00702 #ifdef STUB
00703         G4double Bx = 0.0;
00704         if(mapBx) Bx =
00705                 mapBx[m]*fx*fy*fz +
00706                 mapBx[m+1]*(1.0-fx)*fy*fz +
00707                 mapBx[m+nX]*fx*(1.0-fy)*fz +
00708                 mapBx[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00709                 mapBx[m+nY*nX]*fx*fy*(1.0-fz) +
00710                 mapBx[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00711                 mapBx[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00712                 mapBx[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00713         G4double By = 0.0;
00714         if(mapBy) By =
00715                 mapBy[m]*fx*fy*fz +
00716                 mapBy[m+1]*(1.0-fx)*fy*fz +
00717                 mapBy[m+nX]*fx*(1.0-fy)*fz +
00718                 mapBy[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00719                 mapBy[m+nY*nX]*fx*fy*(1.0-fz) +
00720                 mapBy[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00721                 mapBy[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00722                 mapBy[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00723         G4double Bz = 0.0;
00724         if(mapBz) Bz =
00725                 mapBz[m]*fx*fy*fz +
00726                 mapBz[m+1]*(1.0-fx)*fy*fz +
00727                 mapBz[m+nX]*fx*(1.0-fy)*fz +
00728                 mapBz[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00729                 mapBz[m+nY*nX]*fx*fy*(1.0-fz) +
00730                 mapBz[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00731                 mapBz[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00732                 mapBz[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00733 
00734         G4double Ex = 0.0;
00735         if(mapEx) Ex =
00736                 mapEx[m]*fx*fy*fz +
00737                 mapEx[m+1]*(1.0-fx)*fy*fz +
00738                 mapEx[m+nX]*fx*(1.0-fy)*fz +
00739                 mapEx[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00740                 mapEx[m+nY*nX]*fx*fy*(1.0-fz) +
00741                 mapEx[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00742                 mapEx[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00743                 mapEx[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00744         G4double Ey = 0.0;
00745         if(mapEy) Ey =
00746                 mapEy[m]*fx*fy*fz +
00747                 mapEy[m+1]*(1.0-fx)*fy*fz +
00748                 mapEy[m+nX]*fx*(1.0-fy)*fz +
00749                 mapEy[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00750                 mapEy[m+nY*nX]*fx*fy*(1.0-fz) +
00751                 mapEy[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00752                 mapEy[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00753                 mapEy[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00754         G4double Ez = 0.0;
00755         if(mapEz) Ez =
00756                 mapEz[m]*fx*fy*fz +
00757                 mapEz[m+1]*(1.0-fx)*fy*fz +
00758                 mapEz[m+nX]*fx*(1.0-fy)*fz +
00759                 mapEz[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00760                 mapEz[m+nY*nX]*fx*fy*(1.0-fz) +
00761                 mapEz[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00762                 mapEz[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00763                 mapEz[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00764 #endif //STUB
00765 
00766         field[0] = Bx * factor[0];
00767         field[1] = By * factor[1];
00768         field[2] = Bz * factor[2];
00769         field[3] = Ex * factor[3];
00770         field[4] = Ey * factor[4];
00771         field[5] = Ez * factor[5];
00772 }

bool GridImpl::handleCommand	(	InputFile &	in,
		BLArgumentVector &	argv,
		BLArgumentMap &	namedArgs
	)			`[virtual]`

Implements FieldMapImpl.

References bits(), extendX, extendXbits, extendY, extendYbits, extendZ, extendZbits, InputFile::getline(), InputFile::linenumber(), BLCommand::printError(), InputFile::repeatLine(), and setField().

00776 {
00777         if(argv[0] == "extendX") {
00778                 extendX = true;
00779                 extendXbits = bits(namedArgs["flip"]);
00780                 return true;
00781         } else if(argv[0] == "extendY") {
00782                 extendY = true;
00783                 extendYbits = bits(namedArgs["flip"]);
00784                 return true;
00785         } else if(argv[0] == "extendZ") {
00786                 extendZ = true;
00787                 extendZbits = bits(namedArgs["flip"]);
00788                 return true;
00789         } else if(argv[0] == "Bx") {
00790                 BLCommand::printError("BLFieldMap: Bx not implemented");
00791                 return true;
00792         } else if(argv[0] == "By") {
00793                 BLCommand::printError("BLFieldMap: By not implemented");
00794                 return true;
00795         } else if(argv[0] == "Bz") {
00796                 BLCommand::printError("BLFieldMap: Bz not implemented");
00797                 return true;
00798         } else if(argv[0] == "Ex") {
00799                 BLCommand::printError("BLFieldMap: Ex not implemented");
00800                 return true;
00801         } else if(argv[0] == "Ey") {
00802                 BLCommand::printError("BLFieldMap: Ey not implemented");
00803                 return true;
00804         } else if(argv[0] == "Ez") {
00805                 BLCommand::printError("BLFieldMap: Ez not implemented");
00806                 return true;
00807         } else if(argv[0] == "data") {
00808                 for(char *line=0; (line=in.getline())!=0; ) {
00809                         if(isalpha(line[0])) {
00810                                 in.repeatLine();
00811                                 break;
00812                         }
00813                         int n;
00814                         float X=0.0,Y=0.0,Z=0.0,Bx=0.0,By=0.0,Bz=0.0,
00815                                                         Ex=0.0,Ey=0.0,Ez=0.0;
00816                         for(char *p=line; (p=strchr(p,','))!=0;) *p = ' ';
00817                         n = sscanf(line,"%f%f%f%f%f%f%f%f%f",&X,&Y,&Z,
00818                                                 &Bx,&By,&Bz,&Ex,&Ey,&Ez);
00819                         if(n <= 3) {
00820                                 BLCommand::printError("BLFieldMap: invalid line %d\n",in.linenumber());
00821                                 continue;
00822                         }
00823                         setField(X,Y,Z,Bx*tesla,By*tesla,Bz*tesla,
00824                                 Ex*megavolt/meter,Ey*megavolt/meter,
00825                                 Ez*megavolt/meter,in.linenumber());
00826                 }
00827                 return true;
00828         } else {
00829                 return false;
00830         }
00831 }

void GridImpl::getBoundingPoint	(	int	i,
		G4double	point[4]
	)			`[virtual]`

Implements FieldMapImpl.

References dX, dY, dZ, extendX, extendY, extendZ, nX, nY, nZ, X0, Y0, and Z0.

00834 {
00835         if(extendX)
00836                 point[0] = (i&1 ? -(nX-1)*dX : (nX-1)*dX) + X0;
00837         else
00838                 point[0] = (i&1 ? 0.0 : (nX-1)*dX) + X0;
00839         if(extendY)
00840                 point[1] = (i&2 ? -(nY-1)*dY : (nY-1)*dY) + Y0;
00841         else
00842                 point[1] = (i&2 ? 0.0 : (nY-1)*dY) + Y0;
00843         if(extendZ)
00844                 point[2] = (i&4 ? -(nZ-1)*dZ : (nZ-1)*dZ) + Z0;
00845         else
00846                 point[2] = (i&4 ? 0.0 : (nZ-1)*dZ) + Z0;
00847 }

virtual bool GridImpl::hasB ( ) [inline, virtual]

Implements FieldMapImpl.

References mapBx, mapBy, and mapBz.

00132 { return mapBx!=0 || mapBy!=0 || mapBz!=0; }

virtual bool GridImpl::hasE ( ) [inline, virtual]

Implements FieldMapImpl.

References mapEx, mapEy, and mapEz.

Referenced by writeFile().

00133 { return mapEx!=0 || mapEy!=0 || mapEz!=0; }

int GridImpl::bits ( G4String s )

Referenced by handleCommand().

00850 {
00851         int v=0;
00852         if(s.find("Bx") < s.size()) v |= 1;
00853         if(s.find("By") < s.size()) v |= 2;
00854         if(s.find("Bz") < s.size()) v |= 4;
00855         if(s.find("Ex") < s.size()) v |= 8;
00856         if(s.find("Ey") < s.size()) v |= 16;
00857         if(s.find("Ez") < s.size()) v |= 32;
00858         return v;
00859 }

const char * GridImpl::bits2str ( int v )

Referenced by writeFile().

00862 {
00863         static char retval[32];
00864         retval[0] = '\0';
00865         if(v & 1) strcat(retval,"Bx,");
00866         if(v & 2) strcat(retval,"By,");
00867         if(v & 4) strcat(retval,"Bz,");
00868         if(v & 8) strcat(retval,"Ex,");
00869         if(v & 16) strcat(retval,"Ey,");
00870         if(v & 32) strcat(retval,"Ez,");
00871         return retval;
00872 }

bool GridImpl::setField	(	G4double	X,
		G4double	Y,
		G4double	Z,
		G4double	Bx,
		G4double	By,
		G4double	Bz,
		G4double	Ex,
		G4double	Ey,
		G4double	Ez,
		int	linenumber
	)

References BLAssert, dX, dY, dZ, mapBx, mapBy, mapBz, mapEx, mapEy, mapEz, nX, nY, nZ, BLCommand::printError(), tolerance, X0, Y0, and Z0.

Referenced by BLFieldMap::createGridMap(), and handleCommand().

00877 {
00878         if(!mapBx && Bx != 0.0) {
00879                 mapBx = new float[nX*nY*nZ];
00880                 BLAssert(mapBx != 0);
00881                 for(int i=0; i<nX*nY*nZ; ++i) 
00882                         mapBx[i] = 0.0;
00883         }
00884         if(!mapBy && By != 0.0) {
00885                 mapBy = new float[nX*nY*nZ];
00886                 BLAssert(mapBy != 0);
00887                 for(int i=0; i<nX*nY*nZ; ++i) 
00888                         mapBy[i] = 0.0;
00889         }
00890         if(!mapBz && Bz != 0.0) {
00891                 mapBz = new float[nX*nY*nZ];
00892                 BLAssert(mapBz != 0);
00893                 for(int i=0; i<nX*nY*nZ; ++i) 
00894                         mapBz[i] = 0.0;
00895         }
00896         if(!mapEx && Ex != 0.0) {
00897                 mapEx = new float[nX*nY*nZ];
00898                 BLAssert(mapEx != 0);
00899                 for(int i=0; i<nX*nY*nZ; ++i) 
00900                         mapEx[i] = 0.0;
00901         }
00902         if(!mapEy && Ey != 0.0) {
00903                 mapEy = new float[nX*nY*nZ];
00904                 BLAssert(mapEy != 0);
00905                 for(int i=0; i<nX*nY*nZ; ++i) 
00906                         mapEy[i] = 0.0;
00907         }
00908         if(!mapEz && Ez != 0.0) {
00909                 mapEz = new float[nX*nY*nZ];
00910                 BLAssert(mapEz != 0);
00911                 for(int i=0; i<nX*nY*nZ; ++i) 
00912                         mapEz[i] = 0.0;
00913         }
00914         X -= X0;
00915         Y -= Y0;
00916         Z -= Z0;
00917         int i = (int)floor((X/dX) + 0.5);
00918         if(i<0 || fabs(i*dX-X)>tolerance || i >= nX) {
00919                 BLCommand::printError("BLFieldMap: ERROR point off"
00920                         " grid  X=%.2f line=%d\n",X, linenumber);
00921                 return false;
00922         }
00923         int j = (int)floor((Y/dY) + 0.5);
00924         if(j<0 || fabs(j*dY-Y)>tolerance || j >= nY) {
00925                 BLCommand::printError("BLFieldMap: ERROR point off"
00926                         " grid Y=%.2f line=%d\n",Y, linenumber);
00927                 return false;
00928         }
00929         int k = (int)floor((Z/dZ) + 0.5);
00930         if(k<0 || fabs(k*dZ-Z)>tolerance || k >= nZ) {
00931                 BLCommand::printError("BLFieldMap: ERROR point off"
00932                         " grid X=%.2f line=%d\n",Z, linenumber);
00933                 return false;
00934         }
00935         int m = k*nY*nX + j*nX + i;
00936         BLAssert(m >= 0 && m < nX*nY*nZ);
00937         if(mapBx) mapBx[m] = Bx;
00938         if(mapBy) mapBy[m] = By;
00939         if(mapBz) mapBz[m] = Bz;
00940         if(mapEx) mapEx[m] = Ex;
00941         if(mapEy) mapEy[m] = Ey;
00942         if(mapEz) mapEz[m] = Ez;
00943 
00944         return true;
00945 }

bool GridImpl::writeFile ( FILE * f ) [virtual]

Implements FieldMapImpl.

References bits2str(), BLAssert, dX, dY, dZ, extendX, extendXbits, extendY, extendYbits, extendZ, extendZbits, hasE(), mapBx, mapBy, mapBz, mapEx, mapEy, mapEz, nX, nY, nZ, X0, Y0, and Z0.

00948 {
00949         fprintf(f,"grid X0=%g Y0=%g Z0=%g nX=%d nY=%d nZ=%d dX=%g dY=%g dZ=%g\n",
00950                 X0,Y0,Z0,nX,nY,nZ,dX,dY,dZ);
00951         if(extendX) fprintf(f,"extendX flip=%s\n",bits2str(extendXbits));
00952         if(extendY) fprintf(f,"extendY flip=%s\n",bits2str(extendYbits));
00953         if(extendZ) fprintf(f,"extendZ flip=%s\n",bits2str(extendZbits));
00954         fprintf(f,"data\n");
00955 
00956         for(int i=0; i<nX; ++i) {
00957                 G4double X = X0 + i*dX;
00958                 for(int j=0; j<nY; ++j) {
00959                         G4double Y = Y0 + j*dY;
00960                         for(int k=0; k<nZ; ++k) {
00961                                 G4double Z = Z0 + k*dZ;
00962                                 int m = k*nY*nX + j*nX + i;
00963                                 BLAssert(m >= 0 && m < nX*nY*nZ);
00964                                 G4double Bx = (mapBx ? mapBx[m] : 0.0);
00965                                 G4double By = (mapBy ? mapBy[m] : 0.0);
00966                                 G4double Bz = (mapBz ? mapBz[m] : 0.0);
00967                                 fprintf(f,"%.1f %.1f %.1f %g %g %g",
00968                                         X,Y,Z,Bx/tesla,By/tesla,Bz/tesla);
00969                                 if(hasE()) {
00970                                         G4double Ex = (mapEx ? mapEx[m] : 0.0);
00971                                         G4double Ey = (mapEy ? mapEy[m] : 0.0);
00972                                         G4double Ez = (mapEz ? mapEz[m] : 0.0);
00973                                         fprintf(f," %g %g %g",
00974                                                 Ex/(megavolt/meter),
00975                                                 Ey/(megavolt/meter),
00976                                                 Ez/(megavolt/meter));
00977                                 }
00978                                 fprintf(f,"\n");
00979                         }
00980                 }
00981         }
00982         return true;
00983 }