OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
FM2DMagnetoStatic.cpp
Go to the documentation of this file.
1 #include "Fields/FM2DMagnetoStatic.h"
2 #include "Fields/Fieldmap.hpp"
3 #include "Physics/Units.h"
4 #include "Utilities/GeneralClassicException.h"
5 #include "Utilities/Util.h"
6 
7 #include <fstream>
8 #include <ios>
9 #include <cmath>
10 
11 _FM2DMagnetoStatic::_FM2DMagnetoStatic(const std::string& filename):
12  _Fieldmap(filename),
13  FieldstrengthBz_m(nullptr),
14  FieldstrengthBr_m(nullptr) {
15  std::ifstream file;
16  std::string tmpString;
17  double tmpDouble;
18 
19  Type = T2DMagnetoStatic;
20 
21  // open field map, parse it and disable element on error
22  file.open(Filename_m.c_str());
23  if(file.good()) {
24  bool parsing_passed = true;
25  try {
26  parsing_passed = interpretLine<std::string, std::string>(file,
27  tmpString,
28  tmpString);
29  } catch (GeneralClassicException &e) {
30  parsing_passed = interpretLine<std::string, std::string, std::string>(file,
31  tmpString,
32  tmpString,
33  tmpString);
34 
35  tmpString = Util::toUpper(tmpString);
36  if (tmpString != "TRUE" &&
37  tmpString != "FALSE")
38  throw GeneralClassicException("_FM2DMagnetoStatic::_FM2DMagnetoStatic",
39  "The third string on the first line of 2D field "
40  "maps has to be either TRUE or FALSE");
41 
42  normalize_m = (tmpString == "TRUE");
43  }
44 
45  if(tmpString == "ZX") {
46  swap_m = true;
47  parsing_passed = parsing_passed &&
48  interpretLine<double, double, int>(file, rbegin_m, rend_m, num_gridpr_m);
49  parsing_passed = parsing_passed &&
50  interpretLine<double, double, int>(file, zbegin_m, zend_m, num_gridpz_m);
51  } else if(tmpString == "XZ") {
52  swap_m = false;
53  parsing_passed = parsing_passed &&
54  interpretLine<double, double, int>(file, zbegin_m, zend_m, num_gridpz_m);
55  parsing_passed = parsing_passed &&
56  interpretLine<double, double, int>(file, rbegin_m, rend_m, num_gridpr_m);
57  } else {
58  std::cerr << "unknown orientation of 2D magnetostatic fieldmap" << std::endl;
59  parsing_passed = false;
60  }
61 
62  for(long i = 0; (i < (num_gridpz_m + 1) * (num_gridpr_m + 1)) && parsing_passed; ++ i) {
63  parsing_passed = parsing_passed && interpretLine<double, double>(file, tmpDouble, tmpDouble);
64  }
65 
66  parsing_passed = parsing_passed &&
67  interpreteEOF(file);
68 
69  file.close();
70  lines_read_m = 0;
71 
72  if(!parsing_passed) {
73  disableFieldmapWarning();
74  zend_m = zbegin_m - 1e-3;
75  throw GeneralClassicException("_FM2DMagnetoStatic::_FM2DMagnetoStatic",
76  "An error occured when reading the fieldmap '" + Filename_m + "'");
77  } else {
78  // conversion from cm to m
79  rbegin_m *= Units::cm2m;
80  rend_m *= Units::cm2m;
81  zbegin_m *= Units::cm2m;
82  zend_m *= Units::cm2m;
83 
84  hr_m = (rend_m - rbegin_m) / num_gridpr_m;
85  hz_m = (zend_m - zbegin_m) / num_gridpz_m;
86 
87  // num spacings -> num grid points
88  num_gridpr_m++;
89  num_gridpz_m++;
90  }
91  } else {
92  noFieldmapWarning();
93  zbegin_m = 0.0;
94  zend_m = -1e-3;
95  }
96 }
97 
98 _FM2DMagnetoStatic::~_FM2DMagnetoStatic() {
99  freeMap();
100 }
101 
102 FM2DMagnetoStatic _FM2DMagnetoStatic::create(const std::string& filename)
103 {
104  return FM2DMagnetoStatic(new _FM2DMagnetoStatic(filename));
105 }
106 
107 void _FM2DMagnetoStatic::readMap() {
108  if(FieldstrengthBz_m == nullptr) {
109  // declare variables and allocate memory
110  std::ifstream in;
111  std::string tmpString;
112 
113  FieldstrengthBz_m = new double[num_gridpz_m * num_gridpr_m];
114  FieldstrengthBr_m = new double[num_gridpz_m * num_gridpr_m];
115 
116  // read in and parse field map
117  in.open(Filename_m.c_str());
118  getLine(in, tmpString);
119  getLine(in, tmpString);
120  getLine(in, tmpString);
121 
122  if(swap_m) {
123  for(int i = 0; i < num_gridpz_m; i++) {
124  for(int j = 0; j < num_gridpr_m; j++) {
125  interpretLine<double, double>(in,
126  FieldstrengthBr_m[i + j * num_gridpz_m],
127  FieldstrengthBz_m[i + j * num_gridpz_m]);
128  }
129  }
130  } else {
131  for(int j = 0; j < num_gridpr_m; j++) {
132  for(int i = 0; i < num_gridpz_m; i++) {
133  interpretLine<double, double>(in,
134  FieldstrengthBz_m[i + j * num_gridpz_m],
135  FieldstrengthBr_m[i + j * num_gridpz_m]);
136  }
137  }
138  }
139  in.close();
140 
141  if (normalize_m) {
142  double Bzmax = 0.0;
143  // find maximum field
144  for(int i = 0; i < num_gridpz_m; ++ i) {
145  if(std::abs(FieldstrengthBz_m[i]) > Bzmax) {
146  Bzmax = std::abs(FieldstrengthBz_m[i]);
147  }
148  }
149 
150  // normalize field
151  for(int i = 0; i < num_gridpr_m * num_gridpz_m; ++ i) {
152  FieldstrengthBz_m[i] /= Bzmax;
153  FieldstrengthBr_m[i] /= Bzmax;
154  }
155  }
156 
157  INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << endl);
158  }
159 }
160 
161 void _FM2DMagnetoStatic::freeMap() {
162  if(FieldstrengthBz_m != nullptr) {
163  delete[] FieldstrengthBz_m;
164  delete[] FieldstrengthBr_m;
165 
166  FieldstrengthBz_m = nullptr;
167  FieldstrengthBr_m = nullptr;
168  }
169 }
170 
171 bool _FM2DMagnetoStatic::getFieldstrength(const Vector_t &R, Vector_t &/*E*/, Vector_t &B) const {
172  // do bi-linear interpolation
173  const double RR = std::sqrt(R(0) * R(0) + R(1) * R(1));
174 
175  const int indexr = std::abs((int)std::floor(RR / hr_m));
176  const double leverr = (RR / hr_m) - indexr;
177 
178  const int indexz = std::abs((int)std::floor((R(2) - zbegin_m) / hz_m));
179  const double leverz = (R(2) - zbegin_m) / hz_m - indexz;
180 
181  if((indexz < 0) || (indexz + 2 > num_gridpz_m))
182  return false;
183 
184  if(indexr + 2 > num_gridpr_m)
185  return true;
186 
187  const int index1 = indexz + indexr * num_gridpz_m;
188  const int index2 = index1 + num_gridpz_m;
189  const double BfieldR = (1.0 - leverz) * (1.0 - leverr) * FieldstrengthBr_m[index1]
190  + leverz * (1.0 - leverr) * FieldstrengthBr_m[index1 + 1]
191  + (1.0 - leverz) * leverr * FieldstrengthBr_m[index2]
192  + leverz * leverr * FieldstrengthBr_m[index2 + 1];
193 
194  const double BfieldZ = (1.0 - leverz) * (1.0 - leverr) * FieldstrengthBz_m[index1]
195  + leverz * (1.0 - leverr) * FieldstrengthBz_m[index1 + 1]
196  + (1.0 - leverz) * leverr * FieldstrengthBz_m[index2]
197  + leverz * leverr * FieldstrengthBz_m[index2 + 1];
198 
199  if(RR != 0) {
200  B(0) += BfieldR * R(0) / RR;
201  B(1) += BfieldR * R(1) / RR;
202  }
203  B(2) += BfieldZ;
204  return false;
205 }
206 
207 bool _FM2DMagnetoStatic::getFieldDerivative(const Vector_t &R, Vector_t &/*E*/, Vector_t &B, const DiffDirection &dir) const {
208 
209  double BfieldR, BfieldZ;
210 
211  const double RR = std::sqrt(R(0) * R(0) + R(1) * R(1));
212 
213  const int indexr = (int)std::floor(RR / hr_m);
214  const double leverr = (RR / hr_m) - indexr;
215 
216  const int indexz = (int)std::floor((R(2)) / hz_m);
217  const double leverz = (R(2) / hz_m) - indexz;
218 
219  if((indexz < 0) || (indexz + 2 > num_gridpz_m))
220  return false;
221 
222  if(indexr + 2 > num_gridpr_m)
223  return true;
224 
225  const int index1 = indexz + indexr * num_gridpz_m;
226  const int index2 = index1 + num_gridpz_m;
227  if(dir == DZ) {
228  if(indexz > 0) {
229  if(indexz < num_gridpz_m - 1) {
230  BfieldR = (1.0 - leverr) * ((FieldstrengthBr_m[index1 - 1] - 2. * FieldstrengthBr_m[index1] + FieldstrengthBr_m[index1 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
231  + (FieldstrengthBr_m[index1 + 1] - FieldstrengthBr_m[index1 - 1]) / (2.*hz_m))
232  + leverr * ((FieldstrengthBr_m[index2 - 1] - 2. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index2 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
233  + (FieldstrengthBr_m[index2 + 1] - FieldstrengthBr_m[index2 - 1]) / (2.*hz_m));
234 
235  BfieldZ = (1.0 - leverr) * ((FieldstrengthBz_m[index1 - 1] - 2. * FieldstrengthBz_m[index1] + FieldstrengthBz_m[index1 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
236  + (FieldstrengthBz_m[index1 + 1] - FieldstrengthBz_m[index1 - 1]) / (2.*hz_m))
237  + leverr * ((FieldstrengthBz_m[index2 - 1] - 2. * FieldstrengthBz_m[index2] + FieldstrengthBz_m[index2 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
238  + (FieldstrengthBz_m[index2 + 1] - FieldstrengthBz_m[index2 - 1]) / (2.*hz_m));
239  } else {
240  BfieldR = (1.0 - leverr) * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBr_m[index1 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
241  - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBr_m[index1 - 2]) / (2.*hz_m))
242  + leverr * ((FieldstrengthBr_m[index2] - 2. * FieldstrengthBr_m[index2 - 1] + FieldstrengthBr_m[index2 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
243  - (3. * FieldstrengthBr_m[index2] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBr_m[index2 - 2]) / (2.*hz_m));
244 
245  BfieldZ = (1.0 - leverr) * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index1 - 1] + FieldstrengthBz_m[index1 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
246  - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBz_m[index1 - 2]) / (2.*hz_m))
247  + leverr * ((FieldstrengthBz_m[index2] - 2. * FieldstrengthBz_m[index2 - 1] + FieldstrengthBz_m[index2 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
248  - (3. * FieldstrengthBz_m[index2] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBz_m[index2 - 2]) / (2.*hz_m));
249  }
250  } else {
251  BfieldR = (1.0 - leverr) * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index1 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
252  - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index1 + 2]) / (2.*hz_m))
253  + leverr * ((FieldstrengthBr_m[index2] - 2. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBr_m[index2 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
254  - (3. * FieldstrengthBr_m[index2] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index2 + 2]) / (2.*hz_m));
255 
256  BfieldZ = (1.0 - leverr) * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index1 + 1] + FieldstrengthBz_m[index1 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
257  - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBz_m[index1 + 2]) / (2.*hz_m))
258  + leverr * ((FieldstrengthBz_m[index2] - 2. * FieldstrengthBz_m[index2 + 1] + FieldstrengthBz_m[index2 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
259  - (3. * FieldstrengthBz_m[index2] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBz_m[index2 + 2]) / (2.*hz_m));
260  }
261  } else {
262  if(indexr > 0) {
263  const int index_1 = index1 - num_gridpz_m;
264  if(indexr < num_gridpr_m - 1) {
265  BfieldR = (1.0 - leverz) * ((FieldstrengthBr_m[index_1 + 1] - 2. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index2 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
266  + (FieldstrengthBr_m[index2 + 1] - FieldstrengthBr_m[index_1 + 1]) / (2.*hr_m))
267  + leverz * ((FieldstrengthBr_m[index_1] - 2. * FieldstrengthBr_m[index1] + FieldstrengthBr_m[index2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
268  + (FieldstrengthBr_m[index2] - FieldstrengthBr_m[index_1]) / (2.*hr_m));
269 
270  BfieldZ = (1.0 - leverz) * ((FieldstrengthBz_m[index_1 + 1] - 2. * FieldstrengthBz_m[index1 + 1] + FieldstrengthBz_m[index2 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
271  + (FieldstrengthBz_m[index2 + 1] - FieldstrengthBz_m[index_1 + 1]) / (2.*hr_m))
272  + leverz * ((FieldstrengthBz_m[index_1] - 2. * FieldstrengthBz_m[index1] + FieldstrengthBz_m[index2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
273  + (FieldstrengthBz_m[index2] - FieldstrengthBz_m[index_1]) / (2.*hr_m));
274  } else {
275  const int index_2 = index_1 - num_gridpz_m;
276  BfieldR = (1.0 - leverz) * ((FieldstrengthBr_m[index1 + 1] - 2. * FieldstrengthBr_m[index_1 + 1] + FieldstrengthBr_m[index_2 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
277  - (3. * FieldstrengthBr_m[index1 + 1] - 4. * FieldstrengthBr_m[index_1 + 1] + FieldstrengthBr_m[index_2 + 1]) / (2.*hr_m))
278  + leverz * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index_1] + FieldstrengthBr_m[index_2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
279  - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index_1] + FieldstrengthBr_m[index_2]) / (2.*hr_m));
280 
281  BfieldZ = (1.0 - leverz) * ((FieldstrengthBz_m[index1 + 1] - 2. * FieldstrengthBz_m[index_1 + 1] + FieldstrengthBz_m[index_2 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
282  - (3. * FieldstrengthBz_m[index1 + 1] - 4. * FieldstrengthBr_m[index_1 + 1] + FieldstrengthBz_m[index_2 + 1]) / (2.*hr_m))
283  + leverz * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index_1] + FieldstrengthBz_m[index_2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
284  - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index_1] + FieldstrengthBz_m[index_2]) / (2.*hr_m));
285  }
286  } else {
287  const int index3 = index2 + num_gridpz_m;
288  BfieldR = (1.0 - leverz) * ((FieldstrengthBr_m[index1 + 1] - 2. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBr_m[index3 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
289  - (3. * FieldstrengthBr_m[index1 + 1] - 4. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBr_m[index3 + 1]) / (2.*hr_m))
290  + leverz * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index3]) * (RR - indexr * hr_m) / (hr_m * hr_m)
291  - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index3]) / (2.*hr_m));
292 
293  BfieldZ = (1.0 - leverz) * ((FieldstrengthBz_m[index1 + 1] - 2. * FieldstrengthBz_m[index2 + 1] + FieldstrengthBz_m[index3 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
294  - (3. * FieldstrengthBz_m[index1 + 1] - 4. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBz_m[index3 + 1]) / (2.*hr_m))
295  + leverz * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index2] + FieldstrengthBz_m[index3]) * (RR - indexr * hr_m) / (hr_m * hr_m)
296  - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index2] + FieldstrengthBz_m[index3]) / (2.*hr_m));
297  }
298 
299  }
300 
301  if(RR > 1.e-12) {
302  B(0) += BfieldR * R(0) / RR;
303  B(1) += BfieldR * R(1) / RR;
304  }
305  B(2) += BfieldZ;
306  return false;
307 }
308 
309 void _FM2DMagnetoStatic::getFieldDimensions(double &zBegin, double &zEnd) const {
310  zBegin = zbegin_m;
311  zEnd = zend_m;
312 }
313 
314 void _FM2DMagnetoStatic::getFieldDimensions(double &/*xIni*/, double &/*xFinal*/, double &/*yIni*/, double &/*yFinal*/, double &/*zIni*/, double &/*zFinal*/) const {}
315 
316 void _FM2DMagnetoStatic::swap() {
317  if(swap_m) swap_m = false;
318  else swap_m = true;
319 }
320 
321 void _FM2DMagnetoStatic::getInfo(Inform *msg) {
322  (*msg) << Filename_m << " (2D magnetostatic); zini= " << zbegin_m << " m; zfinal= " << zend_m << " m;" << endl;
323 }
324 
325 double _FM2DMagnetoStatic::getFrequency() const {
326  return 0.0;
327 }
328 
329 void _FM2DMagnetoStatic::setFrequency(double /*freq*/)
330 { ;}
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
DZ
Definition: Fieldmap.h:58
_Fieldmap::Type
MapType Type
Definition: Fieldmap.h:115
_FM2DMagnetoStatic::getFieldDimensions
virtual void getFieldDimensions(double &zBegin, double &zEnd) const
Definition: FM2DMagnetoStatic.cpp:309
_FM2DMagnetoStatic::setFrequency
virtual void setFrequency(double freq)
Definition: FM2DMagnetoStatic.cpp:329
DiffDirection
DiffDirection
Definition: Fieldmap.h:55
_FM2DMagnetoStatic::getFieldstrength
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
Definition: FM2DMagnetoStatic.cpp:171
_Fieldmap::noFieldmapWarning
void noFieldmapWarning()
Definition: Fieldmap.cpp:618
_FM2DMagnetoStatic::getFrequency
virtual double getFrequency() const
Definition: FM2DMagnetoStatic.cpp:325
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
Units::cm2m
constexpr double cm2m
Definition: Units.h:35
FM2DMagnetoStatic
std::shared_ptr< _FM2DMagnetoStatic > FM2DMagnetoStatic
Definition: Definitions.h:57
_Fieldmap::lines_read_m
int lines_read_m
Definition: Fieldmap.h:119
_Fieldmap::typeset_msg
static std::string typeset_msg(const std::string &msg, const std::string &title)
Definition: Fieldmap.cpp:649
Util::toUpper
std::string toUpper(const std::string &str)
Definition: Util.cpp:147
INFOMSG
#define INFOMSG(msg)
Definition: IpplInfo.h:348
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
_Fieldmap::interpreteEOF
bool interpreteEOF(std::ifstream &in)
Definition: Fieldmap.cpp:555
_Fieldmap::normalize_m
bool normalize_m
Definition: Fieldmap.h:121
Inform
Definition: Inform.h:42
_FM2DMagnetoStatic::create
static FM2DMagnetoStatic create(const std::string &filename)
Definition: FM2DMagnetoStatic.cpp:102
_FM2DMagnetoStatic::getInfo
virtual void getInfo(Inform *msg)
Definition: FM2DMagnetoStatic.cpp:321
_Fieldmap::disableFieldmapWarning
void disableFieldmapWarning()
Definition: Fieldmap.cpp:610
_Fieldmap::getLine
void getLine(std::ifstream &in, std::string &buffer)
Definition: Fieldmap.h:122
level3
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
Physics::e
constexpr double e
The value of .
Definition: Physics.h:39
_Fieldmap::Filename_m
std::string Filename_m
Definition: Fieldmap.h:118
_FM2DMagnetoStatic::getFieldDerivative
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
Definition: FM2DMagnetoStatic.cpp:207
_FM2DMagnetoStatic::_FM2DMagnetoStatic
_FM2DMagnetoStatic(const std::string &filename)
Definition: FM2DMagnetoStatic.cpp:11
floor
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Definition: PETE.h:733
Generated on Tue Jul 2 2024 15:05:21 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5