11 #include <geometry/bfieldmap/BFieldComponentRadial.h>
12 #include <geometry/bfieldmap/BFieldComponentBeamline.h>
14 #include <framework/utilities/FileSystem.h>
15 #include <framework/logging/Logger.h>
17 #include <boost/iostreams/filtering_stream.hpp>
18 #include <boost/iostreams/device/file.hpp>
19 #include <boost/iostreams/filter/gzip.hpp>
25 namespace io = boost::iostreams;
27 void BFieldComponentRadial::initialize()
29 if (m_mapFilename.empty()) {
30 B2ERROR(
"The filename for the radial magnetic field component is empty !");
34 string fullPath = FileSystem::findFile(
"/data/" + m_mapFilename);
36 if (!FileSystem::fileExists(fullPath)) {
37 B2ERROR(
"The radial magnetic field map file '" << m_mapFilename <<
"' could not be found !");
42 io::filtering_istream fieldMapFile;
43 fieldMapFile.push(io::gzip_decompressor());
44 fieldMapFile.push(io::file_source(fullPath));
47 B2DEBUG(10,
"Loading the radial magnetic field from file '" << m_mapFilename <<
"' in to the memory...");
49 m_mapBuffer.reserve(m_mapSize[0]*m_mapSize[1]);
52 for (
int i = 0; i < m_mapSize[0]; ++i) {
53 for (
int j = 0; j < m_mapSize[1]; j++) {
54 fieldMapFile >> r >> z >> Br >> Bz;
56 m_mapBuffer.push_back({Br, Bz});
60 m_igridPitchR = 1 / m_gridPitchR;
61 m_igridPitchZ = 1 / m_gridPitchZ;
62 m_Layer = m_gapHeight + m_ironPlateThickness;
63 m_iLayer = 1 / m_Layer;
65 B2DEBUG(10,
"... loaded " << m_mapSize[0] <<
"x" << m_mapSize[1] <<
" (r,z) elements.");
76 return {v.r * a, v.z * a};
83 return {u.r + v.r, u.z + v.z};
94 if (BFieldComponentBeamline::Instance().isInRange(point))
return B;
99 if ((z < m_mapRegionZ[0]) || (z > m_mapRegionZ[1]))
return B;
102 double r2 = point.
Perp2();
104 if ((r2 < m_mapRegionR[0]*m_mapRegionR[0]) || (r2 > m_mapRegionR[1]*m_mapRegionR[1]))
return B;
108 double wr = (r - m_mapRegionR[0]) * m_igridPitchR;
109 double wz = (z - m_mapRegionZ[0]) * m_igridPitchZ;
111 int ir =
static_cast<int>(wr);
112 int iz =
static_cast<int>(wz);
114 double za = z - m_mapOffset;
115 double dz_eklm = abs(za) - (m_endyokeZMin - m_gapHeight);
116 if (dz_eklm > 0 && r > m_slotRMin) {
117 double wl = dz_eklm * m_iLayer;
118 int il =
static_cast<int>(wl);
120 double L = m_Layer * il;
121 double l = dz_eklm - L;
122 int ingap = l < m_gapHeight;
123 ir += ingap & (r < m_slotRMin + m_gridPitchR);
125 double zlg = L + m_endyokeZMin, zl0 = zlg - m_gapHeight, zl1 = zlg + m_ironPlateThickness;
131 double zoff = m_mapOffset - m_mapRegionZ[0];
132 int izl0 =
static_cast<int>((zl0 + zoff) * m_igridPitchZ);
133 int izlg =
static_cast<int>((zlg + zoff) * m_igridPitchZ);
134 int izl1 =
static_cast<int>((zl1 + zoff) * m_igridPitchZ);
136 int idz = (izl0 == iz) - (izl1 == iz) + (ingap ? -ig : ig);
137 iz += (za > 0) ? idz : -idz;
142 ir = min(m_mapSize[0] - 2, ir);
143 iz = min(m_mapSize[1] - 2, iz);
148 double wz0 = 1 - wz, wr0 = 1 - wr;
150 const unsigned int strideZ = m_mapSize[1];
151 const unsigned int i00 = ir * strideZ + iz, i01 = i00 + 1, i10 = i00 + strideZ, i11 = i01 + strideZ;
153 double w00 = wz0 * wr0, w01 = wz * wr0, w10 = wz0 * wr, w11 = wz * wr;
154 BFieldPoint Brz = H[i00] * w00 + H[i01] * w01 + H[i10] * w10 + H[i11] * w11;
157 double norm = Brz.
r / r;
158 B.SetXYZ(point.
X()*norm, point.
Y()*norm, Brz.
z);
166 void BFieldComponentRadial::terminate()