// ------------------------------------------------------------------------
// $RCSfile: Cyclotron.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Definitions for class: Cyclotron
//   Defines the abstract interface for a sector bend magnet.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2007/08/01 $
// $Author: Yang, Adelmann $
//
// ------------------------------------------------------------------------

#include "AbsBeamline/Cyclotron.h"
#include "AbsBeamline/BeamlineVisitor.h"
#include "Physics/Physics.h"
#include "Fields/Fieldmap.hh"
#include <fstream>

extern Inform *gmsg;

using Physics::pi;
using namespace std;

// Class Cyclotron
// ------------------------------------------------------------------------

Cyclotron::Cyclotron():
    Component() {
    Bfield.bfld = 0;
    Bfield.dbt = 0;
    Bfield.dbtt = 0;
    Bfield.dbttt = 0;
    BP.rarr = 0;
    Bfield.dbr = 0;
    Bfield.dbrr = 0;
    Bfield.dbrrr = 0;
    Bfield.dbrt = 0;
    Bfield.dbrrt = 0;
    Bfield.dbrtt = 0;
    Bfield.f2 = 0;
    Bfield.f3 = 0;
    Bfield.g3 = 0;
}


Cyclotron::Cyclotron(const Cyclotron &right):
    Component(right),
    fmapfn_m(right.fmapfn_m),
    rffrequ_m(right.rffrequ_m),
    rfphi_m(right.rfphi_m),
    escale_m(right.escale_m),
    symmetry_m(right.symmetry_m),
    rinit_m(right.rinit_m),
    prinit_m(right.prinit_m),
    phiinit_m(right.phiinit_m),
    type_m(right.type_m),
    harm_m(right.harm_m),
    bscale_m(right.bscale_m),
    tcr1_m(right.tcr1_m),
    tcr2_m(right.tcr2_m),
    mbtc_m(right.mbtc_m),
    slptc_m(right.slptc_m),
    superpose_m(right.superpose_m),
    RFfilename_m(right.RFfilename_m) {
    Bfield.bfld = 0;
    Bfield.dbt = 0;
    Bfield.dbtt = 0;
    Bfield.dbttt = 0;
    BP.rarr = 0;
    Bfield.dbr = 0;
    Bfield.dbrr = 0;
    Bfield.dbrrr = 0;
    Bfield.dbrt = 0;
    Bfield.dbrrt = 0;
    Bfield.dbrtt = 0;
    Bfield.f2 = 0;
    Bfield.f3 = 0;
    Bfield.g3 = 0;
}


Cyclotron::Cyclotron(const string &name):
    Component(name) {
    Bfield.bfld = 0;
    Bfield.dbt = 0;
    Bfield.dbtt = 0;
    Bfield.dbttt = 0;
    BP.rarr = 0;
    Bfield.dbr = 0;
    Bfield.dbrr = 0;
    Bfield.dbrrr = 0;
    Bfield.dbrt = 0;
    Bfield.dbrrt = 0;
    Bfield.dbrtt = 0;
    Bfield.f2 = 0;
    Bfield.f3 = 0;
    Bfield.g3 = 0;
}


Cyclotron::~Cyclotron() {

    if(Bfield.bfld) delete[] Bfield.bfld;
    if(Bfield.dbt) delete[] Bfield.dbt;
    if(Bfield.dbtt) delete[] Bfield.dbtt;
    if(Bfield.dbttt) delete[] Bfield.dbttt;

    if(BP.rarr) delete[] BP.rarr;

    if(Bfield.dbr) delete[] Bfield.dbr;
    if(Bfield.dbrr) delete[] Bfield.dbrr;
    if(Bfield.dbrrr) delete[] Bfield.dbrrr;

    if(Bfield.dbrt) delete[] Bfield.dbrt;
    if(Bfield.dbrrt) delete[] Bfield.dbrrt;
    if(Bfield.dbrtt) delete[] Bfield.dbrtt;

    if(Bfield.f2) delete[] Bfield.f2;
    if(Bfield.f3) delete[] Bfield.f3;
    if(Bfield.g3) delete[] Bfield.g3;
}


void Cyclotron::accept(BeamlineVisitor &visitor) const {
    visitor.visitCyclotron(*this);
}

void Cyclotron::setRinit(double rinit) {
    rinit_m = rinit;
}

double Cyclotron::getRinit() const {
    return rinit_m;
}


void Cyclotron::setPRinit(double prinit) {
    prinit_m = prinit;
}

double Cyclotron::getPRinit() const {
    return prinit_m;
}

void Cyclotron::setPHIinit(double phiinit) {
    phiinit_m = phiinit;
}

double Cyclotron::getPHIinit() const {
    return phiinit_m;
}

void Cyclotron::setFieldMapFN(string f) {
    fmapfn_m = f;
}

string Cyclotron::getFieldMapFN() const {
    return fmapfn_m;
}

void Cyclotron::setRfFieldMapFN(vector<string> f) {
    RFfilename_m = f;
}

void Cyclotron::setRfPhi(vector<double> f) {
    rfphi_m = f;
}

void Cyclotron::setRfFrequ(vector<double> f) {
    rffrequ_m = f;
}

void Cyclotron::setSuperpose(bool flag) {
    superpose_m = flag;
}

bool Cyclotron::getSuperpose() const {
    return superpose_m;
}

void Cyclotron::setSymmetry(double s) {
    symmetry_m = s;
}

double Cyclotron::getSymmetry() const {
    return symmetry_m;
}


void Cyclotron::setType(string t) {
    type_m = t;
}

const string &Cyclotron::getType() const {
    return type_m;
}

void Cyclotron::setCyclHarm(double h) {
    harm_m = h;
}

void Cyclotron::setBScale(double s) {
    bscale_m = s;
}

double Cyclotron::getBScale() const {
    return bscale_m;
}

void Cyclotron::setEScale(vector<double> s) {
    escale_m = s;
}

double Cyclotron::getCyclHarm() const {
    return harm_m;
}

double Cyclotron::getRmin() const {
    return BP.rmin;
}


double Cyclotron::getRmax() const {
    return BP.rmin + (Bfield.nrad - 1) * BP.delr;
}

// This function aims at obtaining magentic field at any given location R by interpolation.
// arguments t is useless here.
// arguments E is set to zero.


