11 #include <geometry/dbobjects/MagneticFieldComponent3D.h>
29 template<
int ORDER = 4> constexpr
double fast_atan2_minimax(
double y,
double x)
31 static_assert(ORDER >= 1 && ORDER <= 6,
"fast_atan2_minimax: Only orders 1-6 are supported");
32 constexpr
double pi4 = M_PI / 4, pi2 = M_PI / 2;
33 const double ax = std::abs(x), ay = std::abs(y);
34 const double z = (ax - ay) / (ay + ax);
35 const double v = std::abs(z),
v2 = v * v;
39 }
else if (ORDER == 3) {
42 }
else if (ORDER == 4) {
45 constexpr
double c4[] = {
51 p0 = c4[0] +
v2 * c4[2];
52 p1 = c4[1] +
v2 * c4[3];
53 }
else if (ORDER == 5) {
56 constexpr
double c5[] = {
63 p0 = c5[0] +
v2 * (c5[2] +
v2 * c5[4]);
64 p1 = c5[1] +
v2 * (c5[3]);
65 }
else if (ORDER == 6) {
68 constexpr
double c6[] = {
76 p0 = c6[0] +
v2 * (c6[2] +
v2 * (c6[4]));
77 p1 = c6[1] +
v2 * (c6[3] +
v2 * (c6[5]));
79 double phi = (z + 1) * pi4 - (z * (v - 1)) * (p0 + v * p1);
80 phi = pi2 - copysign(phi, x);
81 return copysign(phi, y);
108 auto getIndexWeight = [
this](
double value,
int coordinate) -> std::tuple<unsigned int, double> {
110 auto index =
static_cast<unsigned int>(weight);
111 index = std::min(index,
static_cast<unsigned int>(
m_mapSize[coordinate] - 2));
113 return std::make_tuple(index, weight);
116 const double z = pos.z();
117 const double r2 = pos.Perp2();
121 const auto iz = getIndexWeight(
m_maxZ - z, 2);
124 if (r2 == 0)
return interpolate(0, 0, get<0>(iz), 0, 0, get<1>(iz));
126 const double r = std::sqrt(r2);
127 const auto ir = getIndexWeight(r -
m_minR, 0);
130 const double ay = std::abs(pos.y());
131 const auto iphi = getIndexWeight(fast_atan2_minimax<4>(ay, pos.x()), 1);
134 B2Vector3D b =
interpolate(get<0>(ir), get<0>(iphi), get<0>(iz), get<1>(ir), get<1>(iphi), get<1>(iz));
136 const double norm = 1 / r;
137 const double s = ay * norm;
138 const double c = pos.x() * norm;
140 const double sgny = (pos.y() >= 0) - (pos.y() < 0);
142 return B2Vector3D(-(b.x() * c - b.y() * s), -sgny * (b.x() * s + b.y() * c), b.z());
146 double wr1,
double wphi1,
double wz1)
const
149 const unsigned int strideR =
m_mapSize[1];
151 const double wz0 = 1 - wz1;
152 const double wr0 = 1 - wr1;
153 const double wphi0 = 1 - wphi1;
154 const unsigned int j000 = iz * strideZ + ir * strideR + iphi;
155 const unsigned int j001 = j000 + 1;
156 const unsigned int j010 = j000 + strideR;
157 const unsigned int j011 = j001 + strideR;
158 const unsigned int j100 = j000 + strideZ;
159 const unsigned int j101 = j001 + strideZ;
160 const unsigned int j110 = j010 + strideZ;
161 const unsigned int j111 = j011 + strideZ;
162 const double w00 = wphi0 * wr0;
163 const double w10 = wphi0 * wr1;
164 const double w01 = wphi1 * wr0;
165 const double w11 = wphi1 * wr1;
166 const std::vector<B2Vector3F>& B =
m_bmap;
168 (B[j000] * w00 + B[j001] * w01 + B[j010] * w10 + B[j011] * w11) * wz0 +
169 (B[j100] * w00 + B[j101] * w01 + B[j110] * w10 + B[j111] * w11) * wz1;