00001 #include "BFit.h"
00002 #include <cstring>
00003
00004 using namespace std;
00005 using namespace magfieldparam;
00006
00007
00008
00009
00010 #ifdef BFit_PW
00011 const double BFit::Z_nom[4] = { -2.37615687260664e-2,
00012 -1.86400109250045e-2,
00013 -1.80502358070104e-2,
00014 -1.60470955291956e-2 };
00015
00016 const double BFit::B_nom[4] = { 2.02156567013928,
00017 3.51622117206486,
00018 3.81143026675623,
00019 4.01242188708911 };
00020
00021 const double BFit::C_nom[4][16] = {{ 1.0, -3.61278802720839e-3,
00022 6.36561393690475e-6, 8.32541914664693e-5,
00023 -2.42108313492765e-6, -1.87295909297299e-5,
00024 3.06832709074461e-7, 1.91827319271226e-6,
00025 -2.15392717311725e-8, -1.25266203359502e-7,
00026 3.87507522135914e-10, 4.85518568040635e-9,
00027 4.42080729840719e-11, -8.83065447433858e-11,
00028 -2.41380148377896e-12, 0.0 },
00029 { 1.0, -5.04020236643808e-3,
00030 2.03224205921125e-6, 6.79444854179620e-5,
00031 -1.98082200052911e-6, -1.93324798138490e-5,
00032 3.15120940544812e-7, 1.82623212354924e-6,
00033 -3.30483297560429e-8, -1.13251951654739e-7,
00034 1.96974144659278e-9, 4.25153392971594e-9,
00035 -6.12986034064675e-11, -7.59031334826116e-11,
00036 6.40295019219590e-13, 0.0 },
00037 { 1.0, -5.23012318846739e-3,
00038 8.80302231241395e-7, 6.51341641212249e-5,
00039 -1.68564063895995e-6, -1.93693613146655e-5,
00040 2.58178734098114e-7, 1.81311192824207e-6,
00041 -2.79301520182866e-8, -1.11679980224632e-7,
00042 1.72615649164433e-9, 4.17328869038146e-9,
00043 -5.72514160410955e-11, -7.41998111228714e-11,
00044 7.30938527053447e-13, 0.0 },
00045 { 1.0, -5.34172971309074e-3,
00046 2.48943649506081e-7, 6.23054033447814e-5,
00047 -1.60390978074464e-6, -1.92618217244767e-5,
00048 2.42461261622770e-7, 1.78772142159379e-6,
00049 -2.61432416866515e-8, -1.09159464672341e-7,
00050 1.62705377496138e-9, 4.02967933726133e-9,
00051 -5.48168162195020e-11, -7.00249566028285e-11,
00052 8.22254619144001e-13, 0.0 }};
00053 #else
00054 const double BFit::dZ_0 = -2.62328760352034e-2;
00055 const double BFit::dZ_2 = 5.94363870284212e-4;
00056
00057 const double BFit::C_0[16] = { 1.0, -2.52864632909442e-3,
00058 8.76365790071351e-6, 9.19077286315044e-5,
00059 -2.49284256023752e-6, -1.80143891826520e-5,
00060 2.29295162454016e-7, 1.96139195659245e-6,
00061 -3.47342625923464e-9, -1.32147627969588e-7,
00062 -1.50735830442900e-9, 5.17724172101696e-9,
00063 1.54539960459831e-10, -9.30914368388717e-11,
00064 -5.20466591966397e-12, 0.0 };
00065
00066 const double BFit::C_2[16] = { 0.0, -2.96314154618866e-4,
00067 -6.04246295125223e-7, -2.22393436573694e-6,
00068 2.84133631738674e-9, -2.07090716476209e-7,
00069 2.55850963123821e-8, -1.06689136150163e-8,
00070 -5.48842256680751e-9, 1.78987539969165e-9,
00071 5.57809366992069e-10, -8.25055601520632e-11,
00072 -3.18509299957904e-11, 1.11714602344300e-12,
00073 7.90102331886296e-13, 0.0 };
00074
00075 const double BFit::C_4[16] = { 0.0, 7.57194953855834e-6,
00076 4.48169046115052e-9, 2.49606093449927e-8,
00077 3.42264285146368e-9, 7.95338846845187e-9,
00078 -1.57711106312732e-9, 1.02715424120585e-11,
00079 2.57261485255293e-10, -2.41682937761163e-11,
00080 -2.27894837943020e-11, 7.98570801347331e-13,
00081 1.17889573705870e-12, 1.64571374852252e-14,
00082 -2.60212133934707e-14, 0.0 };
00083 #endif
00084
00085 BFit::BFit()
00086 {
00087 dZ = 0.;
00088 memset(C, 0, 16*sizeof(double));
00089 rz_poly *P_base = new rz_poly(16);
00090
00091 Bz_base = new rz_poly(P_base->Diff(1));
00092 Bz_base->SetOFF(1);
00093
00094 Br_base = new rz_poly(P_base->Diff(0));
00095 Br_base->SetOFF(0);
00096
00097 delete P_base;
00098 }
00099
00100
00101 void BFit::SetField(double B)
00102 {
00103
00104
00105 unsigned int jj;
00106
00107 #ifdef BFit_PW
00108 unsigned int kk = 1;
00109 double w_0, w_1;
00110 if (B <= B_nom[0]) {
00111 dZ = Z_nom[0];
00112 for (jj = 0; jj < 16; ++jj) {
00113 C[jj] = B*C_nom[0][jj];
00114 }
00115 } else if (B >= B_nom[3]) {
00116 dZ = Z_nom[3];
00117 for (jj = 0; jj < 16; ++jj) {
00118 C[jj] = B*C_nom[3][jj];
00119 }
00120 } else {
00121 while (B_nom[kk] < B) ++kk;
00122 w_1 = (B - B_nom[kk-1])/(B_nom[kk] - B_nom[kk-1]);
00123 w_0 = 1.0 - w_1;
00124 dZ = Z_nom[kk-1]*w_0 + Z_nom[kk]*w_1;
00125 for (jj = 0; jj < 16; ++jj) {
00126 C[jj] = B*(C_nom[kk-1][jj]*w_0 + C_nom[kk][jj]*w_1);
00127 }
00128 }
00129 #else
00130 double B2 = B*B;
00131 dZ = dZ_0 + dZ_2*B2;
00132 for (jj = 0; jj < 16; ++jj) {
00133 C[jj] = B*((C_4[jj]*B2 + C_2[jj])*B2 + C_0[jj]);
00134 }
00135 #endif
00136 }
00137
00138
00139 void BFit::GetField(double r, double z, double phi,
00140 double &Br, double &Bz, double &Bphi)
00141 {
00142
00143
00144
00145
00146 double zc = z + dZ;
00147
00148 Bz = Bz_base->GetSVal(r, zc, C);
00149 Br = Br_base->GetSVal(r, zc, C+1);
00150 Bphi = 0.;
00151 }