void Cyclotron::setTCr1(double tcr1) {
    tcr1_m = tcr1;
}

double Cyclotron::getTCr1() const {
    return tcr1_m;
}

void Cyclotron::setTCr2(double tcr2) {
    tcr2_m = tcr2;
}

double Cyclotron::getTCr2() const {
    return tcr2_m;
}
void Cyclotron::setMBtc(double mbtc) {
    mbtc_m = mbtc;
}

double Cyclotron::getMBtc() const {
    return mbtc_m;
}

void Cyclotron::setSLPtc(double slptc) {
    slptc_m = slptc;
}

double Cyclotron::getSLPtc() const {
    return slptc_m;
}

bool Cyclotron::apply(const int &i, const double &t, double E[], double B[]) {
    Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
    if(apply(RefPartBunch_m->R[i], Vector_t(0.0), t, Ev, Bv)) return true;

    E[0] = Ev(0);
    E[1] = Ev(1);
    E[2] = Ev(2);
    B[0] = Bv(0);
    B[1] = Bv(1);
    B[2] = Bv(2);

    return false;
}

bool Cyclotron::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
    return apply(RefPartBunch_m->R[i], Vector_t(0.0), t, E, B);
}

bool Cyclotron::apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) {

    const double rad = sqrt(R[0] * R[0] + R[1] * R[1]);
    const double xir = (rad - BP.rmin) / BP.delr;

    // ir : the mumber of path whoes radius is less then the 4 points of cell which surrond particle.
    const int    ir = (int)xir;

    // wr1 : the relative distance to the inner path radius
    const double wr1 = xir - (double)ir;
    // wr2 : the relative distance to the outer path radius
    const double wr2 = 1.0 - wr1;

    const double tempv = atan(R[1] / R[0]);
    double tet = tempv, tet_map, xit;

    /* if((R[0] > 0) && (R[1] >= 0)) tet = tempv;
       else*/
    if((R[0] < 0) && (R[1] >= 0)) tet = pi + tempv;
    else if((R[0] < 0) && (R[1] <= 0)) tet = pi + tempv;
    else if((R[0] > 0) && (R[1] <= 0)) tet = 2.0 * pi + tempv;
    else if((R[0] == 0) && (R[1] > 0)) tet = pi / 2.0;
    else if((R[0] == 0) && (R[1] < 0)) tet = 1.5 * pi;

    double tet_rad = tet;

    // the actual angle of particle
    tet = tet / pi * 180.0;

    // the corresponding angle on the field map
    // Note: this does not work if the start point of field map does not equal zero.
    tet_map = fmod(tet, 360.0 / symmetry_m);

    xit = tet_map / BP.dtet;

    int it = (int) xit;

    //    *gmsg << R << " tet_map= " << tet_map << " ir= " << ir << " it= " << it << " bf= " << Bfield.bfld[idx(ir,it)] << endl;

    const double wt1 = xit - (double)it;
    const double wt2 = 1.0 - wt1;

    // it : the number of point on the inner path whoes angle is less then the particle' corresponding angle.
    // include zero degree point
    it = it + 1;

    int r1t1, r2t1, r1t2, r2t2;
    int ntetS = Bfield.ntet + 1;

    // r1t1 : the index of the "min angle, min radius" point in the 2D field array.
    // considering  the array start with index of zero, minus 1.

    if(myBFieldType_m != FFAGBF) {
        /*
          For FFAG this does not work
        */
        r1t1 = it + ntetS * ir - 1;
        r1t2 = r1t1 + 1;
        r2t1 = r1t1 + ntetS;
        r2t2 = r2t1 + 1 ;

    } else {
        /*
          With t his we have B-field AND this is far more
          intuitive for me ....
        */
        r1t1 = idx(ir, it);
        r2t1 = idx(ir + 1, it);
        r1t2 = idx(ir, it + 1);
        r2t2 = idx(ir + 1, it + 1);
    }

    double bzf = 0.0, bz = 0.0 /*, bzcub = 0.0*/;
    double brf = 0.0, br = 0.0 /*, brcub = 0.0*/;
    double btf = 0.0, bt = 0.0 /*, btcub = 0.0*/;

    if((it >= 0) && (ir >= 0) && (it < Bfield.ntetS) && (ir < Bfield.nrad)) {

        /* Bz */
        bzf = (Bfield.bfld[r1t1] * wr2 * wt2 + Bfield.bfld[r2t1] * wr1 * wt2 +
               Bfield.bfld[r1t2] * wr2 * wt1 + Bfield.bfld[r2t2] * wr1 * wt1);

        // bzcub = (Bfield.f2[r1t1] * wr2 * wt2 +
        //          Bfield.f2[r2t1] * wr1 * wt2 +
        //          Bfield.f2[r1t2] * wr2 * wt1 +
        //          Bfield.f2[r2t2] * wr1 * wt1) * pow(R[2], 2.0);

        // bz = -( bzf - bzcub );
        bz = - bzf ;


        /* Br */
        brf = (Bfield.dbr[r1t1] * wr2 * wt2 +
               Bfield.dbr[r2t1] * wr1 * wt2 +
               Bfield.dbr[r1t2] * wr2 * wt1 +
               Bfield.dbr[r2t2] * wr1 * wt1) * R[2];


        // brcub = (Bfield.f3[r1t1] * wr2 * wt2 +
        //          Bfield.f3[r2t1] * wr1 * wt2 +
        //          Bfield.f3[r1t2] * wr2 * wt1 +
        //          Bfield.f3[r2t2] * wr1 * wt1) * pow(R[2], 3.0);

        // br = -( brf - brcub );
        br = - brf;


        /* Btheta */
        btf = (Bfield.dbt[r1t1] * wr2 * wt2 +
               Bfield.dbt[r2t1] * wr1 * wt2 +
               Bfield.dbt[r1t2] * wr2 * wt1 +
               Bfield.dbt[r2t2] * wr1 * wt1) / rad * R[2];


        // btcub = (Bfield.g3[r1t1] * wr2 * wt2 +
        //          Bfield.g3[r2t1] * wr1 * wt2 +
        //          Bfield.g3[r1t2] * wr2 * wt1 +
        //          Bfield.g3[r2t2] * wr1 * wt1) / rad * pow(R[2], 3.0);

        // bt = -( btf - btcub );
        bt = - btf;

        /* Br Btheta -> Bx By */
        if(tcr1_m != 0.0 && tcr2_m != 0.0) {
            double partR = sqrt(R[0] * R[0] + R[1] * R[1]);
            double Amax1 = 1;
            double Amax2 = 3;
            double Amin = -2;
            double x01 = 4;
            double x02 = 8;
            double h1 = 0.03;
            double h2 = 0.2;
            double ftc = slptc_m;
            double part1 = pow(10.0, (partR / ftc - tcr1_m / ftc - x01) * h1);
            double part2 = pow(10.0, (x02 - partR / ftc + tcr1_m / ftc) * h2);
            double part3 = -(Amax1 - Amin) * h1 * log(10) / ftc / (1 + part1) / (1 + part1) * part1;
            double part4 = (Amax2 - Amin) * h2 * log(10) / ftc / (1 + part2) / (1 + part2) * part2;
            double dr = mbtc_m / 2.78 * (part3 + part4);
            double btr = mbtc_m / 2.78 * (Amin + (Amax1 - Amin) / (1 + part1) + (Amax2 - Amin) / (1 + part2) - 1.0);

            // trim coil
            // double tca1 = 4;
            // double tca2 = 22;
            // bz=bz*(1-btr/20);
            //for (int kk=0;kk<8;++kk) {
            //if  (tet>tca1+kk&&tet<tca2+kk) {

            if(partR < tcr2_m)       {
                bz -= btr;
                br -= dr * R[2];
            }
        }

        B[0] = br * cos(tet_rad) - bt * sin(tet_rad);
        B[1] = br * sin(tet_rad) + bt * cos(tet_rad);
        B[2] = bz;

    } else {
        ERRORMSG("Error!" << getName() << ".getFieldstrength(): out of boudaries (z,r) = (" << R[0] << "," << R[1] << "," << R[2] << ")" << endl
                 << " rad=" << rad << " [mm], theta=" << tet << " [deg], it=" << it << ", ir=" << ir << endl);
        return true;

    }
    if(myBFieldType_m == BANDRF) {
      //The RF field is suppose to be sampled on a cartesian grid
        vector<Fieldmap *>::const_iterator fi  = RFfields_m.begin();
        vector<double>::const_iterator rffi    = rffrequ_m.begin();
        vector<double>::const_iterator rfphii  = rfphi_m.begin();
        vector<double>::const_iterator escali  = escale_m.begin();
        double xBegin(0), xEnd(0), yBegin(0), yEnd(0), zBegin(0), zEnd(0);
        int fcount = 0;
        for(; fi != RFfields_m.end(); ++fi, ++rffi, ++rfphii, ++escali, ++fcount) {
	  (*fi)->getFieldDimensions(xBegin, xEnd, yBegin, yEnd, zBegin, zEnd);
	  if (R(0) >= xBegin && R(0) <= xEnd && R(1) >= yBegin && R(1) <= yEnd && R(2) >= zBegin && R(2) <= zEnd) {
	    Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
	    if(!(*fi)->getFieldstrength(R, tmpE, tmpB)) {
	      double phase = 2.0 * pi * 1E-3 * (*rffi) * t + *rfphii;
	      double ebscale = *escali;
	      E += ebscale * cos(phase) * tmpE;
	      B -= ebscale * sin(phase) * tmpB;
	      //INFOMSG("Field " << fcount << " BANDRF E= " << tmpE << " R= " << R << " phase " << phase << endl);
	    }
	    if (!superpose_m) break;
	  }
    	}
    }
    return false;
}

