BFit.cc

Go to the documentation of this file.

#include "BFit.h"
#include <cstring>

using namespace std;
using namespace magfieldparam;

//_______________________________________________________________________________

#ifdef BFit_PW
const double BFit::Z_nom[4] = {-2.37615687260664e-2, -1.86400109250045e-2, -1.80502358070104e-2, -1.60470955291956e-2};

const double BFit::B_nom[4] = {2.02156567013928,  //Nominal foeld values
                               3.51622117206486,
                               3.81143026675623,
                               4.01242188708911};

const double BFit::C_nom[4][16] = {{1.0,
                                    -3.61278802720839e-3,
                                    6.36561393690475e-6,
                                    8.32541914664693e-5,
                                    -2.42108313492765e-6,
                                    -1.87295909297299e-5,
                                    3.06832709074461e-7,
                                    1.91827319271226e-6,
                                    -2.15392717311725e-8,
                                    -1.25266203359502e-7,
                                    3.87507522135914e-10,
                                    4.85518568040635e-9,
                                    4.42080729840719e-11,
                                    -8.83065447433858e-11,
                                    -2.41380148377896e-12,
                                    0.0},
                                   {1.0,
                                    -5.04020236643808e-3,
                                    2.03224205921125e-6,
                                    6.79444854179620e-5,
                                    -1.98082200052911e-6,
                                    -1.93324798138490e-5,
                                    3.15120940544812e-7,
                                    1.82623212354924e-6,
                                    -3.30483297560429e-8,
                                    -1.13251951654739e-7,
                                    1.96974144659278e-9,
                                    4.25153392971594e-9,
                                    -6.12986034064675e-11,
                                    -7.59031334826116e-11,
                                    6.40295019219590e-13,
                                    0.0},
                                   {1.0,
                                    -5.23012318846739e-3,
                                    8.80302231241395e-7,
                                    6.51341641212249e-5,
                                    -1.68564063895995e-6,
                                    -1.93693613146655e-5,
                                    2.58178734098114e-7,
                                    1.81311192824207e-6,
                                    -2.79301520182866e-8,
                                    -1.11679980224632e-7,
                                    1.72615649164433e-9,
                                    4.17328869038146e-9,
                                    -5.72514160410955e-11,
                                    -7.41998111228714e-11,
                                    7.30938527053447e-13,
                                    0.0},
                                   {1.0,
                                    -5.34172971309074e-3,
                                    2.48943649506081e-7,
                                    6.23054033447814e-5,
                                    -1.60390978074464e-6,
                                    -1.92618217244767e-5,
                                    2.42461261622770e-7,
                                    1.78772142159379e-6,
                                    -2.61432416866515e-8,
                                    -1.09159464672341e-7,
                                    1.62705377496138e-9,
                                    4.02967933726133e-9,
                                    -5.48168162195020e-11,
                                    -7.00249566028285e-11,
                                    8.22254619144001e-13,
                                    0.0}};
#else
const double BFit::dZ_0 = -2.62328760352034e-2;
const double BFit::dZ_2 = 5.94363870284212e-4;

const double BFit::C_0[16] = {1.0,
                              -2.52864632909442e-3,
                              8.76365790071351e-6,
                              9.19077286315044e-5,
                              -2.49284256023752e-6,
                              -1.80143891826520e-5,
                              2.29295162454016e-7,
                              1.96139195659245e-6,
                              -3.47342625923464e-9,
                              -1.32147627969588e-7,
                              -1.50735830442900e-9,
                              5.17724172101696e-9,
                              1.54539960459831e-10,
                              -9.30914368388717e-11,
                              -5.20466591966397e-12,
                              0.0};

const double BFit::C_2[16] = {0.0,
                              -2.96314154618866e-4,
                              -6.04246295125223e-7,
                              -2.22393436573694e-6,
                              2.84133631738674e-9,
                              -2.07090716476209e-7,
                              2.55850963123821e-8,
                              -1.06689136150163e-8,
                              -5.48842256680751e-9,
                              1.78987539969165e-9,
                              5.57809366992069e-10,
                              -8.25055601520632e-11,
                              -3.18509299957904e-11,
                              1.11714602344300e-12,
                              7.90102331886296e-13,
                              0.0};

const double BFit::C_4[16] = {0.0,
                              7.57194953855834e-6,
                              4.48169046115052e-9,
                              2.49606093449927e-8,
                              3.42264285146368e-9,
                              7.95338846845187e-9,
                              -1.57711106312732e-9,
                              1.02715424120585e-11,
                              2.57261485255293e-10,
                              -2.41682937761163e-11,
                              -2.27894837943020e-11,
                              7.98570801347331e-13,
                              1.17889573705870e-12,
                              1.64571374852252e-14,
                              -2.60212133934707e-14,
                              0.0};
#endif
//_______________________________________________________________________________
BFit::BFit() {
  dZ = 0.;
  memset(C, 0, 16 * sizeof(double));
  rz_poly *P_base = new rz_poly(16);  //Potential basis

  Bz_base = new rz_poly(P_base->Diff(1));  //Bz basis
  Bz_base->SetOFF(1);                      //Switch off linear term

  Br_base = new rz_poly(P_base->Diff(0));  //Br basis is shifted, so
  Br_base->SetOFF(0);                      //"0" term is ignored

  delete P_base;
}

//_______________________________________________________________________________
void BFit::SetField(double B) {
  //Set nominal field [Tesla]
  //
  unsigned int jj;

#ifdef BFit_PW
  unsigned int kk = 1;
  double w_0, w_1;
  if (B <= B_nom[0]) {
    dZ = Z_nom[0];
    for (jj = 0; jj < 16; ++jj) {
      C[jj] = B * C_nom[0][jj];
    }
  } else if (B >= B_nom[3]) {
    dZ = Z_nom[3];
    for (jj = 0; jj < 16; ++jj) {
      C[jj] = B * C_nom[3][jj];
    }
  } else {
    while (B_nom[kk] < B)
      ++kk;
    w_1 = (B - B_nom[kk - 1]) / (B_nom[kk] - B_nom[kk - 1]);
    w_0 = 1.0 - w_1;
    dZ = Z_nom[kk - 1] * w_0 + Z_nom[kk] * w_1;
    for (jj = 0; jj < 16; ++jj) {
      C[jj] = B * (C_nom[kk - 1][jj] * w_0 + C_nom[kk][jj] * w_1);
    }
  }
#else
  double B2 = B * B;
  dZ = dZ_0 + dZ_2 * B2;
  for (jj = 0; jj < 16; ++jj) {
    C[jj] = B * ((C_4[jj] * B2 + C_2[jj]) * B2 + C_0[jj]);
  }
#endif
}

//_______________________________________________________________________________
void BFit::GetField(double r, double z, double phi, double &Br, double &Bz, double &Bphi) const {
  //Get field components in the point (r,z,phi). Always return Bphi=0.
  //Parameters phi and Bphi introduced in order to keep interface
  //compatible with future trully 3D version
  //
  double zc = z + dZ;

  Bz = Bz_base->GetSVal(r, zc, C);
  Br = Br_base->GetSVal(r, zc, C + 1);
  Bphi = 0.;
}