// AK20110616 Modified from G4beamline v2.03 BLFieldMap.cc
//
/*
This source file is part of G4beamline, http://g4beamline.muonsinc.com
Copyright (C) 2003,2004,2005,2006 by Tom Roberts, all rights reserved.

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

http://www.gnu.org/copyleft/gpl.html
*/

#include <fstream>
#include <string.h>
#include <ctype.h>
#include <assert.h>

#include "IBLFieldMap.hxx"
#include "IBLFuncs.hxx"
#include "BTSpline1D.hxx"

#include <ICOMETLog.hxx>
#include "IFieldManager.hxx" 

namespace COMET {
/** class InputFile handles input file I/O
 **/
class InputFile {
    std::ifstream in;
    std::string name;
    char *line;
    int maxline;
    int lineno;
    bool repeatline;
    public:
        InputFile(std::string fname, Int_t _maxline=1024) : in(), name() {
            name = fname;
            in.open(name.c_str());
            maxline = _maxline;
            line = new char[maxline];
            assert(line != 0);
            lineno = 0;
            repeatline = false;
        }
        ~InputFile() { close(); }
        const char *filename() { return name.c_str(); }
        bool good() { return in.good(); }
        int linenumber() { return lineno; }
        char *getline() {                             // skips blank lines and comments
            if(repeatline) {
                repeatline = false;
                return line;
            }
            while(in.good()) {
                line[0] = '\0';
                in.getline(line,maxline);
                ++lineno;
                if(line[0] == '*')
                    printf("%s\n",line);
                if(line[0] == '\0' || line[0] == '#' || line[0] == '*')
                    continue;
                return line;
            }
            return 0;
        }
        void setMaxline(int _maxline) {
            maxline = _maxline;
            assert(maxline >= 80);
            if(line) delete[] line;
            line = new char[maxline];
            assert(line != 0);
        }
        void close() {
            if(line) {
                in.close(); delete[] line; line = 0;
            }
        }
        void repeatLine() { repeatline = true; }
};

/** class FieldMapImpl -- base class for a FieldMap implementation
 **/
class FieldMapImpl {
    public:
        FieldMapImpl() { }
        virtual ~FieldMapImpl() { }
        virtual void getFieldValue(const Double_t local[4], Double_t field[6])
            const = 0;
        virtual bool handleCommand(InputFile &in, MyBLArgumentVector &argv,
            MyBLArgumentMap &namedArgs) = 0;
        virtual void getBoundingPoint(int i, Double_t point[4]) = 0;
        virtual void getBoundingBox(BoundingBox_t& aBox) = 0;
        virtual bool hasB() = 0;
        virtual bool hasE() = 0;
        bool readBlock(InputFile &in, float *values, int nRows, int nCols,
            Double_t units);
        virtual bool writeFile(FILE *f) = 0;
};

/** class GridImpl -- class for a Grid FieldMap implementation
 **/
class GridImpl : public FieldMapImpl {
    Int_t nX;
    Int_t nY;
    Int_t nZ;
    Double_t dX;
    Double_t dY;
    Double_t dZ;
    Double_t X0;
    Double_t Y0;
    Double_t Z0;
    Double_t 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;
    public:
        GridImpl(MyBLArgumentVector &argv, MyBLArgumentMap &namedArgs);
        ~GridImpl();
        void getFieldValue(const Double_t local[4], Double_t field[6]) const;
        bool handleCommand(InputFile &in, MyBLArgumentVector &argv,
            MyBLArgumentMap &namedArgs);
        virtual void getBoundingPoint(int i, Double_t point[4]);
        virtual void getBoundingBox(BoundingBox_t& aBox);
        virtual bool hasB() { return mapBx!=0 || mapBy!=0 || mapBz!=0; }
        virtual bool hasE() { return mapEx!=0 || mapEy!=0 || mapEz!=0; }
        int bits(std::string s);
        const char *bits2str(int v);
        bool setField(Double_t X, Double_t Y, Double_t Z, Double_t Bx,
            Double_t By, Double_t Bz, Double_t Ex, Double_t Ey,
            Double_t Ez, int linenumber);
        virtual bool writeFile(FILE *f);
};

/** class CylinderImpl -- class for a Cylinder FieldMap implementation
 **/
class CylinderImpl : public FieldMapImpl {
    Int_t nR;
    Int_t nZ;
    Double_t dR;
    Double_t dZ;
    Double_t Z0;
    Double_t tolerance;
    float *mapBr;
    float *mapBz;
    float *mapEr;
    float *mapEz;
    // axesOrient used to define which axis is parallel to the cylinder and 
    // which are perpendiculat in the local coordinates of the cylinder.
    // This can be used to define a cylinder along x, y or z in global 
    // coordinates without having to rotate when entering the local coordinates.
    // Defaults to a cylinder oriented along Z, as is convention.
    int axesOrient[3];
    bool extendZ;
    float extendBrFactor, extendBzFactor;
    float extendErFactor, extendEzFactor;
    public:
        CylinderImpl(MyBLArgumentVector &argv, MyBLArgumentMap &namedArgs);
        ~CylinderImpl();
        void getFieldValue(const Double_t local[4], Double_t field[6]) const;
        bool handleCommand(InputFile &in, MyBLArgumentVector &argv,
            MyBLArgumentMap &namedArgs);
        virtual void getBoundingPoint(int i, Double_t point[4]);
        virtual void getBoundingBox(BoundingBox_t& aBox);
        virtual bool hasB() { return mapBz != 0 || mapBr != 0; }
        virtual bool hasE() { return mapEz != 0 || mapEr != 0; }
        bool setField(Double_t R, Double_t Z, Double_t Br,
            Double_t Bz, Double_t Er, Double_t Ez, int linenumber);
        virtual bool writeFile(FILE *f);
};

/** class TimeImpl -- class implementation of a time dependence.
 */
class TimeImpl {
    BTSpline1D B;
    BTSpline1D E;
    Double_t tMin;
    Double_t tMax;
    Double_t period;
    public:
        /// Constructor.
        /// n is # points in the arrays, t[] is the list of time values (ns),
        /// b[] and e[] are the factors for B field and E field.
        /// t[], b[], and e[] are copied so the caller can
        /// delete them after calling this routine.
        TimeImpl(int n, Double_t t[], Double_t b[], Double_t e[]=0) :
        B(n,t,b), E(n,t,(e?e:b)) {
            tMin = t[0];
            tMax = t[n-1];
            period = -1.0;
        }
        void setPeriod(Double_t v) { period = v; }
        Double_t factorB(Double_t t) {
            if(period > 0.0) {
                t -= floor(t/period)*period;
            }
            else {
                if(t < tMin) t = tMin;
                if(t > tMax) t = tMax;
            }
            return B(t);
        }
        Double_t factorE(Double_t t) {
            if(period > 0.0) {
                t -= floor(t/period)*period;
            }
            else {
                if(t < tMin) t = tMin;
                if(t > tMax) t = tMax;
            }
            return E(t);
        }
        static TimeImpl *readTime(InputFile &in, MyBLArgumentVector argv,
            MyBLArgumentMap namedArgs);
};


/// argDouble handles an argument with a value of Double_t
static void argDouble(Double_t &var, const char *name, MyBLArgumentMap &namedArgs) {
    if(namedArgs.count(name) == 0) return;
    char *q;
    const char *p=namedArgs[name].c_str();
    Double_t v = strtod(p,&q);
    if(p == q || *q != '\0') {
        COMETError("Invalid value for "<< name<<" in IBLFieldMap");
        return;
    }
    var = v;
}


/// argInt handles an argument with a value of Int_t
static void argInt(Int_t &var, const char *name, MyBLArgumentMap &namedArgs) {
    if(namedArgs.count(name) == 0) return;
    char *q;
    const char *p=namedArgs[name].c_str();
    Int_t v = strtol(p,&q,0);
    if(p == q || *q != '\0') {
      COMETError( "Invalid value for "<< name<<" in IBLFieldMap") ;
        return;
    }
    var = v;
}


/// d2string converts a double to a std::string
static std::string d2string(Double_t v) {
    char tmp[32];
    sprintf(tmp,"%.8g",v);
    return std::string(tmp);
}


/// i2string converts an int to a std::string
static std::string i2string(int v) {
    char tmp[32];
    sprintf(tmp,"%d",v);
    return std::string(tmp);
}


IBLFieldMap::IBLFieldMap() {
    maxline = 1024;
    current = 1.0;
    gradient = 1.0;
    normB = 1.0;
    normE = 1.0;
    impl = 0;
    time = 0;
}


IBLFieldMap::~IBLFieldMap() {
    if(impl) delete impl;
    impl = 0;
    if(time) delete time;
    time = 0;
}


void IBLFieldMap::getFieldValue(const Double_t local[4], Double_t field[6],
                                      Double_t _current, Double_t _gradient) {
    if(!impl)
        throw "IBLFieldMap::getFieldValue called, no implementation";

    Double_t thisField[6];
    impl->getFieldValue(local,thisField);
    Double_t timeB=1.0, timeE=1.0;
    if(time) {
        timeB = time->factorB(local[3]);
        timeE = time->factorE(local[3]);
    }
    field[0] = thisField[0] * normB * timeB * _current/current;
    field[1] = thisField[1] * normB * timeB * _current/current;
    field[2] = thisField[2] * normB * timeB * _current/current;
    field[3] = thisField[3] * normE * timeE * _gradient/gradient;
    field[4] = thisField[4] * normE * timeE * _gradient/gradient;
    field[5] = thisField[5] * normE * timeE * _gradient/gradient;
}


bool IBLFieldMap::readFile(std::string filename) {
    InputFile in(filename,maxline);
    if(!in.good()) {
        COMETError (std::string("IBLFieldMap::readFile: cannot open file ") + in.filename());
        return false;

    }
    printf("IBLFieldMap: reading file '%s'\n",in.filename());

    bool retval = true;

    char *line;
    while((line=in.getline()) != 0) {
        MyBLArgumentVector argv;
        MyBLArgumentMap namedArgs;
        if(BLFuncs.parseArgs(line,argv,namedArgs) < 0)
            goto invalid;
        if(argv[0] == "") continue;
        if(argv[0] == "param") {
            argInt(maxline,"maxline",namedArgs);
            argDouble(current,"current",namedArgs);
            argDouble(gradient,"gradient",namedArgs);
            argDouble(normB,"normB",namedArgs);
            argDouble(normE,"normE",namedArgs);
            in.setMaxline(maxline);
        }
        else if(argv[0] == "grid") {
            if(impl) goto invalid;
            impl = new GridImpl(argv,namedArgs);
        }
        else if(argv[0] == "cylinder") {
            if(impl) goto invalid;
            impl = new CylinderImpl(argv,namedArgs);
        }
        else if(argv[0] == "time") {
            if(time) goto invalid;
            time = TimeImpl::readTime(in,argv,namedArgs);
            if(!time) goto invalid;
        }
        else if(impl) {
            if(!impl->handleCommand(in,argv,namedArgs))
                goto invalid;
        }
        else {
            invalid:        
                COMETError("IBLFieldMap file "  << in.filename()
                           << " line " << in.linenumber()
                           << " invalid command " << argv[0].c_str());
            retval = false;
            break;
        }
    }
    in.close();

    return retval;
}


void IBLFieldMap::getBoundingPoint(int i, Double_t point[4]) {
    impl->getBoundingPoint(i,point);
}

void IBLFieldMap::getBoundingBox(BoundingBox_t& aBox) {
    impl->getBoundingBox(aBox);
}


bool IBLFieldMap::hasB() {
    return impl->hasB();
}


bool IBLFieldMap::hasE() {
    return impl->hasE();
}


TimeImpl *TimeImpl::readTime(InputFile &in, MyBLArgumentVector /*argv*/,
MyBLArgumentMap namedArgs) {
    Double_t period=-1.0;
    argDouble(period,"period",namedArgs);

    std::vector<Double_t> T;
    std::vector<Double_t> B;
    std::vector<Double_t> E;

    char *line;
    while((line=in.getline()) != 0) {
        if(isalpha(line[0])) {
            in.repeatLine();
            break;
        }
        Double_t t,b,e;
        int j = sscanf(line,"%lf%lf%lf",&t,&b,&e);
        if(j < 2 || j > 3) return 0;
        if(j == 2) e = b;
        T.push_back(t);
        B.push_back(b);
        E.push_back(e);
    }

    if(T.size() > 0) {
        TimeImpl *impl =
            new TimeImpl(T.size(),&T[0],&B[0],&E[0]);
        impl->setPeriod(period);
        return impl;
    }

    return 0;
}


bool FieldMapImpl::readBlock(InputFile &in, float *values, int nRows, int nCols,
Double_t units) {
    while(nRows-- > 0) {
        char *line = in.getline();
        if(!line) return false;
        char *p=line;
        for(int i=0; i<nCols; ++i) {
            while(isspace(*p)) ++p;
            if(*p == '\0') return false;
            *values++ = strtod(p,&p) * units;
            if(*p == ',') ++p;
        }
    }
    return true;
}


bool IBLFieldMap::createGridMap(Double_t X0, Double_t Y0, Double_t Z0,
Double_t dX, Double_t dY, Double_t dZ, int nX, int nY, int nZ,
IFieldManager *emField) {
    if(impl) {
        delete impl;
        impl = 0;
    }
    maxline = 128;
    current = 1.0;
    gradient = 1.0;
    normB = 1.0;
    normE = 1.0;

    MyBLArgumentVector argv;
    MyBLArgumentMap args;
    args["X0"] = d2string(X0);
    args["Y0"] = d2string(Y0);
    args["Z0"] = d2string(Z0);
    args["dX"] = d2string(dX);
    args["dY"] = d2string(dY);
    args["dZ"] = d2string(dZ);
    args["nX"] = i2string(nX);
    args["nY"] = i2string(nY);
    args["nZ"] = i2string(nZ);

    GridImpl *grid = new GridImpl(argv,args);
    impl = grid;

    bool retval = true;
    Double_t pos[4], field[6];
    pos[3] = 0.0;
    for(int i=0; i<nX; ++i) {
        pos[0] = X0 + i*dX;
        for(int j=0; j<nY; ++j) {
            pos[1] = Y0 + j*dY;
            for(int k=0; k<nZ; ++k) {
                pos[2] = Z0 + k*dZ;
                emField->GetFieldValue(pos,field);
                if(!grid->setField(pos[0],pos[1],pos[2],
                    field[0],field[1],field[2],
                    field[3],field[4],field[5],0))
                    retval = false;
            }
        }
    }

    return retval;
}


bool IBLFieldMap::createCylinderMap(Double_t Z0, Double_t dR,  Double_t dZ,
int nR, int nZ, IFieldManager *emField) {
    if(impl) {
        delete impl;
        impl = 0;
    }
    maxline = 128;
    current = 1.0;
    gradient = 1.0;
    normB = 1.0;
    normE = 1.0;

    MyBLArgumentVector argv;
    MyBLArgumentMap args;
    args["Z0"] = d2string(Z0);
    args["dR"] = d2string(dR);
    args["dZ"] = d2string(dZ);
    args["nR"] = i2string(nR);
    args["nZ"] = i2string(nZ);

    CylinderImpl *cyl = new CylinderImpl(argv,args);
    impl = cyl;

    // use the Y=0 plane for the R,Z plane
    bool retval = true;
    Double_t pos[4], field[6];
    pos[3] = 0.0;
    for(int i=0; i<nR; ++i) {
        pos[0] = 0.0 + i*dR;
        pos[1] = 0.0;
        for(int k=0; k<nZ; ++k) {
            pos[2] = Z0 + k*dZ;
            emField->GetFieldValue(pos,field);
            if(!cyl->setField(pos[0],pos[2],field[0],field[2],
                field[3],field[5],0))
                retval = false;
        }
    }

    return false;
}


bool IBLFieldMap::writeFile(std::string filename, std::string comment) {
    if(!impl) return false;

    FILE *f = fopen(filename.c_str(),"r");
    if(f) {
        fclose(f);
        COMETError("IBLFieldMap: Output File Exists "<< filename.c_str());
        abort();
    }

    f = fopen(filename.c_str(),"w");
    if(!f) {
        fprintf(stderr,"IBLFieldMap::writeFile CANNOT WRITE file '%s'\n",
            filename.c_str());
        return false;
    }

    fprintf(f,"# %s\n",comment.c_str());
    fprintf(f,"param maxline=%d current=%g gradient=%g normB=%g normE=%g\n",
        maxline,current,gradient,normB,normE);

    bool retval = impl->writeFile(f);

    fclose(f);
    return retval;
}


bool IBLFieldMap::createTimeDependence(int n, Double_t t[], Double_t b[],
Double_t e[], Double_t period) {
    if(time) return false;
    time = new TimeImpl(n,t,b,e);
    if(period > 0.0)
        time->setPeriod(period);
    return true;
}


bool IBLFieldMap::getTimeFactor(Double_t t, Double_t *b, Double_t *e) {
    if(b) *b = (time ? time->factorB(t) : 1.0);
    if(e) *e = (time ? time->factorE(t) : 1.0);
    return true;
}


GridImpl::GridImpl(MyBLArgumentVector &/*argv*/, MyBLArgumentMap &namedArgs)
: FieldMapImpl() {
    nX = 2;
    nY = 2;
    nZ = 2;
    dX = 10.0*unit::mm;
    dY = 10.0*unit::mm;
    dZ = 10.0*unit::mm;
    X0 = 0.0;
    Y0 = 0.0;
    Z0 = 0.0;
    tolerance = 0.01*unit::mm;
    mapBx = 0;
    mapBy = 0;
    mapBz = 0;
    mapEx = 0;
    mapEy = 0;
    mapEz = 0;
    extendX = false;
    extendY = false;
    extendZ = false;
    extendXbits = 0;
    extendYbits = 0;
    extendZbits = 0;
    argInt(nX,"nX",namedArgs);
    argInt(nY,"nY",namedArgs);
    argInt(nZ,"nZ",namedArgs);
    argDouble(dX,"dX",namedArgs);
    argDouble(dY,"dY",namedArgs);
    argDouble(dZ,"dZ",namedArgs);
    argDouble(X0,"X0",namedArgs);
    argDouble(Y0,"Y0",namedArgs);
    argDouble(Z0,"Z0",namedArgs);
    argDouble(tolerance,"tolerance",namedArgs);
}


GridImpl::~GridImpl() {
    if(mapBx) delete[] mapBx;
    if(mapBy) delete[] mapBy;
    if(mapBz) delete[] mapBz;
    if(mapEx) delete[] mapEx;
    if(mapEy) delete[] mapEy;
    if(mapEz) delete[] mapEz;
}


void GridImpl::getFieldValue(const Double_t local[4], Double_t field[6])
const {
    Double_t x = local[0] - X0;
    Double_t y = local[1] - Y0;
    Double_t z = local[2] - Z0;

    // TODO put more guards on this to slim it down a bit
    Double_t factor[6];
    factor[0]=factor[1]=factor[2]=factor[3]=factor[4]=factor[5]=1.0;
    if(extendX && x < 0.0) {
        x = -x;
        for(int i=0; i<6; ++i) {
            if(extendXbits & (1<<i)) factor[i] = -factor[i];
        }
    }
    if(extendY && y < 0.0) {
        y = -y;
        for(int i=0; i<6; ++i) {
            if(extendYbits & (1<<i)) factor[i] = -factor[i];
        }
    }
    if(extendZ && z < 0.0) {
        z = -z;
        for(int i=0; i<6; ++i) {
            if(extendZbits & (1<<i)) factor[i] = -factor[i];
        }
    }

    // We compute a 3D linear average of the 8 surrounding points in the map
    // First, get the X,Y,Z indices into i,j,k
    int i = (int)floor(x/dX);
    int j = (int)floor(y/dY);
    int k = (int)floor(z/dZ);
    if(i < 0 || i >= nX-1 || j < 0 || j >= nY-1 || k < 0 || k >= nZ-1) {
        field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0;
        return;
    }
    // index is the initial index (corner of the cube with minimum X, Y, and Z)
    int index = k*nY*nX + j*nX + i;
    assert(index+nY*nX+nX+1 < nX*nY*nZ);

    // now compute the fractional weighting factors for X, Y, and Z
    float fx = 1.0 - (x - i*dX) / dX;
    assert(fx >= 0.0 && fx <= 1.0);
    float fy = 1.0 - (y - j*dY) / dY;
    assert(fy >= 0.0 && fy <= 1.0);
    float fz = 1.0 - (z - k*dZ) / dZ;
    assert(fz >= 0.0 && fz <= 1.0);

    // now compute the fractional weighting factors for the 8 corners
    float f0 = fx*fy*fz;
    float f1 = (1.0-fx)*fy*fz;
    float f2 = fx*(1.0-fy)*fz;
    float f3 = (1.0-fx)*(1.0-fy)*fz;
    float f4 = fx*fy*(1.0-fz);
    float f5 = (1.0-fx)*fy*(1.0-fz);
    float f6 = fx*(1.0-fy)*(1.0-fz);
    float f7 = (1.0-fx)*(1.0-fy)*(1.0-fz);

    // Finally, compute the components of the field
#define COMPONENT(C) \
        Double_t C = 0.0; \
        if(map##C) C =  map##C[index]*f0 + map##C[index+1]*f1 + \
                map##C[index+nX]*f2 + map##C[index+nX+1]*f3 + \
                map##C[index+nY*nX]*f4 + map##C[index+nY*nX+1]*f5 + \
                map##C[index+nY*nX+nX]*f6 + map##C[index+nY*nX+nX+1]*f7;
        COMPONENT(Bx);
    COMPONENT(By);
    COMPONENT(Bz);
    COMPONENT(Ex);
    COMPONENT(Ey);
    COMPONENT(Ez);

#ifdef STUB
    Double_t Bx = 0.0;
    if(mapBx) Bx =
            mapBx[index]*fx*fy*fz +
            mapBx[index+1]*(1.0-fx)*fy*fz +
            mapBx[index+nX]*fx*(1.0-fy)*fz +
            mapBx[index+nX+1]*(1.0-fx)*(1.0-fy)*fz +
            mapBx[index+nY*nX]*fx*fy*(1.0-fz) +
            mapBx[index+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
            mapBx[index+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
            mapBx[index+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
    Double_t By = 0.0;
    if(mapBy) By =
            mapBy[index]*fx*fy*fz +
            mapBy[index+1]*(1.0-fx)*fy*fz +
            mapBy[index+nX]*fx*(1.0-fy)*fz +
            mapBy[index+nX+1]*(1.0-fx)*(1.0-fy)*fz +
            mapBy[index+nY*nX]*fx*fy*(1.0-fz) +
            mapBy[index+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
            mapBy[index+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
            mapBy[index+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
    Double_t Bz = 0.0;
    if(mapBz) Bz =
            mapBz[index]*fx*fy*fz +
            mapBz[index+1]*(1.0-fx)*fy*fz +
            mapBz[index+nX]*fx*(1.0-fy)*fz +
            mapBz[index+nX+1]*(1.0-fx)*(1.0-fy)*fz +
            mapBz[index+nY*nX]*fx*fy*(1.0-fz) +
            mapBz[index+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
            mapBz[index+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
            mapBz[index+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);

    Double_t Ex = 0.0;
    if(mapEx) Ex =
            mapEx[index]*fx*fy*fz +
            mapEx[index+1]*(1.0-fx)*fy*fz +
            mapEx[index+nX]*fx*(1.0-fy)*fz +
            mapEx[index+nX+1]*(1.0-fx)*(1.0-fy)*fz +
            mapEx[index+nY*nX]*fx*fy*(1.0-fz) +
            mapEx[index+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
            mapEx[index+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
            mapEx[index+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
    Double_t Ey = 0.0;
    if(mapEy) Ey =
            mapEy[index]*fx*fy*fz +
            mapEy[index+1]*(1.0-fx)*fy*fz +
            mapEy[index+nX]*fx*(1.0-fy)*fz +
            mapEy[index+nX+1]*(1.0-fx)*(1.0-fy)*fz +
            mapEy[index+nY*nX]*fx*fy*(1.0-fz) +
            mapEy[index+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
            mapEy[index+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
            mapEy[index+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
    Double_t Ez = 0.0;
    if(mapEz) Ez =
            mapEz[index]*fx*fy*fz +
            mapEz[index+1]*(1.0-fx)*fy*fz +
            mapEz[index+nX]*fx*(1.0-fy)*fz +
            mapEz[index+nX+1]*(1.0-fx)*(1.0-fy)*fz +
            mapEz[index+nY*nX]*fx*fy*(1.0-fz) +
            mapEz[index+nY*nX+1]*(1.0-fx)*fy*(1.0-fz) +
            mapEz[index+nY*nX+nX]*fx*(1.0-fy)*(1.0-fz) +
            mapEz[index+nY*nX+nX+1]*(1.0-fx)*(1.0-fy)*(1.0-fz);
#endif                                            //STUB

    field[0] = Bx * factor[0];
    field[1] = By * factor[1];
    field[2] = Bz * factor[2];
    field[3] = Ex * factor[3];
    field[4] = Ey * factor[4];
    field[5] = Ez * factor[5];
}


bool GridImpl::handleCommand(InputFile &in, MyBLArgumentVector &argv,
MyBLArgumentMap &namedArgs) {
    if(argv[0] == "extendX") {
        extendX = true;
        extendXbits = bits(namedArgs["flip"]);
        return true;
    }
    else if(argv[0] == "extendY") {
        extendY = true;
        extendYbits = bits(namedArgs["flip"]);
        return true;
    }
    else if(argv[0] == "extendZ") {
        extendZ = true;
        extendZbits = bits(namedArgs["flip"]);
        return true;
    }
    else if(argv[0] == "Bx") {
      COMETError( "IBLFieldMap: Bx not implemented") ;
        return true;
    }
    else if(argv[0] == "By") {
      COMETError( "IBLFieldMap: By not implemented") ;
        return true;
    }
    else if(argv[0] == "Bz") {
      COMETError( "IBLFieldMap: Bz not implemented") ;
        return true;
    }
    else if(argv[0] == "Ex") {
      COMETError( "IBLFieldMap: Ex not implemented") ;
        return true;
    }
    else if(argv[0] == "Ey") {
      COMETError( "IBLFieldMap: Ey not implemented");
        return true;
    }
    else if(argv[0] == "Ez") {
        COMETError( "IBLFieldMap: Ez not implemented");
        return true;
    }
    else if(argv[0] == "data") {
        for(char *line=0; (line=in.getline())!=0; ) {
            if(isalpha(line[0])) {
                in.repeatLine();
                break;
            }
            int n;
            float X=0.0,Y=0.0,Z=0.0,Bx=0.0,By=0.0,Bz=0.0,
                Ex=0.0,Ey=0.0,Ez=0.0;
            for(char *p=line; (p=strchr(p,','))!=0;) *p = ' ';
            n = sscanf(line,"%f%f%f%f%f%f%f%f%f",&X,&Y,&Z,
                &Bx,&By,&Bz,&Ex,&Ey,&Ez);
            if(n <= 3) {
                COMETError("IBLFieldMap: invalid line " << in.linenumber());
                continue;
            }
            setField(X,Y,Z,Bx*unit::tesla,By*unit::tesla,Bz*unit::tesla,
                Ex*unit::megavolt/unit::meter,Ey*unit::megavolt/unit::meter,
                Ez*unit::megavolt/unit::meter,in.linenumber());
        }
        return true;
    }
    else {
        return false;
    }
}


void GridImpl::getBoundingPoint(int i, Double_t point[4]) {
    if(extendX)
        point[0] = (i&1 ? -(nX-1)*dX : (nX-1)*dX) + X0;
    else
        point[0] = (i&1 ? 0.0 : (nX-1)*dX) + X0;
    if(extendY)
        point[1] = (i&2 ? -(nY-1)*dY : (nY-1)*dY) + Y0;
    else
        point[1] = (i&2 ? 0.0 : (nY-1)*dY) + Y0;
    if(extendZ)
        point[2] = (i&4 ? -(nZ-1)*dZ : (nZ-1)*dZ) + Z0;
    else
        point[2] = (i&4 ? 0.0 : (nZ-1)*dZ) + Z0;
}

void GridImpl::getBoundingBox(BoundingBox_t& aBox) {
    if(extendX)
      aBox.MinX = X0 -(nX-1)*dX;
    else
      aBox.MinX =  X0;

    if(extendY)
      aBox.MinY = Y0 -(nY-1)*dY;
    else
      aBox.MinY =  Y0;

    if(extendZ)
      aBox.MinZ = Z0 -(nZ-1)*dZ;
    else
      aBox.MinZ =  Z0;
 
    aBox.MaxX= ((nX-1)*dX) + X0;
    aBox.MaxY= ((nY-1)*dY) + Y0;
    aBox.MaxZ= ((nZ-1)*dZ) + Z0;
}

int GridImpl::bits(std::string input) {
    int v=0;
    if(input.find("Bx") < input.size()) v |= 1;
    if(input.find("By") < input.size()) v |= 2;
    if(input.find("Bz") < input.size()) v |= 4;
    if(input.find("Ex") < input.size()) v |= 8;
    if(input.find("Ey") < input.size()) v |= 16;
    if(input.find("Ez") < input.size()) v |= 32;
    return v;
}


const char *GridImpl::bits2str(int v) {
    static char retval[32];
    retval[0] = '\0';
    if(v & 1) strcat(retval,"Bx,");
    if(v & 2) strcat(retval,"By,");
    if(v & 4) strcat(retval,"Bz,");
    if(v & 8) strcat(retval,"Ex,");
    if(v & 16) strcat(retval,"Ey,");
    if(v & 32) strcat(retval,"Ez,");
    return retval;
}


bool GridImpl::setField(Double_t X,  Double_t Y,  Double_t Z, 
                        Double_t Bx, Double_t By, Double_t Bz, 
                        Double_t Ex, Double_t Ey, Double_t Ez, int linenumber){
    if(!mapBx && Bx != 0.0) {
        mapBx = new float[nX*nY*nZ];
        assert(mapBx != 0);
        for(int i=0; i<nX*nY*nZ; ++i)
            mapBx[i] = 0.0;
    }
    if(!mapBy && By != 0.0) {
        mapBy = new float[nX*nY*nZ];
        assert(mapBy != 0);
        for(int i=0; i<nX*nY*nZ; ++i)
            mapBy[i] = 0.0;
    }
    if(!mapBz && Bz != 0.0) {
        mapBz = new float[nX*nY*nZ];
        assert(mapBz != 0);
        for(int i=0; i<nX*nY*nZ; ++i)
            mapBz[i] = 0.0;
    }
    if(!mapEx && Ex != 0.0) {
        mapEx = new float[nX*nY*nZ];
        assert(mapEx != 0);
        for(int i=0; i<nX*nY*nZ; ++i)
            mapEx[i] = 0.0;
    }
    if(!mapEy && Ey != 0.0) {
        mapEy = new float[nX*nY*nZ];
        assert(mapEy != 0);
        for(int i=0; i<nX*nY*nZ; ++i)
            mapEy[i] = 0.0;
    }
    if(!mapEz && Ez != 0.0) {
        mapEz = new float[nX*nY*nZ];
        assert(mapEz != 0);
        for(int i=0; i<nX*nY*nZ; ++i)
            mapEz[i] = 0.0;
    }
    X -= X0;
    Y -= Y0;
    Z -= Z0;

    int i = (int)floor((X/dX) + 0.5);
    int j = (int)floor((Y/dY) + 0.5);
    int k = (int)floor((Z/dZ) + 0.5);

    std::string wrongcoords = "";
    bool coorderror = false;
   
    if(i<0 || fabs(i*dX-X)>tolerance || i >= nX) {
        wrongcoords += " X ";
        coorderror = true;
    }
    if(j<0 || fabs(j*dY-Y)>tolerance || j >= nY) {
        wrongcoords += " Y ";
        coorderror = true;
    }
    if(k<0 || fabs(k*dZ-Z)>tolerance || k >= nZ) {
        wrongcoords += " Z ";
        coorderror = true;
    }
    if(coorderror){
        COMETError( "IBLFieldMap: ERROR coordinates off grid " << "(" << wrongcoords << "):"
                     << " X=" << X << " Y=" << Y << " Z=" << Z << ", line=" << linenumber);
        return false;
    }

    int index = k*nY*nX + j*nX + i;
    assert(index >= 0 && index < nX*nY*nZ);
    if(mapBx) mapBx[index] = Bx;
    if(mapBy) mapBy[index] = By;
    if(mapBz) mapBz[index] = Bz;
    if(mapEx) mapEx[index] = Ex;
    if(mapEy) mapEy[index] = Ey;
    if(mapEz) mapEz[index] = Ez;

    return true;
}


bool GridImpl::writeFile(FILE *f) {
    fprintf(f,"grid X0=%g Y0=%g Z0=%g nX=%d nY=%d nZ=%d dX=%g dY=%g dZ=%g\n",
        X0,Y0,Z0,nX,nY,nZ,dX,dY,dZ);
    if(extendX) fprintf(f,"extendX flip=%s\n",bits2str(extendXbits));
    if(extendY) fprintf(f,"extendY flip=%s\n",bits2str(extendYbits));
    if(extendZ) fprintf(f,"extendZ flip=%s\n",bits2str(extendZbits));
    fprintf(f,"data\n");

    for(int i=0; i<nX; ++i) {
        Double_t X = X0 + i*dX;
        for(int j=0; j<nY; ++j) {
            Double_t Y = Y0 + j*dY;
            for(int k=0; k<nZ; ++k) {
                Double_t Z = Z0 + k*dZ;
                int index = k*nY*nX + j*nX + i;
                assert(index >= 0 && index < nX*nY*nZ);
                Double_t Bx = (mapBx ? mapBx[index] : 0.0);
                Double_t By = (mapBy ? mapBy[index] : 0.0);
                Double_t Bz = (mapBz ? mapBz[index] : 0.0);
                fprintf(f,"%.1f %.1f %.1f %g %g %g",
                    X,Y,Z,Bx/unit::tesla,By/unit::tesla,Bz/unit::tesla);
                if(hasE()) {
                    Double_t Ex = (mapEx ? mapEx[index] : 0.0);
                    Double_t Ey = (mapEy ? mapEy[index] : 0.0);
                    Double_t Ez = (mapEz ? mapEz[index] : 0.0);
                    fprintf(f," %g %g %g",
                        Ex/(unit::megavolt/unit::meter),
                        Ey/(unit::megavolt/unit::meter),
                        Ez/(unit::megavolt/unit::meter));
                }
                fprintf(f,"\n");
            }
        }
    }
    return true;
}


CylinderImpl::CylinderImpl(MyBLArgumentVector &/*argv*/, MyBLArgumentMap &namedArgs)
: FieldMapImpl() {
    nR = 2;
    nZ = 2;
    dR = 10.0*unit::mm;
    dZ = 10.0*unit::mm;
    Z0 = 0.0;
    tolerance = 0.01*unit::mm;
    mapBr = 0;
    mapBz = 0;
    mapEr = 0;
    mapEz = 0;
    extendZ = false;
    extendBrFactor = extendBzFactor = 1.0;
    extendErFactor = extendEzFactor = 1.0;
    argInt(nR,"nR",namedArgs);
    argInt(nZ,"nZ",namedArgs);
    argDouble(dR,"dR",namedArgs);
    argDouble(dZ,"dZ",namedArgs);
    argDouble(Z0,"Z0",namedArgs);
    argDouble(tolerance,"tolerance",namedArgs);
    // Default assigns no intrinisc rotation
    for (int dim = 0; dim < 3; dim++) axesOrient[dim] = dim; 
}


CylinderImpl::~CylinderImpl() {
    if(mapBr) delete[] mapBr;
    if(mapBz) delete[] mapBz;
    if(mapEr) delete[] mapEr;
    if(mapEz) delete[] mapEz;
}


void CylinderImpl::getFieldValue(const Double_t local[4], Double_t field[6])
const {
    Double_t z = local[axesOrient[2]];
    Double_t r = sqrt(local[axesOrient[0]]*local[axesOrient[0]]+
                      local[axesOrient[1]]*local[axesOrient[1]]);
    z -= Z0;

    bool extending = false;
    if(z < 0.0 && extendZ) {
        z = -z;
        extending = true;
    }


    // 2D linear average of the 4 surrounding points in the map
    int i = (int)floor(r/dR);
    int j = (int)floor(z/dZ);
//    if(z < Z0 || i >= nR-1 || j < 0 || j >= nZ-1) {
    if(z < 0 || i >= nR-1 || j < 0 || j >= nZ-1) { // z < Z0 maybe wrong ????
        field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0;
        return;
    }

    float fr = 1.0 - (r - i*dR) / dR;
    assert(fr >= 0.0 && fr <= 1.0);
    float fz = 1.0 - (z - j*dZ) / dZ;
    assert(fz >= 0.0 && fz <= 1.0);

    Double_t Bz = 0.0;
    if(mapBz) Bz =
            mapBz[j*nR+i]*fr*fz +
            mapBz[j*nR+i+1]*(1.0-fr)*fz +
            mapBz[(j+1)*nR+i]*fr*(1.0-fz) +
            mapBz[(j+1)*nR+i+1]*(1.0-fr)*(1.0-fz);
    Double_t Br = 0.0;
    if(mapBr) Br =
            mapBr[j*nR+i]*fr*fz +
            mapBr[j*nR+i+1]*(1.0-fr)*fz +
            mapBr[(j+1)*nR+i]*fr*(1.0-fz) +
            mapBr[(j+1)*nR+i+1]*(1.0-fr)*(1.0-fz);
    Double_t Ez = 0.0;
    if(mapEz) Ez =
            mapEz[j*nR+i]*fr*fz +
            mapEz[j*nR+i+1]*(1.0-fr)*fz +
            mapEz[(j+1)*nR+i]*fr*(1.0-fz) +
            mapEz[(j+1)*nR+i+1]*(1.0-fr)*(1.0-fz);
    Double_t Er = 0.0;
    if(mapEr) Er =
            mapEr[j*nR+i]*fr*fz +
            mapEr[j*nR+i+1]*(1.0-fr)*fz +
            mapEr[(j+1)*nR+i]*fr*(1.0-fz) +
            mapEr[(j+1)*nR+i+1]*(1.0-fr)*(1.0-fz);
    if(extending) {
        Br *= extendBrFactor;
        Bz *= extendBzFactor;
        Er *= extendErFactor;
        Ez *= extendEzFactor;
    }

    Double_t x_dir = r!=0? local[axesOrient[0]]/r:0;
    Double_t y_dir = r!=0? local[axesOrient[1]]/r:0;
    field[axesOrient[0]] = Br * x_dir;
    field[axesOrient[1]] = Br * y_dir;
    field[axesOrient[2]] = Bz;
    field[axesOrient[0]+3] = Er * x_dir;
    field[axesOrient[1]+3] = Er * y_dir;
    field[axesOrient[2]+3] = Ez;
}


bool CylinderImpl::handleCommand(InputFile &in, MyBLArgumentVector &argv,
MyBLArgumentMap &namedArgs) {
    if(argv[0] == "orientZ") {
        const char *p = namedArgs["central"].c_str();
        // Orient cylinder along Z
        if(strstr(p,"Z")){
            // Local X perpendicular to central axis
            axesOrient[0] = 0;
            // Local Y perpendicular to central axis
            axesOrient[1] = 1;
            // Local Z parallel to central axis
            axesOrient[2] = 2;
        }
        // Orient cylinder along Y
        if(strstr(p,"Y")){
            // Local Z perpendicular to central axis
            axesOrient[0] = 2;
            // Local X perpendicular to central axis
            axesOrient[1] = 0;
            // Local Y parallel to central axis
            axesOrient[2] = 1;
        }
        // Orient cylinder along Y
        if(strstr(p,"X")){
            // Local Y perpendicular to central axis
            axesOrient[0] = 1;
            // Local Z perpendicular to central axis
            axesOrient[1] = 2;
            // Local X parallel to central axis
            axesOrient[2] = 0;
        }
        return true;
    }
    else if(argv[0] == "extendZ") {
        extendZ = true;
        const char *p = namedArgs["flip"].c_str();
        if(strstr(p,"Br")) extendBrFactor = -1.0;
        if(strstr(p,"Bz")) extendBzFactor = -1.0;
        if(strstr(p,"Er")) extendErFactor = -1.0;
        if(strstr(p,"Ez")) extendEzFactor = -1.0;
        return true;
    }
    else if(argv[0] == "Br") {
        if(mapBr) return false;
        mapBr = new float[nR*nZ];
        for(int i=0; i<nR*nZ; ++i) mapBr[i] = 0.0;
        return readBlock(in,mapBr,nZ,nR,unit::tesla);
    }
    else if(argv[0] == "Bz") {
        if(mapBz) return false;
        mapBz = new float[nR*nZ];
        for(int i=0; i<nR*nZ; ++i) mapBz[i] = 0.0;
        return readBlock(in,mapBz,nZ,nR,unit::tesla);
    }
    else if(argv[0] == "Er") {
        if(mapEr) return false;
        mapEr = new float[nR*nZ];
        for(int i=0; i<nR*nZ; ++i) mapEr[i] = 0.0;
        return readBlock(in,mapEr,nZ,nR,unit::megavolt/unit::meter);
    }
    else if(argv[0] == "Ez") {
        if(mapEz) return false;
        mapEz = new float[nR*nZ];
        for(int i=0; i<nR*nZ; ++i) mapEz[i] = 0.0;
        return readBlock(in,mapEz,nZ,nR,unit::megavolt/unit::meter);
    }
    else if(argv[0] == "data") {
        for(char *line=0; (line=in.getline())!=0;) {
            if(isalpha(line[0])) {
                in.repeatLine();
                break;
            }
            int n;
            float R=0.0,Z=0.0,Br=0.0,Bz=0.0,Er=0.0,Ez=0.0;
            for(char *p=line; (p=strchr(p,','))!=0;) *p = ' ';
            n = sscanf(line,"%f%f%f%f%f%f",&R,&Z,&Br,&Bz,&Er,&Ez);
            if(n <= 2) continue;
            setField(R,Z,Br*unit::tesla,Bz*unit::tesla,Er*unit::megavolt/unit::meter,
                Ez*unit::megavolt/unit::meter,in.linenumber());
        }
        return true;
    }
    else {
        return false;
    }
}


void CylinderImpl::getBoundingPoint(int i, Double_t point[4]) {
    point[axesOrient[0]] = (i&1 ? 1.0 : -1.0) * (nR-1) * dR;
    point[axesOrient[0]] = (i&2 ? 1.0 : -1.0) * (nR-1) * dR;
    if(extendZ)
// AK20150315: Bug in original G4beamline code, Z0 not included.
//        point[2] = (i&4 ? 1.0 : -1.0) * (nZ-1) * dZ;
        point[axesOrient[0]] = ((i&4 ? 1.0 : -1.0) * (nZ-1) * dZ) + Z0;
    else
        point[axesOrient[0]] = (i&4 ? Z0 : Z0+(nZ-1)*dZ);
}


void CylinderImpl::getBoundingBox(BoundingBox_t& aBox) {
    // Get the central axis for this cylindrical map
    int centralAxis = axesOrient[2];
    // Ensure this is a sensible argument
    assert(centralAxis == 0 || centralAxis == 1 || centralAxis == 2);
    // Conventional orientation with cylinder parallel to +Z
    if (centralAxis == 2){
        // Set the maximum values
        aBox.MaxX = aBox.MaxY = (nR-1) * dR;
        aBox.MaxZ = Z0 + (nZ-1)*dZ;
        // Set the minimum values
        aBox.MinX = aBox.MinY = -(nR-1) * dR;
        aBox.MinZ = Z0;
        // Extend Z if we need to
        if(extendZ) aBox.MinZ -= (nZ-1) * dZ;
    }
    // Cylinder is now parallel to Y, X and Z are perpendicular
    else if (centralAxis == 1){
        // Set the maximum values
        aBox.MaxZ = aBox.MaxX = (nR-1) * dR;
        aBox.MaxY = Z0 + (nZ-1)*dZ;
        // Set the minimum values
        aBox.MinZ = aBox.MinX = -(nR-1) * dR;
        aBox.MinY = Z0;
        // Extend Z if we need to
        if(extendZ) aBox.MinY -= (nZ-1) * dZ;
    }
    // Cylinder is now parallel to X, Y and Z are perpendicular
    else if (axesOrient[2] == 0){
        // Set the maximum values
        aBox.MaxY = aBox.MaxZ = (nR-1) * dR;
        aBox.MaxX = Z0 + (nZ-1)*dZ;
        // Set the minimum values
        aBox.MinY = aBox.MinZ = -(nR-1) * dR;
        aBox.MinX = Z0;
        // Extend Z if we need to
        if(extendZ) aBox.MinX -= (nZ-1) * dZ;
    }
}


bool CylinderImpl::setField(Double_t R, Double_t Z, Double_t Br,
Double_t Bz, Double_t Er, Double_t Ez, int linenumber) {
    if((Br != 0.0 || Bz != 0.0) && !mapBr) {
        mapBr = new float[nR*nZ];
        mapBz = new float[nR*nZ];
        assert(mapBr != 0 && mapBz != 0);
        for(int i=0; i<nR*nZ; ++i)
            mapBr[i] = mapBz[i] = 0.0;
    }
    if((Er != 0.0 || Ez != 0.0) && !mapEr) {
        mapEr = new float[nR*nZ];
        mapEz = new float[nR*nZ];
        assert(mapEr != 0 && mapEz != 0);
        for(int i=0; i<nR*nZ; ++i)
            mapEr[i] = mapEz[i] = 0.0;
    }

    int i = (int)floor((R/dR) + 0.5);
    int j = (int)floor(((Z-Z0)/dZ) + 0.5);
    
    std::string wrongcoords = "";
    bool coorderror = false;
    if(i<0 || fabs(i*dR-R)>tolerance || i >= nR) {
        wrongcoords += " R ";
        coorderror = true;
    }
    if(j<0 || fabs(j*dZ+Z0-Z)>tolerance || j >= nZ) {
        wrongcoords += " Z ";
        coorderror = true;
    }
    if(coorderror){
    	COMETError( "IBLFieldMap: ERROR coordinates off grid " << "(" << wrongcoords << "):"
                     << " R="<< R << " Z=" << Z << ", line=" << linenumber);
    	return false;
    }
    if(mapBr) {
        mapBr[j*nR+i] = Br;
        mapBz[j*nR+i] = Bz;
    }
    if(mapEr) {
        mapEr[j*nR+i] = Er;
        mapEz[j*nR+i] = Ez;
    }

    return true;
}


bool CylinderImpl::writeFile(FILE *f) {
    fprintf(f,"cylinder Z0=%g nR=%d nZ=%d dR=%g dZ=%g\n",Z0,nR,nZ,dR,dZ);
    if(extendZ) {
        fprintf(f,"extendZ flip=");
        if(extendBrFactor < 0.0) fprintf(f,"Br,");
        if(extendBzFactor < 0.0) fprintf(f,"Bz,");
        if(extendErFactor < 0.0) fprintf(f,"Er,");
        if(extendEzFactor < 0.0) fprintf(f,"Ez,");
        fprintf(f,"\n");
    }
    fprintf(f,"data\n");

    for(int i=0; i<nR; ++i) {
        Double_t R = i*dR;
        for(int j=0; j<nZ; ++j) {
            Double_t Z = Z0 + j*dZ;
            int index = j*nR + i;
            assert(index >= 0 && index < nR*nZ);
            Double_t Br = (mapBr ? mapBr[index] : 0.0);
            Double_t Bz = (mapBz ? mapBz[index] : 0.0);
            fprintf(f,"%.1f %.1f %g %g", R,Z,Br/unit::tesla,Bz/unit::tesla);
            if(hasE()) {
                Double_t Er = (mapEr ? mapEr[index] : 0.0);
                Double_t Ez = (mapEz ? mapEz[index] : 0.0);
                fprintf(f," %g %g", Er/(unit::megavolt/unit::meter),
                    Ez/(unit::megavolt/unit::meter));
            }
            fprintf(f,"\n");
        }
    }
    return true;
}
}