void Cyclotron::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {
    RefPartBunch_m = bunch;
    online_m = true;
}

void Cyclotron::finalise() {
    online_m = false;
}

bool Cyclotron::bends() const {
    return true;
}





// calculate derivatives with 5-point lagrange's formula.
double Cyclotron::gutdf5d(double *f, double dx, const int kor, const int krl, const int lpr)

{
    double C[5][5][3], FAC[3];
    double result;
    int j;
    /* CALCULATE DERIVATIVES WITH 5-POINT LAGRANGE FORMULA
     * PARAMETERS:
     * F  STARTADDRESS FOR THE 5 SUPPORT POINTS
     * DX STEPWIDTH FOR ARGUMENT
     * KOR        ORDER OF DERIVATIVE (KOR=1,2,3).
     * KRL        NUMBER OF SUPPORT POINT, WHERE THE DERIVATIVE IS TO BE CALCULATED
     *  (USUALLY 3, USE FOR BOUNDARY 1 ,2, RESP. 4, 5)
     * LPR        DISTANCE OF THE 5 STORAGE POSITIONS (=1 IF THEY ARE NEIGHBORS OR LENGTH
     * OF COLUMNLENGTH OF A MATRIX, IF THE SUPPORT POINTS ARE ON A LINE).
     * ATTENTION! THE INDICES ARE NOW IN C-FORMAT AND NOT IN FORTRAN-FORMAT.*/

    /* COEFFICIENTS FOR THE 1ST DERIVATIVE: */
    C[0][0][0] = -50.0;
    C[1][0][0] = 96.0;
    C[2][0][0] = -72.0;
    C[3][0][0] = 32.0;
    C[4][0][0] = -6.0;
    C[0][1][0] = -6.0;
    C[1][1][0] = -20.0;
    C[2][1][0] = 36.0;
    C[3][1][0] = -12.0;
    C[4][1][0] =  2.0;
    C[0][2][0] =  2.0;
    C[1][2][0] = -16.0;
    C[2][2][0] =  0.0;
    C[3][2][0] = 16.0;
    C[4][2][0] = -2.0;
    C[0][3][0] = -2.0;
    C[1][3][0] = 12.0;
    C[2][3][0] = -36.0;
    C[3][3][0] = 20.0;
    C[4][3][0] =  6.0;
    C[0][4][0] =  6.0;
    C[1][4][0] = -32.0;
    C[2][4][0] = 72.0;
    C[3][4][0] = -96.0;
    C[4][4][0] = 50.0;

    /* COEFFICIENTS FOR THE 2ND DERIVATIVE: */
    C[0][0][1] = 35.0;
    C[1][0][1] = -104;
    C[2][0][1] = 114.0;
    C[3][0][1] = -56.0;
    C[4][0][1] = 11.0;
    C[0][1][1] = 11.0;
    C[1][1][1] = -20.0;
    C[2][1][1] =  6.0;
    C[3][1][1] =  4.0;
    C[4][1][1] = -1.0;
    C[0][2][1] = -1.0;
    C[1][2][1] = 16.0;
    C[2][2][1] = -30.0;
    C[3][2][1] = 16.0;
    C[4][2][1] = -1.0;
    C[0][3][1] = -1.0;
    C[1][3][1] =  4.0;
    C[2][3][1] =  6.0;
    C[3][3][1] = -20.0;
    C[4][3][1] = 11.0;
    C[0][4][1] = 11.0;
    C[1][4][1] = -56.0;
    C[2][4][1] = 114.0;
    C[3][4][1] = -104;
    C[4][4][1] = 35.0;


    /* COEFFICIENTS FOR THE 3RD DERIVATIVE: */
    C[0][0][2] = -10.0;
    C[1][0][2] = 36.0;
    C[2][0][2] = -48.0;
    C[3][0][2] = 28.0;
    C[4][0][2] = -6.0;
    C[0][1][2] = -6.0;
    C[1][1][2] = 20.0;
    C[2][1][2] = -24.0;
    C[3][1][2] = 12.0;
    C[4][1][2] = -2.0;
    C[0][2][2] = -2.0;
    C[1][2][2] =  4.0;
    C[2][2][2] =  0.0;
    C[3][2][2] = -4.0;
    C[4][2][2] =  2.0;
    C[0][3][2] =  2.0;
    C[1][3][2] = -12.0;
    C[2][3][2] = 24.0;
    C[3][3][2] = -20.0;
    C[4][3][2] =  6.0;
    C[0][4][2] =  6.0;
    C[1][4][2] = -28.0;
    C[2][4][2] = 48.0;
    C[3][4][2] = -36.0;
    C[4][4][2] = 10.0;

    /* FACTOR: */
    FAC[0] = 24.0;
    FAC[1] = 12.0;
    FAC[2] = 4.0;

    result = 0.0;
    for(j = 0; j < 5; j++) {
        result += C[j][krl][kor] * *(f + j * lpr);
    }

    return result / (FAC[kor] * pow(dx, (kor + 1)));
}


