ParametrizedEngine/src/BFit3D.cc

0001 #include <cmath>
0002 #ifndef M_PI
0003 #define M_PI 3.1415926535897932385
0004 #endif
0005
0006 #include "BFit3D_data.h"
0007
0008 using namespace magfieldparam;
0009
0010 //_______________________________________________________________________________
0011 void BFit3D::SetCoeff_Linear(const double B) {
0012   unsigned jj;
0013   double w_0, w_1, B_mod = fabs(B);
0014   if (B_mod <= B_nom[0]) {
0015     w_0 = B / B_nom[0];
0016     for (jj = 0; jj < 360; ++jj) {
0017       C[jj] = w_0 * C0[jj][0];
0018     }
0019   } else if (B_mod >= B_nom[3]) {
0020     w_0 = B / B_nom[3];
0021     for (jj = 0; jj < 360; ++jj) {
0022       C[jj] = w_0 * C0[jj][3];
0023     }
0024   } else {
0025     unsigned kk = 1;
0026     [[clang::suppress]]
0027     while (B_nom[kk] < B_mod)
0028       ++kk;  //Simple linear search
0029     w_1 = (B_mod - B_nom[kk - 1]) / (B_nom[kk] - B_nom[kk - 1]);
0030     w_0 = 1.0 - w_1;
0031     if (B < 0.) {
0032       w_0 = -w_0;
0033       w_1 = -w_1;
0034     }
0035     for (jj = 0; jj < 360; ++jj) {
0036       C[jj] = w_0 * C0[jj][kk - 1] + w_1 * C0[jj][kk];
0037     }
0038   }
0039 }
0040
0041 //_______________________________________________________________________________
0042 void BFit3D::SetCoeff_Spline(const double B) {
0043   int jc;
0044   double dB2, dB = fabs(B);
0045   if (dB >= B_nom[3]) {  //Linear extrapolation for a large field
0046     dB -= B_nom[3];
0047     for (jc = 0; jc < 360; ++jc)
0048       C[jc] = C0[jc][3] + C1[jc][4] * dB;
0049   } else {
0050     if (dB < B_nom[0]) {
0051       dB2 = dB * dB / (3. * B_nom[0]);
0052       for (jc = 0; jc < 360; ++jc)
0053         C[jc] = (C2[jc][0] * dB2 + C1[jc][0]) * dB;
0054     } else {
0055       int k1 = 1;
0056       [[clang::suppress]]
0057       while (B_nom[k1] < dB)
0058         ++k1;  //Simple linear search
0059       int k0 = k1 - 1;
0060       dB2 = (dB -= B_nom[k0]) / (3. * (B_nom[k1] - B_nom[k0]));
0061       if (k1 < 3) {  //"Regular" interval
0062         for (jc = 0; jc < 360; ++jc)
0063           C[jc] = (((C2[jc][k1] - C2[jc][k0]) * dB2 + C2[jc][k0]) * dB + C1[jc][k1]) * dB + C0[jc][k0];
0064       } else {  //The last interval
0065         dB2 = (1. - dB2) * dB;
0066         for (jc = 0; jc < 360; ++jc)
0067           C[jc] = (C2[jc][k0] * dB2 + C1[jc][k1]) * dB + C0[jc][k0];
0068       }
0069     }
0070   }
0071   if (B < 0)
0072     for (jc = 0; jc < 360; ++jc)
0073       C[jc] = -C[jc];
0074 }
0075
0076 //_______________________________________________________________________________
0077 void BFit3D::GetField(const double r, const double z, const double phi, double &Br, double &Bz, double &Bphi) {
0078   //Return field components in Br, Bz, Bphi. SetField must be called before use.
0079   //
0080   if (signed_rad || (r >= 0.))
0081     HB->SetPoint(r, z, phi);
0082   else
0083     HB->SetPoint(-r, z, phi + M_PI);
0084   HB->EvalBr();
0085   HB->EvalBz();
0086   HB->EvalBphi();
0087   Br = HB->GetBr(C);
0088   Bz = HB->GetBz(C);
0089   Bphi = HB->GetBphi(C);
0090 }
0091
0092 /////////////////////////////////////////////////////////////////////////////////
0093 //                                                                             //
0094 //                           T E S T  A R E A                                  //
0095 //                                                                             //
0096 /////////////////////////////////////////////////////////////////////////////////