g4sim: src/MyBLFieldMap.cc Source File

00001 // AK20110616 Modified from G4beamline v2.03 BLFieldMap.cc
00002 //
00003 /*
00004 This source file is part of G4beamline, http://g4beamline.muonsinc.com
00005 Copyright (C) 2003,2004,2005,2006 by Tom Roberts, all rights reserved.
00006 
00007 This program is free software; you can redistribute it and/or
00008 modify it under the terms of the GNU General Public License
00009 as published by the Free Software Foundation; either version 2
00010 of the License, or (at your option) any later version.
00011 
00012 This program is distributed in the hope that it will be useful,
00013 but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 GNU General Public License for more details.
00016 
00017 http://www.gnu.org/copyleft/gpl.html
00018 */
00019 
00020 #include <fstream>
00021 #include <string.h>
00022 #include <ctype.h>
00023 #include <assert.h>
00024 
00025 #include "G4ElectroMagneticField.hh"
00026 #include "MyBLFieldMap.hh"
00027 #include "MyBLFuncs.hh"
00028 #include "BTSpline1D.hh"
00029 
00032 class InputFile {
00033     std::ifstream in;
00034     G4String name;
00035     char *line;
00036     int maxline;
00037     int lineno;
00038     bool repeatline;
00039     public:
00040         InputFile(G4String filename, G4int _maxline=1024) : in(), name() {
00041             name = filename;
00042             in.open(name.c_str());
00043             maxline = _maxline;
00044             line = new char[maxline];
00045             assert(line != 0);
00046             lineno = 0;
00047             repeatline = false;
00048         }
00049         ~InputFile() { close(); }
00050         const char *filename() { return name.c_str(); }
00051         bool good() { return in.good(); }
00052         int linenumber() { return lineno; }
00053         char *getline() {                             // skips blank lines and comments
00054             if(repeatline) {
00055                 repeatline = false;
00056                 return line;
00057             }
00058             while(in.good()) {
00059                 line[0] = '\0';
00060                 in.getline(line,maxline);
00061                 ++lineno;
00062                 if(line[0] == '*')
00063                     printf("%s\n",line);
00064                 if(line[0] == '\0' || line[0] == '#' || line[0] == '*')
00065                     continue;
00066                 return line;
00067             }
00068             return 0;
00069         }
00070         void setMaxline(int _maxline) {
00071             maxline = _maxline;
00072             assert(maxline >= 80);
00073             if(line) delete line;
00074             line = new char[maxline];
00075             assert(line != 0);
00076         }
00077         void close() {
00078             if(line) {
00079                 in.close(); delete line; line = 0;
00080             }
00081         }
00082         void repeatLine() { repeatline = true; }
00083 };
00084 
00087 class FieldMapImpl {
00088     public:
00089         FieldMapImpl() { }
00090         virtual ~FieldMapImpl() { }
00091         virtual void getFieldValue(const G4double local[4], G4double field[6])
00092             const = 0;
00093         virtual bool handleCommand(InputFile &in, MyBLArgumentVector &argv,
00094             MyBLArgumentMap &namedArgs) = 0;
00095         virtual void getBoundingPoint(int i, G4double point[4]) = 0;
00096         virtual bool hasB() = 0;
00097         virtual bool hasE() = 0;
00098         bool readBlock(InputFile &in, float *values, int nRows, int nCols,
00099             G4double units);
00100         virtual bool writeFile(FILE *f) = 0;
00101 };
00102 
00105 class GridImpl : public FieldMapImpl {
00106     G4int nX;
00107     G4int nY;
00108     G4int nZ;
00109     G4double dX;
00110     G4double dY;
00111     G4double dZ;
00112     G4double X0;
00113     G4double Y0;
00114     G4double Z0;
00115     G4double tolerance;
00116     float *mapBx;
00117     float *mapBy;
00118     float *mapBz;
00119     float *mapEx;
00120     float *mapEy;
00121     float *mapEz;
00122     bool extendX;
00123     bool extendY;
00124     bool extendZ;
00125     int extendXbits;
00126     int extendYbits;
00127     int extendZbits;
00128     public:
00129         GridImpl(MyBLArgumentVector &argv, MyBLArgumentMap &namedArgs);
00130         ~GridImpl();
00131         void getFieldValue(const G4double local[4], G4double field[6]) const;
00132         bool handleCommand(InputFile &in, MyBLArgumentVector &argv,
00133             MyBLArgumentMap &namedArgs);
00134         virtual void getBoundingPoint(int i, G4double point[4]);
00135         virtual bool hasB() { return mapBx!=0 || mapBy!=0 || mapBz!=0; }
00136         virtual bool hasE() { return mapEx!=0 || mapEy!=0 || mapEz!=0; }
00137         int bits(G4String s);
00138         const char *bits2str(int v);
00139         bool setField(G4double X, G4double Y, G4double Z, G4double Bx,
00140             G4double By, G4double Bz, G4double Ex, G4double Ey,
00141             G4double Ez, int linenumber);
00142         virtual bool writeFile(FILE *f);
00143 };
00144 
00147 class CylinderImpl : public FieldMapImpl {
00148     G4int nR;
00149     G4int nZ;
00150     G4double dR;
00151     G4double dZ;
00152     G4double Z0;
00153     G4double tolerance;
00154     float *mapBr;
00155     float *mapBz;
00156     float *mapEr;
00157     float *mapEz;
00158     bool extendZ;
00159     float extendBrFactor, extendBzFactor;
00160     float extendErFactor, extendEzFactor;
00161     public:
00162         CylinderImpl(MyBLArgumentVector &argv, MyBLArgumentMap &namedArgs);
00163         ~CylinderImpl();
00164         void getFieldValue(const G4double local[4], G4double field[6]) const;
00165         bool handleCommand(InputFile &in, MyBLArgumentVector &argv,
00166             MyBLArgumentMap &namedArgs);
00167         virtual void getBoundingPoint(int i, G4double point[4]);
00168         virtual bool hasB() { return mapBz != 0 || mapBr != 0; }
00169         virtual bool hasE() { return mapEz != 0 || mapEr != 0; }
00170         bool setField(G4double R, G4double Z, G4double Br,
00171             G4double Bz, G4double Er, G4double Ez, int linenumber);
00172         virtual bool writeFile(FILE *f);
00173 };
00174 
00177 class TimeImpl {
00178     BTSpline1D B;
00179     BTSpline1D E;
00180     G4double tMin;
00181     G4double tMax;
00182     G4double period;
00183     public:
00189         TimeImpl(int n, G4double t[], G4double b[], G4double e[]=0) :
00190         B(n,t,b), E(n,t,(e?e:b)) {
00191             tMin = t[0];
00192             tMax = t[n-1];
00193             period = -1.0;
00194         }
00195         void setPeriod(G4double v) { period = v; }
00196         G4double factorB(G4double t) {
00197             if(period > 0.0) {
00198                 t -= floor(t/period)*period;
00199             }
00200             else {
00201                 if(t < tMin) t = tMin;
00202                 if(t > tMax) t = tMax;
00203             }
00204             return B(t);
00205         }
00206         G4double factorE(G4double t) {
00207             if(period > 0.0) {
00208                 t -= floor(t/period)*period;
00209             }
00210             else {
00211                 if(t < tMin) t = tMin;
00212                 if(t > tMax) t = tMax;
00213             }
00214             return E(t);
00215         }
00216         static TimeImpl *readTime(InputFile &in, MyBLArgumentVector argv,
00217             MyBLArgumentMap namedArgs);
00218 };
00219 
00221 static void argDouble(G4double &var, const char *name, MyBLArgumentMap &namedArgs) {
00222     if(namedArgs.count(name) == 0) return;
00223     char *q;
00224     const char *p=namedArgs[name].c_str();
00225     G4double v = strtod(p,&q);
00226     if(p == q || *q != '\0') {
00227         G4cerr << "Invalid value for "<< name<<" in MyBLFieldMap"
00228             << G4endl;
00229 
00230         return;
00231     }
00232     var = v;
00233 }
00234 
00235 
00237 static void argInt(G4int &var, const char *name, MyBLArgumentMap &namedArgs) {
00238     if(namedArgs.count(name) == 0) return;
00239     char *q;
00240     const char *p=namedArgs[name].c_str();
00241     G4int v = strtol(p,&q,0);
00242     if(p == q || *q != '\0') {
00243         G4cerr << "Invalid value for "<< name<<" in MyBLFieldMap"
00244             << G4endl;
00245         return;
00246     }
00247     var = v;
00248 }
00249 
00250 
00252 static G4String d2string(G4double v) {
00253     char tmp[32];
00254     sprintf(tmp,"%.8g",v);
00255     return G4String(tmp);
00256 }
00257 
00258 
00260 static G4String i2string(int v) {
00261     char tmp[32];
00262     sprintf(tmp,"%d",v);
00263     return G4String(tmp);
00264 }
00265 
00266 
00267 MyBLFieldMap::MyBLFieldMap() {
00268     maxline = 1024;
00269     current = 1.0;
00270     gradient = 1.0;
00271     normB = 1.0;
00272     normE = 1.0;
00273     impl = 0;
00274     time = 0;
00275 }
00276 
00277 
00278 MyBLFieldMap::~MyBLFieldMap() {
00279     if(impl) delete impl;
00280     impl = 0;
00281     if(time) delete time;
00282     time = 0;
00283 }
00284 
00285 
00286 void MyBLFieldMap::getFieldValue(const G4double local[4], G4double field[6],
00287 G4double _current, G4double _gradient) {
00288     if(!impl)
00289         throw "MyBLFieldMap::getFieldValue called, no implementation";
00290 
00291     G4double thisField[6];
00292     impl->getFieldValue(local,thisField);
00293     G4double timeB=1.0, timeE=1.0;
00294     if(time) {
00295         timeB = time->factorB(local[3]);
00296         timeE = time->factorE(local[3]);
00297     }
00298     field[0] = thisField[0] * normB * timeB * _current/current;
00299     field[1] = thisField[1] * normB * timeB * _current/current;
00300     field[2] = thisField[2] * normB * timeB * _current/current;
00301     field[3] = thisField[3] * normE * timeE * _gradient/gradient;
00302     field[4] = thisField[4] * normE * timeE * _gradient/gradient;
00303     field[5] = thisField[5] * normE * timeE * _gradient/gradient;
00304 }
00305 
00306 
00307 bool MyBLFieldMap::readFile(G4String filename) {
00308     InputFile in(filename,maxline);
00309     if(!in.good()) {
00310         throw (G4String("MyBLFieldMap::readFile: cannot open file ") + in.filename());
00311 
00312     }
00313     printf("MyBLFieldMap: reading file '%s'\n",in.filename());
00314 
00315     bool retval = true;
00316 
00317     char *line;
00318     while((line=in.getline()) != 0) {
00319         MyBLArgumentVector argv;
00320         MyBLArgumentMap namedArgs;
00321         if(BLFuncs.parseArgs(line,argv,namedArgs) < 0)
00322             goto invalid;
00323         if(argv[0] == "") continue;
00324         if(argv[0] == "param") {
00325             argInt(maxline,"maxline",namedArgs);
00326             argDouble(current,"current",namedArgs);
00327             argDouble(gradient,"gradient",namedArgs);
00328             argDouble(normB,"normB",namedArgs);
00329             argDouble(normE,"normE",namedArgs);
00330             in.setMaxline(maxline);
00331         }
00332         else if(argv[0] == "grid") {
00333             if(impl) goto invalid;
00334             impl = new GridImpl(argv,namedArgs);
00335         }
00336         else if(argv[0] == "cylinder") {
00337             if(impl) goto invalid;
00338             impl = new CylinderImpl(argv,namedArgs);
00339         }
00340         else if(argv[0] == "time") {
00341             if(time) goto invalid;
00342             time = TimeImpl::readTime(in,argv,namedArgs);
00343             if(!time) goto invalid;
00344         }
00345         else if(impl) {
00346             if(!impl->handleCommand(in,argv,namedArgs))
00347                 goto invalid;
00348         }
00349         else {
00350             invalid:        G4cerr << "MyBLFieldMap file "  << in.filename()
00351                 << " line " << in.linenumber()
00352                 << " invalid command " << argv[0].c_str()
00353                 << G4endl;
00354             retval = false;
00355             break;
00356         }
00357     }
00358     in.close();
00359 
00360     return retval;
00361 }
00362 
00363 
00364 void MyBLFieldMap::getBoundingPoint(int i, G4double point[4]) {
00365     impl->getBoundingPoint(i,point);
00366 }
00367 
00368 
00369 bool MyBLFieldMap::hasB() {
00370     return impl->hasB();
00371 }
00372 
00373 
00374 bool MyBLFieldMap::hasE() {
00375     return impl->hasE();
00376 }
00377 
00378 
00379 TimeImpl *TimeImpl::readTime(InputFile &in, MyBLArgumentVector argv,
00380 MyBLArgumentMap namedArgs) {
00381     G4double period=-1.0;
00382     argDouble(period,"period",namedArgs);
00383 
00384     std::vector<G4double> T;
00385     std::vector<G4double> B;
00386     std::vector<G4double> E;
00387 
00388     char *line;
00389     while((line=in.getline()) != 0) {
00390         if(isalpha(line[0])) {
00391             in.repeatLine();
00392             break;
00393         }
00394         G4double t,b,e;
00395         int j = sscanf(line,"%lf%lf%lf",&t,&b,&e);
00396         if(j < 2 || j > 3) return 0;
00397         if(j == 2) e = b;
00398         T.push_back(t);
00399         B.push_back(b);
00400         E.push_back(e);
00401     }
00402 
00403     if(T.size() > 0) {
00404         TimeImpl *impl =
00405             new TimeImpl(T.size(),&T[0],&B[0],&E[0]);
00406         impl->setPeriod(period);
00407         return impl;
00408     }
00409 
00410     return 0;
00411 }
00412 
00413 
00414 bool FieldMapImpl::readBlock(InputFile &in, float *values, int nRows, int nCols,
00415 G4double units) {
00416     while(nRows-- > 0) {
00417         char *line = in.getline();
00418         if(!line) return false;
00419         char *p=line;
00420         for(int i=0; i<nCols; ++i) {
00421             while(isspace(*p)) ++p;
00422             if(*p == '\0') return false;
00423             *values++ = strtod(p,&p) * units;
00424             if(*p == ',') ++p;
00425         }
00426     }
00427     return true;
00428 }
00429 
00430 
00431 bool MyBLFieldMap::createGridMap(G4double X0, G4double Y0, G4double Z0,
00432 G4double dX, G4double dY, G4double dZ, int nX, int nY, int nZ,
00433 G4ElectroMagneticField *emField) {
00434     if(impl) {
00435         delete impl;
00436         impl = 0;
00437     }
00438     maxline = 128;
00439     current = 1.0;
00440     gradient = 1.0;
00441     normB = 1.0;
00442     normE = 1.0;
00443 
00444     MyBLArgumentVector argv;
00445     MyBLArgumentMap args;
00446     args["X0"] = d2string(X0);
00447     args["Y0"] = d2string(Y0);
00448     args["Z0"] = d2string(Z0);
00449     args["dX"] = d2string(dX);
00450     args["dY"] = d2string(dY);
00451     args["dZ"] = d2string(dZ);
00452     args["nX"] = i2string(nX);
00453     args["nY"] = i2string(nY);
00454     args["nZ"] = i2string(nZ);
00455 
00456     GridImpl *grid = new GridImpl(argv,args);
00457     impl = grid;
00458 
00459     bool retval = true;
00460     G4double pos[4], field[6];
00461     pos[3] = 0.0;
00462     for(int i=0; i<nX; ++i) {
00463         pos[0] = X0 + i*dX;
00464         for(int j=0; j<nY; ++j) {
00465             pos[1] = Y0 + j*dY;
00466             for(int k=0; k<nZ; ++k) {
00467                 pos[2] = Z0 + k*dZ;
00468                 emField->GetFieldValue(pos,field);
00469                 if(!grid->setField(pos[0],pos[1],pos[2],
00470                     field[0],field[1],field[2],
00471                     field[3],field[4],field[5],0))
00472                     retval = false;
00473             }
00474         }
00475     }
00476 
00477     return retval;
00478 }
00479 
00480 
00481 bool MyBLFieldMap::createCylinderMap(G4double Z0, G4double dR,  G4double dZ,
00482 int nR, int nZ, G4ElectroMagneticField *emField) {
00483     if(impl) {
00484         delete impl;
00485         impl = 0;
00486     }
00487     maxline = 128;
00488     current = 1.0;
00489     gradient = 1.0;
00490     normB = 1.0;
00491     normE = 1.0;
00492 
00493     MyBLArgumentVector argv;
00494     MyBLArgumentMap args;
00495     args["Z0"] = d2string(Z0);
00496     args["dR"] = d2string(dR);
00497     args["dZ"] = d2string(dZ);
00498     args["nR"] = i2string(nR);
00499     args["nZ"] = i2string(nZ);
00500 
00501     CylinderImpl *cyl = new CylinderImpl(argv,args);
00502     impl = cyl;
00503 
00504     // use the Y=0 plane for the R,Z plane
00505     bool retval = true;
00506     G4double pos[4], field[6];
00507     pos[3] = 0.0;
00508     for(int i=0; i<nR; ++i) {
00509         pos[0] = 0.0 + i*dR;
00510         pos[1] = 0.0;
00511         for(int k=0; k<nZ; ++k) {
00512             pos[2] = Z0 + k*dZ;
00513             emField->GetFieldValue(pos,field);
00514             if(!cyl->setField(pos[0],pos[2],field[0],field[2],
00515                 field[3],field[5],0))
00516                 retval = false;
00517         }
00518     }
00519 
00520     return false;
00521 }
00522 
00523 
00524 bool MyBLFieldMap::writeFile(G4String filename, G4String comment) {
00525     if(!impl) return false;
00526 
00527     FILE *f = fopen(filename.c_str(),"r");
00528     if(f) {
00529         fclose(f);
00530         G4Exception("MyBLFieldMap","Output File Exists",FatalException,
00531             filename.c_str());
00532     }
00533 
00534     f = fopen(filename,"w");
00535     if(!f) {
00536         fprintf(stderr,"MyBLFieldMap::writeFile CANNOT WRITE file '%s'\n",
00537             filename.c_str());
00538         return false;
00539     }
00540 
00541     fprintf(f,"# %s\n",comment.c_str());
00542     fprintf(f,"param maxline=%d current=%g gradient=%g normB=%g normE=%g\n",
00543         maxline,current,gradient,normB,normE);
00544 
00545     bool retval = impl->writeFile(f);
00546 
00547     fclose(f);
00548     return retval;
00549 }
00550 
00551 
00552 bool MyBLFieldMap::createTimeDependence(int n, G4double t[], G4double b[],
00553 G4double e[], G4double period) {
00554     if(time) return false;
00555     time = new TimeImpl(n,t,b,e);
00556     if(period > 0.0)
00557         time->setPeriod(period);
00558     return true;
00559 }
00560 
00561 
00562 bool MyBLFieldMap::getTimeFactor(G4double t, G4double *b, G4double *e) {
00563     if(b) *b = (time ? time->factorB(t) : 1.0);
00564     if(e) *e = (time ? time->factorE(t) : 1.0);
00565     return true;
00566 }
00567 
00568 
00569 GridImpl::GridImpl(MyBLArgumentVector &argv, MyBLArgumentMap &namedArgs)
00570 : FieldMapImpl() {
00571     nX = 2;
00572     nY = 2;
00573     nZ = 2;
00574     dX = 10.0*mm;
00575     dY = 10.0*mm;
00576     dZ = 10.0*mm;
00577     X0 = 0.0;
00578     Y0 = 0.0;
00579     Z0 = 0.0;
00580     tolerance = 0.01*mm;
00581     mapBx = 0;
00582     mapBy = 0;
00583     mapBz = 0;
00584     mapEx = 0;
00585     mapEy = 0;
00586     mapEz = 0;
00587     extendX = false;
00588     extendY = false;
00589     extendZ = false;
00590     extendXbits = 0;
00591     extendYbits = 0;
00592     extendZbits = 0;
00593     argInt(nX,"nX",namedArgs);
00594     argInt(nY,"nY",namedArgs);
00595     argInt(nZ,"nZ",namedArgs);
00596     argDouble(dX,"dX",namedArgs);
00597     argDouble(dY,"dY",namedArgs);
00598     argDouble(dZ,"dZ",namedArgs);
00599     argDouble(X0,"X0",namedArgs);
00600     argDouble(Y0,"Y0",namedArgs);
00601     argDouble(Z0,"Z0",namedArgs);
00602     argDouble(tolerance,"tolerance",namedArgs);
00603 }
00604 
00605 
00606 GridImpl::~GridImpl() {
00607     if(mapBx) delete mapBx;
00608     if(mapBy) delete mapBy;
00609     if(mapBz) delete mapBz;
00610     if(mapEx) delete mapEx;
00611     if(mapEy) delete mapEy;
00612     if(mapEz) delete mapEz;
00613 }
00614 
00615 
00616 void GridImpl::getFieldValue(const G4double local[4], G4double field[6])
00617 const {
00618     G4double x = local[0];
00619     G4double y = local[1];
00620     G4double z = local[2];
00621 
00622     x -= X0;
00623     y -= Y0;
00624     z -= Z0;
00625 
00626     G4double factor[6];
00627     factor[0]=factor[1]=factor[2]=factor[3]=factor[4]=factor[5]=1.0;
00628     if(extendX && x < 0.0) {
00629         x = -x;
00630         for(int i=0; i<6; ++i) {
00631             if(extendXbits & (1<<i)) factor[i] = -factor[i];
00632         }
00633     }
00634     if(extendY && y < 0.0) {
00635         y = -y;
00636         for(int i=0; i<6; ++i) {
00637             if(extendYbits & (1<<i)) factor[i] = -factor[i];
00638         }
00639     }
00640     if(extendZ && z < 0.0) {
00641         z = -z;
00642         for(int i=0; i<6; ++i) {
00643             if(extendZbits & (1<<i)) factor[i] = -factor[i];
00644         }
00645     }
00646 
00647     // We compute a 3D linear average of the 8 surrounding points in the map
00648     // First, get the X,Y,Z indices into i,j,k
00649     int i = (int)floor(x/dX);
00650     int j = (int)floor(y/dY);
00651     int k = (int)floor(z/dZ);
00652     if(i < 0 || i >= nX-1 || j < 0 || j >= nY-1 || k < 0 || k >= nZ-1) {
00653         field[0] = field[1] = field[2] = field[3] = field[4] =
00654             field[5] = 0.0;
00655         return;
00656     }
00657     // m is the initial index (corner of the cube with minimum X, Y, and Z)
00658     int m = k*nY*nX + j*nX + i;
00659     assert(m+nY*nX+nX+1 < nX*nY*nZ);
00660 
00661     // now compute the fractional weighting factors for X, Y, and Z
00662     float fx = 1.0 - (x - i*dX) / dX;
00663     assert(fx >= 0.0 && fx <= 1.0);
00664     float fy = 1.0 - (y - j*dY) / dY;
00665     assert(fy >= 0.0 && fy <= 1.0);
00666     float fz = 1.0 - (z - k*dZ) / dZ;
00667     assert(fz >= 0.0 && fz <= 1.0);
00668 
00669     // now compute the fractional weighting factors for the 8 corners
00670     float f0 = fx*fy*fz;
00671     float f1 = (1.0-fx)*fy*fz;
00672     float f2 = fx*(1.0-fy)*fz;
00673     float f3 = (1.0-fx)*(1.0-fy)*fz;
00674     float f4 = fx*fy*(1.0-fz);
00675     float f5 = (1.0-fx)*fy*(1.0-fz);
00676     float f6 = fx*(1.0-fy)*(1.0-fz);
00677     float f7 = (1.0-fx)*(1.0-fy)*(1.0-fz);
00678 
00679     // Finally, compute the components of the field
00680 #define COMPONENT(C) \
00681         G4double C = 0.0; \
00682         if(map##C) C =  map##C[m]*f0 + map##C[m+1]*f1 + \
00683                 map##C[m+nX]*f2 + map##C[m+nX+1]*f3 + \
00684                 map##C[m+nY*nX]*f4 + map##C[m+nY*nX+1]*f5 + \
00685                 map##C[m+nY*nX+nX]*f6 + map##C[m+nY*nX+nX+1]*f7;
00686         COMPONENT(Bx);
00687     COMPONENT(By);
00688     COMPONENT(Bz);
00689     COMPONENT(Ex);
00690     COMPONENT(Ey);
00691     COMPONENT(Ez);
00692 
00693 #ifdef STUB
00694     G4double Bx = 0.0;
00695     if(mapBx) Bx =
00696             mapBx[m]*fx*fy*fz +
00697             mapBx[m+1]*(1.0-fx)*fy*fz +
00698             mapBx[m+nX]*fx*(1.0-fy)*fz +
00699             mapBx[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00700             mapBx[m+nY*nX]*fx*fy*(1.0-fz) +
00701             mapBx[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00702             mapBx[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00703             mapBx[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00704     G4double By = 0.0;
00705     if(mapBy) By =
00706             mapBy[m]*fx*fy*fz +
00707             mapBy[m+1]*(1.0-fx)*fy*fz +
00708             mapBy[m+nX]*fx*(1.0-fy)*fz +
00709             mapBy[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00710             mapBy[m+nY*nX]*fx*fy*(1.0-fz) +
00711             mapBy[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00712             mapBy[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00713             mapBy[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00714     G4double Bz = 0.0;
00715     if(mapBz) Bz =
00716             mapBz[m]*fx*fy*fz +
00717             mapBz[m+1]*(1.0-fx)*fy*fz +
00718             mapBz[m+nX]*fx*(1.0-fy)*fz +
00719             mapBz[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00720             mapBz[m+nY*nX]*fx*fy*(1.0-fz) +
00721             mapBz[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00722             mapBz[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00723             mapBz[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00724 
00725     G4double Ex = 0.0;
00726     if(mapEx) Ex =
00727             mapEx[m]*fx*fy*fz +
00728             mapEx[m+1]*(1.0-fx)*fy*fz +
00729             mapEx[m+nX]*fx*(1.0-fy)*fz +
00730             mapEx[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00731             mapEx[m+nY*nX]*fx*fy*(1.0-fz) +
00732             mapEx[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00733             mapEx[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00734             mapEx[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00735     G4double Ey = 0.0;
00736     if(mapEy) Ey =
00737             mapEy[m]*fx*fy*fz +
00738             mapEy[m+1]*(1.0-fx)*fy*fz +
00739             mapEy[m+nX]*fx*(1.0-fy)*fz +
00740             mapEy[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00741             mapEy[m+nY*nX]*fx*fy*(1.0-fz) +
00742             mapEy[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00743             mapEy[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00744             mapEy[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00745     G4double Ez = 0.0;
00746     if(mapEz) Ez =
00747             mapEz[m]*fx*fy*fz +
00748             mapEz[m+1]*(1.0-fx)*fy*fz +
00749             mapEz[m+nX]*fx*(1.0-fy)*fz +
00750             mapEz[m+nX+1]*(1.0-fx)*(1.0-fy)*fz +
00751             mapEz[m+nY*nX]*fx*fy*(1.0-fz) +
00752             mapEz[m+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
00753             mapEz[m+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
00754             mapEz[m+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
00755 #endif                                            //STUB
00756 
00757     field[0] = Bx * factor[0];
00758     field[1] = By * factor[1];
00759     field[2] = Bz * factor[2];
00760     field[3] = Ex * factor[3];
00761     field[4] = Ey * factor[4];
00762     field[5] = Ez * factor[5];
00763 }
00764 
00765 
00766 bool GridImpl::handleCommand(InputFile &in, MyBLArgumentVector &argv,
00767 MyBLArgumentMap &namedArgs) {
00768     if(argv[0] == "extendX") {
00769         extendX = true;
00770         extendXbits = bits(namedArgs["flip"]);
00771         return true;
00772     }
00773     else if(argv[0] == "extendY") {
00774         extendY = true;
00775         extendYbits = bits(namedArgs["flip"]);
00776         return true;
00777     }
00778     else if(argv[0] == "extendZ") {
00779         extendZ = true;
00780         extendZbits = bits(namedArgs["flip"]);
00781         return true;
00782     }
00783     else if(argv[0] == "Bx") {
00784         G4cerr << "MyBLFieldMap: Bx not implemented" << G4endl;
00785         return true;
00786     }
00787     else if(argv[0] == "By") {
00788         G4cerr << "MyBLFieldMap: By not implemented" << G4endl;
00789         return true;
00790     }
00791     else if(argv[0] == "Bz") {
00792         G4cerr << "MyBLFieldMap: Bz not implemented" << G4endl;
00793         return true;
00794     }
00795     else if(argv[0] == "Ex") {
00796         G4cerr << "MyBLFieldMap: Ex not implemented" << G4endl;
00797         return true;
00798     }
00799     else if(argv[0] == "Ey") {
00800         G4cerr << "MyBLFieldMap: Ey not implemented" << G4endl;
00801         return true;
00802     }
00803     else if(argv[0] == "Ez") {
00804         G4cerr << "MyBLFieldMap: Ez not implemented" << G4endl;
00805         return true;
00806     }
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                 G4cerr << "MyBLFieldMap: invalid line "
00821                     << in.linenumber() << G4endl;
00822                 continue;
00823             }
00824             setField(X,Y,Z,Bx*tesla,By*tesla,Bz*tesla,
00825                 Ex*megavolt/meter,Ey*megavolt/meter,
00826                 Ez*megavolt/meter,in.linenumber());
00827         }
00828         return true;
00829     }
00830     else {
00831         return false;
00832     }
00833 }
00834 
00835 
00836 void GridImpl::getBoundingPoint(int i, G4double point[4]) {
00837     if(extendX)
00838         point[0] = (i&1 ? -(nX-1)*dX : (nX-1)*dX) + X0;
00839     else
00840         point[0] = (i&1 ? 0.0 : (nX-1)*dX) + X0;
00841     if(extendY)
00842         point[1] = (i&2 ? -(nY-1)*dY : (nY-1)*dY) + Y0;
00843     else
00844         point[1] = (i&2 ? 0.0 : (nY-1)*dY) + Y0;
00845     if(extendZ)
00846         point[2] = (i&4 ? -(nZ-1)*dZ : (nZ-1)*dZ) + Z0;
00847     else
00848         point[2] = (i&4 ? 0.0 : (nZ-1)*dZ) + Z0;
00849 }
00850 
00851 
00852 int GridImpl::bits(G4String s) {
00853     int v=0;
00854     if(s.find("Bx") < s.size()) v |= 1;
00855     if(s.find("By") < s.size()) v |= 2;
00856     if(s.find("Bz") < s.size()) v |= 4;
00857     if(s.find("Ex") < s.size()) v |= 8;
00858     if(s.find("Ey") < s.size()) v |= 16;
00859     if(s.find("Ez") < s.size()) v |= 32;
00860     return v;
00861 }
00862 
00863 
00864 const char *GridImpl::bits2str(int v) {
00865     static char retval[32];
00866     retval[0] = '\0';
00867     if(v & 1) strcat(retval,"Bx,");
00868     if(v & 2) strcat(retval,"By,");
00869     if(v & 4) strcat(retval,"Bz,");
00870     if(v & 8) strcat(retval,"Ex,");
00871     if(v & 16) strcat(retval,"Ey,");
00872     if(v & 32) strcat(retval,"Ez,");
00873     return retval;
00874 }
00875 
00876 
00877 bool GridImpl::setField(G4double X, G4double Y, G4double Z, G4double Bx,
00878 G4double By, G4double Bz, G4double Ex, G4double Ey, G4double Ez,
00879 int linenumber) {
00880     if(!mapBx && Bx != 0.0) {
00881         mapBx = new float[nX*nY*nZ];
00882         assert(mapBx != 0);
00883         for(int i=0; i<nX*nY*nZ; ++i)
00884             mapBx[i] = 0.0;
00885     }
00886     if(!mapBy && By != 0.0) {
00887         mapBy = new float[nX*nY*nZ];
00888         assert(mapBy != 0);
00889         for(int i=0; i<nX*nY*nZ; ++i)
00890             mapBy[i] = 0.0;
00891     }
00892     if(!mapBz && Bz != 0.0) {
00893         mapBz = new float[nX*nY*nZ];
00894         assert(mapBz != 0);
00895         for(int i=0; i<nX*nY*nZ; ++i)
00896             mapBz[i] = 0.0;
00897     }
00898     if(!mapEx && Ex != 0.0) {
00899         mapEx = new float[nX*nY*nZ];
00900         assert(mapEx != 0);
00901         for(int i=0; i<nX*nY*nZ; ++i)
00902             mapEx[i] = 0.0;
00903     }
00904     if(!mapEy && Ey != 0.0) {
00905         mapEy = new float[nX*nY*nZ];
00906         assert(mapEy != 0);
00907         for(int i=0; i<nX*nY*nZ; ++i)
00908             mapEy[i] = 0.0;
00909     }
00910     if(!mapEz && Ez != 0.0) {
00911         mapEz = new float[nX*nY*nZ];
00912         assert(mapEz != 0);
00913         for(int i=0; i<nX*nY*nZ; ++i)
00914             mapEz[i] = 0.0;
00915     }
00916     X -= X0;
00917     Y -= Y0;
00918     Z -= Z0;
00919     int i = (int)floor((X/dX) + 0.5);
00920     if(i<0 || fabs(i*dX-X)>tolerance || i >= nX) {
00921         G4cerr << "MyBLFieldMap: ERROR point off grid X="
00922             << X << " line=" << linenumber << G4endl;
00923         return false;
00924     }
00925     int j = (int)floor((Y/dY) + 0.5);
00926     if(j<0 || fabs(j*dY-Y)>tolerance || j >= nY) {
00927         G4cerr << "MyBLFieldMap: ERROR point off grid Y="
00928             << Y << " line=" << linenumber << G4endl;
00929         return false;
00930     }
00931     int k = (int)floor((Z/dZ) + 0.5);
00932     if(k<0 || fabs(k*dZ-Z)>tolerance || k >= nZ) {
00933         G4cerr << "MyBLFieldMap: ERROR point off grid Z="
00934             << Z << " line=" << linenumber << G4endl;
00935         return false;
00936     }
00937     int m = k*nY*nX + j*nX + i;
00938     assert(m >= 0 && m < nX*nY*nZ);
00939     if(mapBx) mapBx[m] = Bx;
00940     if(mapBy) mapBy[m] = By;
00941     if(mapBz) mapBz[m] = Bz;
00942     if(mapEx) mapEx[m] = Ex;
00943     if(mapEy) mapEy[m] = Ey;
00944     if(mapEz) mapEz[m] = Ez;
00945 
00946     return true;
00947 }
00948 
00949 
00950 bool GridImpl::writeFile(FILE *f) {
00951     fprintf(f,"grid X0=%g Y0=%g Z0=%g nX=%d nY=%d nZ=%d dX=%g dY=%g dZ=%g\n",
00952         X0,Y0,Z0,nX,nY,nZ,dX,dY,dZ);
00953     if(extendX) fprintf(f,"extendX flip=%s\n",bits2str(extendXbits));
00954     if(extendY) fprintf(f,"extendY flip=%s\n",bits2str(extendYbits));
00955     if(extendZ) fprintf(f,"extendZ flip=%s\n",bits2str(extendZbits));
00956     fprintf(f,"data\n");
00957 
00958     for(int i=0; i<nX; ++i) {
00959         G4double X = X0 + i*dX;
00960         for(int j=0; j<nY; ++j) {
00961             G4double Y = Y0 + j*dY;
00962             for(int k=0; k<nZ; ++k) {
00963                 G4double Z = Z0 + k*dZ;
00964                 int m = k*nY*nX + j*nX + i;
00965                 assert(m >= 0 && m < nX*nY*nZ);
00966                 G4double Bx = (mapBx ? mapBx[m] : 0.0);
00967                 G4double By = (mapBy ? mapBy[m] : 0.0);
00968                 G4double Bz = (mapBz ? mapBz[m] : 0.0);
00969                 fprintf(f,"%.1f %.1f %.1f %g %g %g",
00970                     X,Y,Z,Bx/tesla,By/tesla,Bz/tesla);
00971                 if(hasE()) {
00972                     G4double Ex = (mapEx ? mapEx[m] : 0.0);
00973                     G4double Ey = (mapEy ? mapEy[m] : 0.0);
00974                     G4double Ez = (mapEz ? mapEz[m] : 0.0);
00975                     fprintf(f," %g %g %g",
00976                         Ex/(megavolt/meter),
00977                         Ey/(megavolt/meter),
00978                         Ez/(megavolt/meter));
00979                 }
00980                 fprintf(f,"\n");
00981             }
00982         }
00983     }
00984     return true;
00985 }
00986 
00987 
00988 CylinderImpl::CylinderImpl(MyBLArgumentVector &argv, MyBLArgumentMap &namedArgs)
00989 : FieldMapImpl() {
00990     nR = 2;
00991     nZ = 2;
00992     dR = 10.0*mm;
00993     dZ = 10.0*mm;
00994     Z0 = 0.0;
00995     tolerance = 0.01*mm;
00996     mapBr = 0;
00997     mapBz = 0;
00998     mapEr = 0;
00999     mapEz = 0;
01000     extendZ = false;
01001     extendBrFactor = extendBzFactor = 1.0;
01002     extendErFactor = extendEzFactor = 1.0;
01003     argInt(nR,"nR",namedArgs);
01004     argInt(nZ,"nZ",namedArgs);
01005     argDouble(dR,"dR",namedArgs);
01006     argDouble(dZ,"dZ",namedArgs);
01007     argDouble(Z0,"Z0",namedArgs);
01008     argDouble(tolerance,"tolerance",namedArgs);
01009 }
01010 
01011 
01012 CylinderImpl::~CylinderImpl() {
01013     if(mapBr) delete mapBr;
01014     if(mapBz) delete mapBz;
01015     if(mapEr) delete mapEr;
01016     if(mapEz) delete mapEz;
01017 }
01018 
01019 
01020 void CylinderImpl::getFieldValue(const G4double local[4], G4double field[6])
01021 const {
01022     G4double z = local[2];
01023     G4double r = sqrt(local[0]*local[0]+local[1]*local[1]);
01024 
01025     bool extending = false;
01026     if(z < 0.0 && extendZ) {
01027         z = -z;
01028         extending = true;
01029     }
01030 
01031     z -= Z0;
01032 
01033     // 2D linear average of the 4 surrounding points in the map
01034     int i = (int)floor(r/dR);
01035     int j = (int)floor(z/dZ);
01036     if(z < Z0 || i >= nR-1 || j < 0 || j >= nZ-1) {
01037         field[0] = field[1] = field[2] = field[3] = field[4] =
01038             field[5] = 0.0;
01039         return;
01040     }
01041 
01042     float fr = 1.0 - (r - i*dR) / dR;
01043     assert(fr >= 0.0 && fr <= 1.0);
01044     float fz = 1.0 - (z - j*dZ) / dZ;
01045     assert(fz >= 0.0 && fz <= 1.0);
01046 
01047     G4double Bz = 0.0;
01048     if(mapBz) Bz =
01049             mapBz[j*nR+i]*fr*fz +
01050             mapBz[j*nR+i+1]*(1.0-fr)*fz +
01051             mapBz[(j+1)*nR+i]*fr*(1.0-fz) +
01052             mapBz[(j+1)*nR+i+1]*(1.0-fr)*(1.0-fz);
01053     G4double Br = 0.0;
01054     if(mapBr) Br =
01055             mapBr[j*nR+i]*fr*fz +
01056             mapBr[j*nR+i+1]*(1.0-fr)*fz +
01057             mapBr[(j+1)*nR+i]*fr*(1.0-fz) +
01058             mapBr[(j+1)*nR+i+1]*(1.0-fr)*(1.0-fz);
01059     G4double Ez = 0.0;
01060     if(mapEz) Ez =
01061             mapEz[j*nR+i]*fr*fz +
01062             mapEz[j*nR+i+1]*(1.0-fr)*fz +
01063             mapEz[(j+1)*nR+i]*fr*(1.0-fz) +
01064             mapEz[(j+1)*nR+i+1]*(1.0-fr)*(1.0-fz);
01065     G4double Er = 0.0;
01066     if(mapEr) Er =
01067             mapEr[j*nR+i]*fr*fz +
01068             mapEr[j*nR+i+1]*(1.0-fr)*fz +
01069             mapEr[(j+1)*nR+i]*fr*(1.0-fz) +
01070             mapEr[(j+1)*nR+i+1]*(1.0-fr)*(1.0-fz);
01071     if(extending) {
01072         Br *= extendBrFactor;
01073         Bz *= extendBzFactor;
01074         Er *= extendErFactor;
01075         Ez *= extendEzFactor;
01076     }
01077 
01078     G4double phi = atan2(local[1], local[0]);
01079     field[0] = Br * cos(phi);
01080     field[1] = Br * sin(phi);
01081     field[2] = Bz;
01082     field[3] = Er * cos(phi);
01083     field[4] = Er * sin(phi);
01084     field[5] = Ez;
01085 }
01086 
01087 
01088 bool CylinderImpl::handleCommand(InputFile &in, MyBLArgumentVector &argv,
01089 MyBLArgumentMap &namedArgs) {
01090     if(argv[0] == "extendZ") {
01091         extendZ = true;
01092         const char *p = namedArgs["flip"].c_str();
01093         if(strstr(p,"Br")) extendBrFactor = -1.0;
01094         if(strstr(p,"Bz")) extendBzFactor = -1.0;
01095         if(strstr(p,"Er")) extendErFactor = -1.0;
01096         if(strstr(p,"Ez")) extendEzFactor = -1.0;
01097         return true;
01098     }
01099     else if(argv[0] == "Br") {
01100         if(mapBr) return false;
01101         mapBr = new float[nR*nZ];
01102         for(int i=0; i<nR*nZ; ++i) mapBr[i] = 0.0;
01103         return readBlock(in,mapBr,nZ,nR,tesla);
01104     }
01105     else if(argv[0] == "Bz") {
01106         if(mapBz) return false;
01107         mapBz = new float[nR*nZ];
01108         for(int i=0; i<nR*nZ; ++i) mapBz[i] = 0.0;
01109         return readBlock(in,mapBz,nZ,nR,tesla);
01110     }
01111     else if(argv[0] == "Er") {
01112         if(mapEr) return false;
01113         mapEr = new float[nR*nZ];
01114         for(int i=0; i<nR*nZ; ++i) mapEr[i] = 0.0;
01115         return readBlock(in,mapEr,nZ,nR,megavolt/meter);
01116     }
01117     else if(argv[0] == "Ez") {
01118         if(mapEz) return false;
01119         mapEz = new float[nR*nZ];
01120         for(int i=0; i<nR*nZ; ++i) mapEz[i] = 0.0;
01121         return readBlock(in,mapEz,nZ,nR,megavolt/meter);
01122     }
01123     else if(argv[0] == "data") {
01124         for(char *line=0; (line=in.getline())!=0;) {
01125             if(isalpha(line[0])) {
01126                 in.repeatLine();
01127                 break;
01128             }
01129             int n;
01130             float R=0.0,Z=0.0,Br=0.0,Bz=0.0,Er=0.0,Ez=0.0;
01131             for(char *p=line; (p=strchr(p,','))!=0;) *p = ' ';
01132             n = sscanf(line,"%f%f%f%f%f%f",&R,&Z,&Br,&Bz,&Er,&Ez);
01133             if(n <= 2) continue;
01134             setField(R,Z,Br*tesla,Bz*tesla,Er*megavolt/meter,
01135                 Ez*megavolt/meter,in.linenumber());
01136         }
01137         return true;
01138     }
01139     else {
01140         return false;
01141     }
01142 }
01143 
01144 
01145 void CylinderImpl::getBoundingPoint(int i, G4double point[4]) {
01146     point[0] = (i&1 ? 1.0 : -1.0) * (nR-1) * dR;
01147     point[1] = (i&2 ? 1.0 : -1.0) * (nR-1) * dR;
01148     if(extendZ)
01149         point[2] = (i&4 ? 1.0 : -1.0) * (nZ-1) * dZ;
01150     else
01151         point[2] = (i&4 ? Z0 : Z0+(nZ-1)*dZ);
01152 }
01153 
01154 
01155 bool CylinderImpl::setField(G4double R, G4double Z, G4double Br,
01156 G4double Bz, G4double Er, G4double Ez, int linenumber) {
01157     if((Br != 0.0 || Bz != 0.0) && !mapBr) {
01158         mapBr = new float[nR*nZ];
01159         mapBz = new float[nR*nZ];
01160         assert(mapBr != 0 && mapBz != 0);
01161         for(int i=0; i<nR*nZ; ++i)
01162             mapBr[i] = mapBz[i] = 0.0;
01163     }
01164     if((Er != 0.0 || Ez != 0.0) && !mapEr) {
01165         mapEr = new float[nR*nZ];
01166         mapEz = new float[nR*nZ];
01167         assert(mapEr != 0 && mapEz != 0);
01168         for(int i=0; i<nR*nZ; ++i)
01169             mapEr[i] = mapEz[i] = 0.0;
01170     }
01171     int i = (int)floor((R/dR) + 0.5);
01172     if(i<0 || fabs(i*dR-R)>tolerance || i >= nR) {
01173         G4cerr << "MyBLFieldMap: ERROR point off grid R="
01174             << R << " line=" << linenumber << G4endl;
01175         return false;
01176     }
01177     int j = (int)floor(((Z-Z0)/dZ) + 0.5);
01178     if(j<0 || fabs(j*dZ+Z0-Z)>tolerance || j >= nZ) {
01179         G4cerr << "MyBLFieldMap: ERROR point off grid Z="
01180             << Z << " line=" << linenumber << G4endl;
01181         return false;
01182     }
01183     if(mapBr) {
01184         mapBr[j*nR+i] = Br;
01185         mapBz[j*nR+i] = Bz;
01186     }
01187     if(mapEr) {
01188         mapEr[j*nR+i] = Er;
01189         mapEz[j*nR+i] = Ez;
01190     }
01191 
01192     return true;
01193 }
01194 
01195 
01196 bool CylinderImpl::writeFile(FILE *f) {
01197     fprintf(f,"cylinder Z0=%g nR=%d nZ=%d dR=%g dZ=%g\n",Z0,nR,nZ,dR,dZ);
01198     if(extendZ) {
01199         fprintf(f,"extendZ flip=");
01200         if(extendBrFactor < 0.0) fprintf(f,"Br,");
01201         if(extendBzFactor < 0.0) fprintf(f,"Bz,");
01202         if(extendErFactor < 0.0) fprintf(f,"Er,");
01203         if(extendEzFactor < 0.0) fprintf(f,"Ez,");
01204         fprintf(f,"\n");
01205     }
01206     fprintf(f,"data\n");
01207 
01208     for(int i=0; i<nR; ++i) {
01209         G4double R = i*dR;
01210         for(int j=0; j<nZ; ++j) {
01211             G4double Z = Z0 + j*dZ;
01212             int m = j*nR + i;
01213             assert(m >= 0 && m < nR*nZ);
01214             G4double Br = (mapBr ? mapBr[m] : 0.0);
01215             G4double Bz = (mapBz ? mapBz[m] : 0.0);
01216             fprintf(f,"%.1f %.1f %g %g", R,Z,Br/tesla,Bz/tesla);
01217             if(hasE()) {
01218                 G4double Er = (mapEr ? mapEr[m] : 0.0);
01219                 G4double Ez = (mapEz ? mapEz[m] : 0.0);
01220                 fprintf(f," %g %g", Er/(megavolt/meter),
01221                     Ez/(megavolt/meter));
01222             }
01223             fprintf(f,"\n");
01224         }
01225     }
01226     return true;
01227 }