// evaulate other derivative of magnetic field.
void Cyclotron::getdiffs() {

    if(Bfield.dbr) delete[] Bfield.dbr;
    Bfield.dbr   = new double[Bfield.ntot];
    if(Bfield.dbrr) delete[] Bfield.dbrr;
    Bfield.dbrr  = new double[Bfield.ntot];
    if(Bfield.dbrrr) delete[] Bfield.dbrrr;
    Bfield.dbrrr = new double[Bfield.ntot];

    if(Bfield.dbrt) delete[] Bfield.dbrt;
    Bfield.dbrt  = new double[Bfield.ntot];
    if(Bfield.dbrrt) delete[] Bfield.dbrrt;
    Bfield.dbrrt = new double[Bfield.ntot];
    if(Bfield.dbrtt) delete[] Bfield.dbrtt;
    Bfield.dbrtt = new double[Bfield.ntot];

    if(Bfield.f2) delete[] Bfield.f2;
    Bfield.f2    = new double[Bfield.ntot];
    if(Bfield.f3) delete[] Bfield.f3;
    Bfield.f3    = new double[Bfield.ntot];
    if(Bfield.g3) delete[] Bfield.g3;
    Bfield.g3    = new double[Bfield.ntot];

    for(int i = 0; i < Bfield.nrad; i++) {

        for(int k = 0; k < Bfield.ntet; k++) {

            double dtheta = pi / 180.0 * BP.dtet;

            int kEdge;

            kEdge = max(k - 2, 0);
            kEdge = min(kEdge, Bfield.ntet - 5);

            int dkFromEdge = k - kEdge;
            int index = idx(i, k);
            int indexkEdge = idx(i, kEdge);


            Bfield.dbt[index]    = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 0, dkFromEdge, 1);
            Bfield.dbtt[index]   = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 1, dkFromEdge, 1);
            Bfield.dbttt[index]  = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 2, dkFromEdge, 1);
        }
    }



    for(int k = 0; k < Bfield.ntet; k++) {
        // inner loop varies R
        for(int i = 0; i < Bfield.nrad; i++) {
            double rac = BP.rarr[i];
            // define iredg, the reference index for radial interpolation
            // standard: i-2 minimal: 0 (not negative!)  maximal: nrad-4
            int iredg = max(i - 2, 0);
            iredg = min(iredg, Bfield.nrad - 5);
            int irtak = i - iredg;
            int index = idx(i, k);
            int indexredg = idx(iredg, k);


            Bfield.dbr[index]    = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 0, irtak, Bfield.ntetS);
            Bfield.dbrr[index]   = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 1, irtak, Bfield.ntetS);
            Bfield.dbrrr[index]  = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 2, irtak, Bfield.ntetS);

            Bfield.dbrt[index]   = gutdf5d(&Bfield.dbt[indexredg], BP.delr, 0, irtak, Bfield.ntetS);
            Bfield.dbrrt[index]  = gutdf5d(&Bfield.dbt[indexredg], BP.delr, 1, irtak, Bfield.ntetS);
            Bfield.dbrtt[index]  = gutdf5d(&Bfield.dbtt[indexredg], BP.delr, 0, irtak, Bfield.ntetS);

            // fehlt noch!! f2,f3,g3,
            Bfield.f2[index] = (Bfield.dbrr[index]
                                + Bfield.dbr[index] / rac
                                + Bfield.dbtt[index] / rac / rac) / 2.0;

            Bfield.f3[index] = (Bfield.dbrrr[index]
                                + Bfield.dbrr[index] / rac
                                + (Bfield.dbrtt[index] - Bfield.dbr[index]) / rac / rac
                                - 2.0 * Bfield.dbtt[index] / rac / rac / rac) / 6.0;

            Bfield.g3[index] = (Bfield.dbrrt[index]
                                + Bfield.dbrt[index] / rac
                                + Bfield.dbttt[index] / rac / rac) / 6.0;
        } // Radius Loop
    } // Azimuth loop

    // copy 1st azimuth to last + 1 to always yield an interval
    for(int i = 0; i < Bfield.nrad; i++) {
        int iend = idx(i, Bfield.ntet);
        int istart = idx(i, 0);

        Bfield.bfld[iend]   = Bfield.bfld[istart];
        Bfield.dbt[iend]    = Bfield.dbt[istart];
        Bfield.dbtt[iend]   = Bfield.dbtt[istart];
        Bfield.dbttt[iend]  = Bfield.dbttt[istart];

        Bfield.dbr[iend]    = Bfield.dbr[istart];
        Bfield.dbrr[iend]   = Bfield.dbrr[istart];
        Bfield.dbrrr[iend]  = Bfield.dbrrr[istart];

        Bfield.dbrt[iend]   = Bfield.dbrt[istart];
        Bfield.dbrtt[iend]  = Bfield.dbrtt[istart];
        Bfield.dbrrt[iend]  = Bfield.dbrrt[istart];

        Bfield.f2[iend]     = Bfield.f2[istart];
        Bfield.f3[iend]     = Bfield.f3[istart];
        Bfield.g3[iend]     = Bfield.g3[istart];

    }
    /*
      debug

    for(int i = 0; i< Bfield.nrad; i++){
      for(int j = 0; j< Bfield.ntetS; j++){
    int index = idx(i,j);
    double x = i*BP.delr * sin(j);
    double y = i*BP.delr * cos(j);
    *gmsg<<"x= "<<x<<" y= "<<y<<" B= "<<Bfield.bfld[index]<<endl;
      }
    }
    */
}

