MagneticField.cpp 4.5 KB
Newer Older
1
#include "MagneticField.h"
frey_m's avatar
frey_m committed
2

3 4
MagneticField::MagneticField(const std::string fmapfn,
                             const double& symmetry) :
frey_m's avatar
frey_m committed
5 6 7 8 9 10 11
    Cyclotron(),
    nullGeom_m(NullGeometry()),
    nullField_m(NullField())
{
    this->setFieldMapFN(fmapfn);
    this->setSymmetry(symmetry);
}
frey_m's avatar
frey_m committed
12

13
void MagneticField::read(const std::string& type,
14
                         const double& scaleFactor) {
frey_m's avatar
frey_m committed
15
    
frey_m's avatar
frey_m committed
16 17 18 19 20 21 22 23 24 25 26 27 28 29
    if ( type == "CARBONCYCL" )
        this->getFieldFromFile_Carbon(scaleFactor);
    else if ( type == "CYCIAE" )
        this->getFieldFromFile_CYCIAE(scaleFactor);
    else if ( type == "AVFEQ" )
        this->getFieldFromFile_AVFEQ(scaleFactor);
    else if ( type == "FFAG" )
        this->getFieldFromFile_FFAG(scaleFactor);
    else if ( type == "BANDRF" )
        this->getFieldFromFile_BandRF(scaleFactor);
    else if ( type == "SYNCHROCYCLOTRON" )
        this->getFieldFromFile_Synchrocyclotron(scaleFactor);
    else
        this->getFieldFromFile(scaleFactor);
frey_m's avatar
frey_m committed
30
    
frey_m's avatar
frey_m committed
31 32 33
    if ( BP.rarr.empty() ) {
        // calculate the radii of initial grid.
        this->initR(BP.rmin, BP.delr, Bfield.nrad);
frey_m's avatar
frey_m committed
34 35
    }
    
frey_m's avatar
frey_m committed
36 37 38
    if ( Bfield.dbr.empty() || Bfield.dbt.empty() ) {
        // calculate the remaining derivatives
        this->getdiffs();
frey_m's avatar
frey_m committed
39 40 41 42
    }
}


43
void MagneticField::interpolate(double& bint,
44 45 46 47 48
                                double& brint,
                                double& btint,
                                const double& rad,
                                const double& theta
                                )
frey_m's avatar
frey_m committed
49 50 51 52
{
    // x horizontal
    // y longitudinal
    // z is vertical
53
    const double xir = (rad - BP.rmin) / (BP.delr);
frey_m's avatar
frey_m committed
54

55
    // ir : the number of path whose radius is less than the 4 points of cell which surround the particle.
frey_m's avatar
frey_m committed
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
    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;
    

    double tet_rad = theta;

    // the actual angle of particle
//     tet_rad = theta / Physics::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.
frey_m's avatar
frey_m committed
71
    double tet_map = fmod(tet_rad, 360.0 / this->getSymmetry());
frey_m's avatar
frey_m committed
72 73 74 75 76 77 78 79

    double xit = tet_map / BP.dtet;

    int it = (int) xit;

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

80
    // it : the number of point on the inner path whose angle is less than the particle' corresponding angle.
frey_m's avatar
frey_m committed
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
    // 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 this 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);
//     }

    bint = 0.0;
    brint = 0.0;
    btint = 0.0;

    if((it >= 0) && (ir >= 0) && (it < Bfield.ntetS) && (ir < Bfield.nrad)) {
        
        // dB_{z}/dr
        brint = (Bfield.dbr[r1t1] * wr2 * wt2 +
                Bfield.dbr[r2t1] * wr1 * wt2 +
                Bfield.dbr[r1t2] * wr2 * wt1 +
120
                Bfield.dbr[r2t2] * wr1 * wt1);
frey_m's avatar
frey_m committed
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
        
        // dB_{z}/dtheta
        btint = Bfield.dbt[r1t1] * wr2 * wt2 +
                Bfield.dbt[r2t1] * wr1 * wt2 +
                Bfield.dbt[r1t2] * wr2 * wt1 +
                Bfield.dbt[r2t2] * wr1 * wt1;
        
        // B_{z}
        bint = Bfield.bfld[r1t1] * wr2 * wt2 +
               Bfield.bfld[r2t1] * wr1 * wt2 +
               Bfield.bfld[r1t2] * wr2 * wt1 +
               Bfield.bfld[r2t2] * wr1 * wt1;
    }
}

136
double MagneticField::getSlices() const {
frey_m's avatar
frey_m committed
137
    return -1.0;
frey_m's avatar
frey_m committed
138 139
}

140
double MagneticField::getStepsize() const {
frey_m's avatar
frey_m committed
141
    return -1.0;
frey_m's avatar
frey_m committed
142 143
}

144
const EMField &MagneticField::getField() const {
frey_m's avatar
frey_m committed
145 146
    return nullField_m;
}
frey_m's avatar
frey_m committed
147

148
EMField &MagneticField::getField() {
frey_m's avatar
frey_m committed
149 150
    return nullField_m;
}
frey_m's avatar
frey_m committed
151

152
ElementBase* MagneticField::clone() const {
frey_m's avatar
frey_m committed
153 154
    return nullptr;
}
frey_m's avatar
frey_m committed
155

156
const BGeometryBase &MagneticField::getGeometry() const {
frey_m's avatar
frey_m committed
157 158
    return nullGeom_m;
}
frey_m's avatar
frey_m committed
159

160
BGeometryBase &MagneticField::getGeometry() {
frey_m's avatar
frey_m committed
161
    return nullGeom_m;
frey_m's avatar
frey_m committed
162
}