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