// read field map from external file.
void Cyclotron::getFieldFromFile(const double &scaleFactor) {

    FILE *f = NULL;
    int lpar;
    char fout[100];
    double dtmp;

    *gmsg << "* ----------------------------------------------" << endl;
    *gmsg << "*             READ IN RING FIELD MAP            " << endl;
    *gmsg << "*      (The first data block is useless)       " << endl;
    *gmsg << "* ----------------------------------------------" << endl;

    BP.Bfact = scaleFactor;

    if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
        ERRORMSG("Error in Cyclotron::getFieldFromFile()!" << endl);
        ERRORMSG(" Cannot open file, please check if it really exists." << endl);
        exit(1);
    }

    fscanf(f, "%lf", &BP.rmin);
    *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;

    fscanf(f, "%lf", &BP.delr);
    *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;

    fscanf(f, "%lf", &BP.tetmin);
    *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;

    fscanf(f, "%lf", &BP.dtet);
    //if the value is nagtive, the actual value is its reciprocal.
    if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
    *gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg.]" << endl;

    for(int i = 0; i < 13; i++)fscanf(f, "%s", fout);

    fscanf(f, "%d", &Bfield.nrad);
    *gmsg << "* Index in radial direction: " << Bfield.nrad << endl;

    fscanf(f, "%d", &Bfield.ntet);
    *gmsg << "* Index in azimuthal direction: " << Bfield.ntet << endl;

    Bfield.ntetS = Bfield.ntet + 1;
    *gmsg << "* Accordingly, total grid point along azimuth:  " << Bfield.ntetS << endl;

    for(int i = 0; i < 5; i++) {
        fscanf(f, "%s", fout);
    }
    fscanf(f, "%d", &lpar);
    // msg<< "READ"<<lpar<<" DATA ENTRIES"<<endl;

    for(int i = 0; i < 4; i++) {
        fscanf(f, "%s", fout);
    }

    for(int i = 0; i < lpar; i++) {
        fscanf(f, "%16lE", &dtmp);
    }
    for(int i = 0; i < 6; i++) {
        fscanf(f, "%s", fout);
    }
    //*gmsg << "* READ FILE DESCRIPTION..." <<endl;
    for(int i = 0; i < 10000; i++) {
        fscanf(f, "%s", fout);
        if(strcmp(fout, "LREC=") == 0)break;
    }

    for(int i = 0; i < 5; i++) {
        fscanf(f, "%s", fout);
    }
    Bfield.ntot = idx(Bfield.nrad - 1, Bfield.ntet) + 1;
    //jjyang
    *gmsg << "* Total stored grid point number ( ntetS * nrad ) : " << Bfield.ntot << endl;
    if(Bfield.bfld) delete[] Bfield.bfld;
    Bfield.bfld  = new double[Bfield.ntot];
    if(Bfield.dbt) delete[] Bfield.dbt;
    Bfield.dbt   = new double[Bfield.ntot];
    if(Bfield.dbtt) delete[] Bfield.dbtt;
    Bfield.dbtt  = new double[Bfield.ntot];
    if(Bfield.dbttt) delete[] Bfield.dbttt;
    Bfield.dbttt = new double[Bfield.ntot];

    *gmsg << "* read -in loop one block per radius" << endl;
    *gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;
    for(int i = 0; i < Bfield.nrad; i++) {

        if(i > 0) {
            for(int dummy = 0; dummy < 6; dummy++) {
                fscanf(f, "%s", fout); // INFO-LINE
            }
        }
        for(int k = 0; k < Bfield.ntet; k++) {
            fscanf(f, "%16lE", &(Bfield.bfld[idx(i, k)]));
            Bfield.bfld[idx(i, k)] *= BP.Bfact;
        }
        for(int k = 0; k < Bfield.ntet; k++) {
            fscanf(f, "%16lE", &(Bfield.dbt[idx(i, k)]));
            Bfield.dbt[idx(i, k)] *= BP.Bfact;
        }
        for(int k = 0; k < Bfield.ntet; k++) {
            fscanf(f, "%16lE", &(Bfield.dbtt[idx(i, k)]));
            Bfield.dbtt[idx(i, k)] *= BP.Bfact;
        }
        for(int k = 0; k < Bfield.ntet; k++) {
            fscanf(f, "%16lE", &(Bfield.dbttt[idx(i, k)]));
            Bfield.dbttt[idx(i, k)] *= BP.Bfact;
        }
    }
    fclose(f);


    *gmsg << "* Field Map read successfully!" << endl << endl;
}

// Calculates Radiae of initial grid.
// dimensions in [mm]!
void Cyclotron::initR(double rmin, double dr, int nrad) {
    BP.rarr = new double[nrad];
    for(int i = 0; i < nrad; i++) {
        BP.rarr[i] = rmin + i * dr;
    }
    BP.delr = dr;
}


