OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
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 
9 extern Inform *gmsg;
10 
11 FM2DMagnetoStatic::FM2DMagnetoStatic(std::string aFilename):
12  Fieldmap(aFilename),
13  FieldstrengthBz_m(NULL),
14  FieldstrengthBr_m(NULL) {
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 = interpreteLine<std::string, std::string>(file,
27  tmpString,
28  tmpString);
29  } catch (GeneralClassicException &e) {
30  parsing_passed = interpreteLine<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  interpreteLine<double, double, int>(file, rbegin_m, rend_m, num_gridpr_m);
49  parsing_passed = parsing_passed &&
50  interpreteLine<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  interpreteLine<double, double, int>(file, zbegin_m, zend_m, num_gridpz_m);
55  parsing_passed = parsing_passed &&
56  interpreteLine<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 && interpreteLine<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 /= 100.;
80  rend_m /= 100.;
81  zbegin_m /= 100.;
82  zend_m /= 100.;
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 void FM2DMagnetoStatic::readMap() {
103  if(FieldstrengthBz_m == NULL) {
104  // declare variables and allocate memory
105  std::ifstream in;
106  std::string tmpString;
107 
108  FieldstrengthBz_m = new double[num_gridpz_m * num_gridpr_m];
109  FieldstrengthBr_m = new double[num_gridpz_m * num_gridpr_m];
110 
111  // read in and parse field map
112  in.open(Filename_m.c_str());
113  getLine(in, tmpString);
114  getLine(in, tmpString);
115  getLine(in, tmpString);
116 
117  if(swap_m) {
118  for(int i = 0; i < num_gridpz_m; i++) {
119  for(int j = 0; j < num_gridpr_m; j++) {
120  interpreteLine<double, double>(in,
121  FieldstrengthBr_m[i + j * num_gridpz_m],
122  FieldstrengthBz_m[i + j * num_gridpz_m]);
123  }
124  }
125  } else {
126  for(int j = 0; j < num_gridpr_m; j++) {
127  for(int i = 0; i < num_gridpz_m; i++) {
128  interpreteLine<double, double>(in,
129  FieldstrengthBz_m[i + j * num_gridpz_m],
130  FieldstrengthBr_m[i + j * num_gridpz_m]);
131  }
132  }
133  }
134  in.close();
135 
136  if (normalize_m) {
137  double Bzmax = 0.0;
138  // find maximum field
139  for(int i = 0; i < num_gridpz_m; ++ i) {
140  if(std::abs(FieldstrengthBz_m[i]) > Bzmax) {
141  Bzmax = std::abs(FieldstrengthBz_m[i]);
142  }
143  }
144 
145  // normalize field
146  for(int i = 0; i < num_gridpr_m * num_gridpz_m; ++ i) {
147  FieldstrengthBz_m[i] /= Bzmax;
148  FieldstrengthBr_m[i] /= Bzmax;
149  }
150  }
151 
152  INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << endl);
153  }
154 }
155 
156 void FM2DMagnetoStatic::freeMap() {
157  if(FieldstrengthBz_m != NULL) {
158  delete[] FieldstrengthBz_m;
159  delete[] FieldstrengthBr_m;
160 
161  FieldstrengthBz_m = NULL;
162  FieldstrengthBr_m = NULL;
163 
164  INFOMSG(level3 << typeset_msg("freed fieldmap '" + Filename_m + "'", "info") << endl);
165  }
166 }
167 
168 bool FM2DMagnetoStatic::getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const {
169  // do bi-linear interpolation
170  const double RR = sqrt(R(0) * R(0) + R(1) * R(1));
171 
172  const int indexr = abs((int)floor(RR / hr_m));
173  const double leverr = (RR / hr_m) - indexr;
174 
175  const int indexz = abs((int)floor((R(2)) / hz_m));
176  const double leverz = (R(2) / hz_m) - indexz;
177 
178  if((indexz < 0) || (indexz + 2 > num_gridpz_m))
179  return false;
180 
181  if(indexr + 2 > num_gridpr_m)
182  return true;
183 
184  const int index1 = indexz + indexr * num_gridpz_m;
185  const int index2 = index1 + num_gridpz_m;
186  const double BfieldR = (1.0 - leverz) * (1.0 - leverr) * FieldstrengthBr_m[index1]
187  + leverz * (1.0 - leverr) * FieldstrengthBr_m[index1 + 1]
188  + (1.0 - leverz) * leverr * FieldstrengthBr_m[index2]
189  + leverz * leverr * FieldstrengthBr_m[index2 + 1];
190 
191  const double BfieldZ = (1.0 - leverz) * (1.0 - leverr) * FieldstrengthBz_m[index1]
192  + leverz * (1.0 - leverr) * FieldstrengthBz_m[index1 + 1]
193  + (1.0 - leverz) * leverr * FieldstrengthBz_m[index2]
194  + leverz * leverr * FieldstrengthBz_m[index2 + 1];
195 
196  if(RR != 0) {
197  B(0) += BfieldR * R(0) / RR;
198  B(1) += BfieldR * R(1) / RR;
199  }
200  B(2) += BfieldZ;
201  return false;
202 }
203 
204 bool FM2DMagnetoStatic::getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const {
205 
206  double BfieldR, BfieldZ;
207 
208  const double RR = sqrt(R(0) * R(0) + R(1) * R(1));
209 
210  const int indexr = (int)floor(RR / hr_m);
211  const double leverr = (RR / hr_m) - indexr;
212 
213  const int indexz = (int)floor((R(2)) / hz_m);
214  const double leverz = (R(2) / hz_m) - indexz;
215 
216  if((indexz < 0) || (indexz + 2 > num_gridpz_m))
217  return false;
218 
219  if(indexr + 2 > num_gridpr_m)
220  return true;
221 
222  const int index1 = indexz + indexr * num_gridpz_m;
223  const int index2 = index1 + num_gridpz_m;
224  if(dir == DZ) {
225  if(indexz > 0) {
226  if(indexz < num_gridpz_m - 1) {
227  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)
228  + (FieldstrengthBr_m[index1 + 1] - FieldstrengthBr_m[index1 - 1]) / (2.*hz_m))
229  + leverr * ((FieldstrengthBr_m[index2 - 1] - 2. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index2 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
230  + (FieldstrengthBr_m[index2 + 1] - FieldstrengthBr_m[index2 - 1]) / (2.*hz_m));
231 
232  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)
233  + (FieldstrengthBz_m[index1 + 1] - FieldstrengthBz_m[index1 - 1]) / (2.*hz_m))
234  + leverr * ((FieldstrengthBz_m[index2 - 1] - 2. * FieldstrengthBz_m[index2] + FieldstrengthBz_m[index2 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
235  + (FieldstrengthBz_m[index2 + 1] - FieldstrengthBz_m[index2 - 1]) / (2.*hz_m));
236  } else {
237  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)
238  - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBr_m[index1 - 2]) / (2.*hz_m))
239  + leverr * ((FieldstrengthBr_m[index2] - 2. * FieldstrengthBr_m[index2 - 1] + FieldstrengthBr_m[index2 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
240  - (3. * FieldstrengthBr_m[index2] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBr_m[index2 - 2]) / (2.*hz_m));
241 
242  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)
243  - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBz_m[index1 - 2]) / (2.*hz_m))
244  + leverr * ((FieldstrengthBz_m[index2] - 2. * FieldstrengthBz_m[index2 - 1] + FieldstrengthBz_m[index2 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
245  - (3. * FieldstrengthBz_m[index2] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBz_m[index2 - 2]) / (2.*hz_m));
246  }
247  } else {
248  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)
249  - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index1 + 2]) / (2.*hz_m))
250  + leverr * ((FieldstrengthBr_m[index2] - 2. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBr_m[index2 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
251  - (3. * FieldstrengthBr_m[index2] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index2 + 2]) / (2.*hz_m));
252 
253  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)
254  - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBz_m[index1 + 2]) / (2.*hz_m))
255  + leverr * ((FieldstrengthBz_m[index2] - 2. * FieldstrengthBz_m[index2 + 1] + FieldstrengthBz_m[index2 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
256  - (3. * FieldstrengthBz_m[index2] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBz_m[index2 + 2]) / (2.*hz_m));
257  }
258  } else {
259  if(indexr > 0) {
260  const int index_1 = index1 - num_gridpz_m;
261  if(indexr < num_gridpr_m - 1) {
262  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)
263  + (FieldstrengthBr_m[index2 + 1] - FieldstrengthBr_m[index_1 + 1]) / (2.*hr_m))
264  + leverz * ((FieldstrengthBr_m[index_1] - 2. * FieldstrengthBr_m[index1] + FieldstrengthBr_m[index2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
265  + (FieldstrengthBr_m[index2] - FieldstrengthBr_m[index_1]) / (2.*hr_m));
266 
267  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)
268  + (FieldstrengthBz_m[index2 + 1] - FieldstrengthBz_m[index_1 + 1]) / (2.*hr_m))
269  + leverz * ((FieldstrengthBz_m[index_1] - 2. * FieldstrengthBz_m[index1] + FieldstrengthBz_m[index2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
270  + (FieldstrengthBz_m[index2] - FieldstrengthBz_m[index_1]) / (2.*hr_m));
271  } else {
272  const int index_2 = index_1 - num_gridpz_m;
273  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)
274  - (3. * FieldstrengthBr_m[index1 + 1] - 4. * FieldstrengthBr_m[index_1 + 1] + FieldstrengthBr_m[index_2 + 1]) / (2.*hr_m))
275  + leverz * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index_1] + FieldstrengthBr_m[index_2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
276  - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index_1] + FieldstrengthBr_m[index_2]) / (2.*hr_m));
277 
278  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)
279  - (3. * FieldstrengthBz_m[index1 + 1] - 4. * FieldstrengthBr_m[index_1 + 1] + FieldstrengthBz_m[index_2 + 1]) / (2.*hr_m))
280  + leverz * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index_1] + FieldstrengthBz_m[index_2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
281  - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index_1] + FieldstrengthBz_m[index_2]) / (2.*hr_m));
282  }
283  } else {
284  const int index3 = index2 + num_gridpz_m;
285  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)
286  - (3. * FieldstrengthBr_m[index1 + 1] - 4. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBr_m[index3 + 1]) / (2.*hr_m))
287  + leverz * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index3]) * (RR - indexr * hr_m) / (hr_m * hr_m)
288  - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index3]) / (2.*hr_m));
289 
290  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)
291  - (3. * FieldstrengthBz_m[index1 + 1] - 4. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBz_m[index3 + 1]) / (2.*hr_m))
292  + leverz * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index2] + FieldstrengthBz_m[index3]) * (RR - indexr * hr_m) / (hr_m * hr_m)
293  - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index2] + FieldstrengthBz_m[index3]) / (2.*hr_m));
294  }
295 
296  }
297 
298  if(RR > 1.e-12) {
299  B(0) += BfieldR * R(0) / RR;
300  B(1) += BfieldR * R(1) / RR;
301  }
302  B(2) += BfieldZ;
303  return false;
304 }
305 
306 void FM2DMagnetoStatic::getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const {
307  zBegin = zbegin_m;
308  zEnd = zend_m;
309  rBegin = rbegin_m;
310  rEnd = rend_m;
311 }
312 
313 void FM2DMagnetoStatic::getFieldDimensions(double &xIni, double &xFinal, double &yIni, double &yFinal, double &zIni, double &zFinal) const {}
314 
315 void FM2DMagnetoStatic::swap() {
316  if(swap_m) swap_m = false;
317  else swap_m = true;
318 }
319 
320 void FM2DMagnetoStatic::getInfo(Inform *msg) {
321  (*msg) << Filename_m << " (2D magnetostatic); zini= " << zbegin_m << " m; zfinal= " << zend_m << " m;" << endl;
322 }
323 
324 double FM2DMagnetoStatic::getFrequency() const {
325  return 0.0;
326 }
327 
328 void FM2DMagnetoStatic::setFrequency(double freq)
329 { ;}
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
Physics::e
constexpr double e
The value of .
Definition: Physics.h:40
FM2DMagnetoStatic::getFrequency
virtual double getFrequency() const
Definition: FM2DMagnetoStatic.cpp:324
DiffDirection
DiffDirection
Definition: Fieldmap.h:54
FM2DMagnetoStatic::getFieldDimensions
virtual void getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const
Definition: FM2DMagnetoStatic.cpp:306
gmsg
Inform * gmsg
Definition: Main.cpp:21
Util::toUpper
std::string toUpper(const std::string &str)
Definition: Util.cpp:130
Fieldmap::disableFieldmapWarning
void disableFieldmapWarning()
Definition: Fieldmap.cpp:610
FM2DMagnetoStatic::getFieldDerivative
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
Definition: FM2DMagnetoStatic.cpp:204
Fieldmap::typeset_msg
static std::string typeset_msg(const std::string &msg, const std::string &title)
Definition: Fieldmap.cpp:647
Fieldmap::noFieldmapWarning
void noFieldmapWarning()
Definition: Fieldmap.cpp:618
FM2DMagnetoStatic::freeMap
virtual void freeMap()
Definition: FM2DMagnetoStatic.cpp:156
Fieldmap::normalize_m
bool normalize_m
Definition: Fieldmap.h:113
Fieldmap::lines_read_m
int lines_read_m
Definition: Fieldmap.h:111
FM2DMagnetoStatic::getInfo
virtual void getInfo(Inform *msg)
Definition: FM2DMagnetoStatic.cpp:320
Fieldmap::interpreteEOF
bool interpreteEOF(std::ifstream &in)
Definition: Fieldmap.cpp:547
INFOMSG
#define INFOMSG(msg)
Definition: IpplInfo.h:397
FM2DMagnetoStatic::getFieldstrength
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
Definition: FM2DMagnetoStatic.cpp:168
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
FM2DMagnetoStatic::setFrequency
virtual void setFrequency(double freq)
Definition: FM2DMagnetoStatic.cpp:328
DZ
Definition: Fieldmap.h:57
Fieldmap::getLine
void getLine(std::ifstream &in, std::string &buffer)
Definition: Fieldmap.h:198
Fieldmap::Filename_m
std::string Filename_m
Definition: Fieldmap.h:110
Fieldmap::Type
MapType Type
Definition: Fieldmap.h:107
level3
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
Inform
Definition: Inform.h:41
FM2DMagnetoStatic::FM2DMagnetoStatic
FM2DMagnetoStatic(std::string aFilename)
Definition: FM2DMagnetoStatic.cpp:11
FM2DMagnetoStatic::readMap
virtual void readMap()
Definition: FM2DMagnetoStatic.cpp:102
floor
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Definition: PETE.h:816
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
Generated on Sat Jan 4 2020 14:47:54 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5