void Cyclotron::initialise(PartBunch *bunch, const int &fieldflag, const double &scaleFactor) {
    RefPartBunch_m = bunch;
    //    PSIBF, AVFEQBF, ANSYSBF, FFAGBF
    // for your own format field, you should add your own getFieldFromFile() function by yourself.

    if(fieldflag == 1) {
        //*gmsg<<"Read field data from PSI format field map file."<<endl;
        myBFieldType_m = PSIBF;
        getFieldFromFile(scaleFactor);

    } else if(fieldflag == 2) {
        // *gmsg<<"Read data from 450MeV Carbon cyclotron field file"<<endl
        myBFieldType_m = CARBONBF;
        getFieldFromFile_Carbon(scaleFactor);

    } else if(fieldflag == 3) {
        // *gmsg<<"Read data from 100MeV H- cyclotron CYCIAE-100 field file"<<endl;
        myBFieldType_m = ANSYSBF;
        getFieldFromFile_CYCIAE(scaleFactor);

    } else if(fieldflag == 4) {
        *gmsg << "* Read AVFEQ data (Riken) use bfield scale factor bs= " << getBScale() << endl;
        myBFieldType_m = AVFEQBF;
        getFieldFromFile_AVFEQ(scaleFactor);

    } else if(fieldflag == 5) {
        *gmsg << "Read FFAG data MSU/FNAL " << getBScale() << endl;
        myBFieldType_m = FFAGBF;
        getFieldFromFile_FFAG(scaleFactor);

    } else if(fieldflag == 6) {
        *gmsg << "* Read both median plane B field map and 3D E field map of RF cavity for compact cyclotron" << getBScale() << endl;
        myBFieldType_m = BANDRF;
        getFieldFromFile_BandRF(scaleFactor);

    } else
        ERRORMSG("The field reading function of this TYPE of CYCLOTRON has not implemented yet.!" << endl);

    // calculate the radii of initial grid.
    initR(BP.rmin, BP.delr, Bfield.nrad);

    // calculate the remaining derivatives
    getdiffs();

}


void Cyclotron::getFieldFromFile_FFAG(const double &scaleFactor) {

    /*
      Field is read in from ascci file (COSY output) in the oder:
      R(m) theta(Deg) x(m) y(m) Bz(T).

      Theta is the fast varing variable

      2.0000   0.0  2.0000  0.0000      0.0000000000000000
      2.0000   1.0  1.9997  0.0349      0.0000000000000000
      2.0000   2.0  1.9988  0.0698      0.0000000000000000
      2.0000   3.0  1.9973  0.1047      0.0000000000000000

      ......
      <blank line>

      2.1000   0.0  2.1000  0.0000      0.0000000000000000
      2.1000   1.0  2.0997  0.0367      0.0000000000000000
      2.1000   2.0  2.0987  0.0733      0.0000000000000000
    */



    vector<double> rv;
    vector<double> thv;
    vector<double> xv;
    vector<double> yv;
    vector<double> bzv;
    vector<double>::iterator vit;

    *gmsg << "* ----------------------------------------------" << endl;
    *gmsg << "*    READ IN FFAG FIELD MAP     " << endl;
    *gmsg << "* ----------------------------------------------" << endl;

    BP.Bfact = -10.0; // T->kG and H- for the current FNAL FFAG

    ifstream file_to_read(fmapfn_m.c_str());
    const int max_num_of_char_in_a_line = 128;
    const int num_of_header_lines = 1;

    // STEP2: SKIP ALL THE HEADER LINES
    for(int i = 0; i < num_of_header_lines; ++i)
        file_to_read.ignore(max_num_of_char_in_a_line, '\n');

    while(!file_to_read.eof()) {
        double r, th, x, y, bz;
        file_to_read >> r >> th >> x >> y >> bz;
        if((int)th != 360) {
            rv.push_back(r * 1000.0);
            thv.push_back(th);
            xv.push_back(x * 1000.0);
            yv.push_back(y * 1000.0);
            bzv.push_back(bz);
        }
    }

    double maxtheta = 360.0;
    BP.dtet = thv[1] - thv[0];
    BP.rmin = *(rv.begin());
    double rmax = rv.back();

    // find out dR
    for(vit = rv.begin(); *vit <= BP.rmin; vit++) {}
    BP.delr = *vit - BP.rmin;

    BP.tetmin = thv[0];

    Bfield.ntet = (int)((maxtheta - thv[0]) / BP.dtet);
    Bfield.nrad  = (int)(rmax - BP.rmin) / BP.delr + 1;
    Bfield.ntetS  = Bfield.ntet + 1;
    *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;
    *gmsg << "* Maximal radius of measured field map: " << rmax << " [mm]" << endl;
    *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
    *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;
    *gmsg << "* Maximal angle of measured field map: " << maxtheta << " [deg.]" << endl;

    //if the value is negtive, the actual value is its reciprocal.
    if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
    *gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg.]" << endl;
    *gmsg << "* Total grid point along azimuth:  " << Bfield.ntetS << endl;
    *gmsg << "* Total grid point along radius: " << Bfield.nrad << endl;

    Bfield.ntot = Bfield.ntetS * Bfield.nrad;
    *gmsg << "* Total stored grid point number ( ntetS * nrad ) : " << Bfield.ntot << endl;

    if(Bfield.bfld) delete[] Bfield.bfld;
    Bfield.bfld  = new double[Bfield.ntot];
    if(Bfield.dbt) delete[] Bfield.dbt;
    Bfield.dbt   = new double[Bfield.ntot];
    if(Bfield.dbtt) delete[] Bfield.dbtt;
    Bfield.dbtt  = new double[Bfield.ntot];
    if(Bfield.dbttt) delete[] Bfield.dbttt;
    Bfield.dbttt = new double[Bfield.ntot];

    *gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;

    int count = 0;
    if(Ippl::getNodes() == 1) {
        fstream fp;
        fp.open("gnu.out", ios::out);

        for(int r = 0; r < Bfield.nrad; r++) {
            for(int k = 0; k < Bfield.ntet; k++) {
                Bfield.bfld[idx(r, k)] = bzv[count] * BP.Bfact;
                fp << BP.rmin + (r * BP.delr) << " \t " << k*(BP.tetmin + BP.dtet) << " \t " << Bfield.bfld[idx(r, k)] << endl;
                count++;
            }
        }
        fp.close();
    }

    *gmsg << "* Field Map read successfully nelem= " << count << endl << endl;
}

void Cyclotron::getFieldFromFile_AVFEQ(const double &scaleFactor) {

    FILE *f = NULL;
    *gmsg << "* ----------------------------------------------" << endl;
    *gmsg << "*    READ IN AVFEQ CYCLOTRON FIELD MAP     " << endl;
    *gmsg << "* ----------------------------------------------" << endl;

    /*  From Hiroki-san
        The first line tells r minimum (500mm),
                             r maximum(4150mm),
                             r step(50mm),
                             theta minimum(0deg),
                             theta maximum(90deg)
                             theta step(0.5deg).

        From the next line data repeat the block for a given r which the first line of the block tells.
        Each block consists of the data Bz from theta minimum (0deg) to theta maximum(90deg) with theta step(0.5deg).
    */

    BP.Bfact = scaleFactor / 1000.;

    if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
        ERRORMSG("Error in Cyclotron::getFieldFromFile_AVFEQ()!" << endl);
        ERRORMSG(" Cannot open file, please check if it really exists." << endl);
        exit(1);
    }

    fscanf(f, "%lf", &BP.rmin);
    *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;

    double rmax;
    fscanf(f, "%lf", &rmax);
    *gmsg << "* Maximal radius of measured field map: " << rmax << " [mm]" << endl;

    fscanf(f, "%lf", &BP.delr);
    *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;

    fscanf(f, "%lf", &BP.tetmin);
    *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;

    double tetmax;
    fscanf(f, "%lf", &tetmax);
    *gmsg << "* Maximal angle of measured field map: " << tetmax << " [deg.]" << endl;

    fscanf(f, "%lf", &BP.dtet);
    //if the value is nagtive, the actual value is its reciprocal.

    if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
    *gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg.]" << endl;

    Bfield.ntetS = (int)((tetmax - BP.tetmin) / BP.dtet + 1);
    *gmsg << "* Total grid point along azimuth:  " << Bfield.ntetS << endl;

    Bfield.nrad = (int)(rmax - BP.rmin) / BP.delr;


    int ntotidx = idx(Bfield.nrad, Bfield.ntetS) + 1;

    Bfield.ntot = Bfield.ntetS * Bfield.nrad;
    *gmsg << "* Total stored grid point number ( ntetS * nrad ) : " << Bfield.ntot << " ntot-idx= " << ntotidx << endl;

    if(Bfield.bfld) delete[] Bfield.bfld;
    Bfield.bfld  = new double[Bfield.ntot];
    if(Bfield.dbt) delete[] Bfield.dbt;
    Bfield.dbt   = new double[Bfield.ntot];
    if(Bfield.dbtt) delete[] Bfield.dbtt;
    Bfield.dbtt  = new double[Bfield.ntot];
    if(Bfield.dbttt) delete[] Bfield.dbttt;
    Bfield.dbttt = new double[Bfield.ntot];

    *gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;

    fstream fp;
    if(Ippl::getNodes() == 1)
        fp.open("gnu.out", ios::out);

    double tmp;
    int count = 0;

    for(int r = 0; r < Bfield.nrad; r++) {
        fscanf(f, "%16lE", &tmp);   // over read
        for(int k = 0; k < Bfield.ntetS; k++) {
            fscanf(f, "%16lE", &(Bfield.bfld[idx(r, k)]));
            Bfield.bfld[idx(r, k)] *= BP.Bfact;
            if(Ippl::getNodes() == 1)
                fp << BP.rmin + (r * BP.delr) << " \t " << k*(BP.tetmin + BP.dtet) << " \t " << Bfield.bfld[idx(r, k)] << " idx= " << idx(r, k)  << endl;
            count++;
        }
    }
    if(Ippl::getNodes() == 1)
        fp.close();
    fclose(f);
    *gmsg << "* Field Map read successfully nelem= " << count << endl << endl;
}


// read field map from external file.
void Cyclotron::getFieldFromFile_Carbon(const double &scaleFactor) {

    FILE *f = NULL;
    *gmsg << "* ----------------------------------------------" << endl;
    *gmsg << "*      READ IN CARBON CYCLOTRON FIELD MAP       " << endl;
    *gmsg << "* ----------------------------------------------" << endl;

    BP.Bfact = scaleFactor;

    if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
        ERRORMSG("* Error in Cyclotron::getFieldFromFile_Carbon()!" << endl);
        ERRORMSG(" Cannot open file, please check if it really exists." << endl);
        exit(1);
    }

    fscanf(f, "%lf", &BP.rmin);
    *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;

    fscanf(f, "%lf", &BP.delr);
    *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;

    fscanf(f, "%lf", &BP.tetmin);
    *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;

    fscanf(f, "%lf", &BP.dtet);
    //if the value is nagtive, the actual value is its reciprocal.
    if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
    *gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg.]" << endl;

    fscanf(f, "%d", &Bfield.ntet);
    *gmsg << "* Index in azimuthal direction: " << Bfield.ntet << endl;

    fscanf(f, "%d", &Bfield.nrad);
    *gmsg << "* Index in radial direction: " << Bfield.nrad << endl;

    Bfield.ntetS = Bfield.ntet + 1;
    *gmsg << "* Accordingly, total grid point along azimuth:  " << Bfield.ntetS << endl;

    Bfield.ntot = idx(Bfield.nrad - 1, Bfield.ntet) + 1;

    *gmsg << "* Total stored grid point number ( ntetS * nrad ) : " << Bfield.ntot << endl;
    if(Bfield.bfld) delete[] Bfield.bfld;
    Bfield.bfld  = new double[Bfield.ntot];
    if(Bfield.dbt) delete[] Bfield.dbt;
    Bfield.dbt   = new double[Bfield.ntot];
    if(Bfield.dbtt) delete[] Bfield.dbtt;
    Bfield.dbtt  = new double[Bfield.ntot];
    if(Bfield.dbttt) delete[] Bfield.dbttt;
    Bfield.dbttt = new double[Bfield.ntot];

    *gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;

    for(int i = 0; i < Bfield.nrad; i++) {
        for(int k = 0; k < Bfield.ntet; k++) {
            fscanf(f, "%16lE", &(Bfield.bfld[idx(i, k)]));
            Bfield.bfld[idx(i, k)] *= BP.Bfact;
        }
    }

    if(Ippl::getNodes() == 1) {
        fstream fp1, fp2;
        fp1.open("gnu.out", ios::out);
        fp2.open("eb.out", ios::out);
        for(int i = 0; i < Bfield.nrad; i++) {
            for(int k = 0; k < Bfield.ntet; k++) {
                fp1 << BP.rmin + (i * BP.delr) << " \t " << k*(BP.tetmin + BP.dtet) << " \t " << Bfield.bfld[idx(i, k)] << endl;

                Vector_t tmpR = Vector_t (BP.rmin + (i * BP.delr), 0.0, k * (BP.tetmin + BP.dtet));
                Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
                tmpR /= 1000.0; // -> mm to m
                vector<Fieldmap *>::const_iterator fi  = RFfields_m.begin();
                vector<double>::const_iterator rffi    = rffrequ_m.begin();
                vector<double>::const_iterator rfphii  = rfphi_m.begin();
                vector<double>::const_iterator escali  = escale_m.begin();
                for(; fi != RFfields_m.end(); ++fi, ++rffi, ++rfphii, ++escali) {
                    Vector_t E(0.0, 0.0, 0.0), B(0.0, 0.0, 0.0);
                    if(!(*fi)->getFieldstrength(tmpR, tmpE, tmpB)) {
                        tmpE += E;
                        tmpB -= B;
                    }
                }
                fp2 << tmpR  <<  " \t E= " << tmpE << "\t B= " << tmpB << endl;
            }
        }
        fp1.close();
        fp2.close();
    }

    fclose(f);

    *gmsg << "* Field Maps read successfully!" << endl << endl;
}


// read field map from external file.
void Cyclotron::getFieldFromFile_CYCIAE(const double &scaleFactor) {

    FILE *f = NULL;
    char fout[100];
    int dtmp;

    *gmsg << "* ----------------------------------------------" << endl;
    *gmsg << "*    READ IN CYCIAE-100 CYCLOTRON FIELD MAP     " << endl;
    *gmsg << "* ----------------------------------------------" << endl;

    BP.Bfact = scaleFactor;

    if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
        ERRORMSG("* Error in Cyclotron::getFieldFromFile_Carbon()!" << endl);
        ERRORMSG(" Cannot open file, please check if it really exists." << endl);
        exit(1);
    }

    fscanf(f, "%lf", &BP.rmin);
    *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;

    fscanf(f, "%lf", &BP.delr);
    *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;

    fscanf(f, "%lf", &BP.tetmin);
    *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;

    fscanf(f, "%lf", &BP.dtet);
    //if the value is nagtive, the actual value is its reciprocal.
    if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
    *gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg.]" << endl;

    fscanf(f, "%d", &Bfield.ntet);
    *gmsg << "* Index in azimuthal direction: " << Bfield.ntet << endl;

    fscanf(f, "%d", &Bfield.nrad);
    *gmsg << "* Index in radial direction: " << Bfield.nrad << endl;

    Bfield.ntetS = Bfield.ntet + 1;
    *gmsg << "* Accordingly, total grid point along azimuth:  " << Bfield.ntetS << endl;

    Bfield.ntot = idx(Bfield.nrad - 1, Bfield.ntet) + 1;

    *gmsg << "* Total stored grid point number ( ntetS * nrad ) : " << Bfield.ntot << endl;
    if(Bfield.bfld) delete[] Bfield.bfld;
    Bfield.bfld  = new double[Bfield.ntot];
    if(Bfield.dbt) delete[] Bfield.dbt;
    Bfield.dbt   = new double[Bfield.ntot];
    if(Bfield.dbtt) delete[] Bfield.dbtt;
    Bfield.dbtt  = new double[Bfield.ntot];
    if(Bfield.dbttt) delete[] Bfield.dbttt;
    Bfield.dbttt = new double[Bfield.ntot];

    *gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;

    int nHalfPoints = Bfield.ntet / 2.0 + 1;

    for(int i = 0; i < Bfield.nrad; i++) {

        for(int ii = 0; ii < 13; ii++)fscanf(f, "%s", fout);
        for(int k = 0; k < nHalfPoints; k++) {
            fscanf(f, "%d", &dtmp);
            fscanf(f, "%d", &dtmp);
            fscanf(f, "%d", &dtmp);
            fscanf(f, "%lf", &(Bfield.bfld[idx(i, k)]));
            Bfield.bfld[idx(i, k)] = Bfield.bfld[idx(i, k)] * (-10.0); //  T --> kGs, minus for minus hydrongen
        }
        for(int k = nHalfPoints; k < Bfield.ntet; k++) {
            Bfield.bfld[idx(i, k)] = Bfield.bfld[idx(i, Bfield.ntet-k)];
        }
    }
    //  for(int i=0; i < 300; i++) msg <<"i="<<i<<", Bfield = "<< Bfield.bfld[i]<<endl;

    fclose(f);

    *gmsg << "* Field Map read successfully!" << endl << endl;
}

void Cyclotron::getFieldFromFile_BandRF(const double &scaleFactor) {

    // read 3D E&B field data file
    Inform msg("RF Field: ");

    vector<string>::const_iterator fm    = RFfilename_m.begin();
    // loop over all field maps and superpose fields
    vector<double>::const_iterator rffi    = rffrequ_m.begin();
    vector<double>::const_iterator rfphii  = rfphi_m.begin();
    vector<double>::const_iterator escali  = escale_m.begin();
    int fcount = 0;
    
    for(; fm != RFfilename_m.end(); ++fm, ++rffi, ++rfphii, ++escali, ++fcount) {
        msg << "Load field map " << fcount << " - " << *fm << endl;
        Fieldmap *f = Fieldmap::getFieldmap(*fm, false);
        if(f == NULL) {
            ERRORMSG("* Error in Cyclotron::getFieldFromFile_BandRF()!" << endl);
            ERRORMSG(" Cannot open file, please check if it really exists." << endl);
            exit(1);
        }
        f->readMap();
        f->getInfo(&msg);
        RFfields_m.push_back(f);
    }
    // read CARBON type B field
    getFieldFromFile_Carbon(scaleFactor);
}

void Cyclotron::getDimensions(double &zBegin, double &zEnd) const
{ }