OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
Astra1DDynamic.cpp
Go to the documentation of this file.
2 #include "Fields/Fieldmap.hpp"
3 #include "Physics/Physics.h"
4 #include "Physics/Units.h"
6 #include "Utilities/Util.h"
7 
8 #include "gsl/gsl_interp.h"
9 #include "gsl/gsl_spline.h"
10 #include "gsl/gsl_fft_real.h"
11 
12 #include <fstream>
13 #include <ios>
14 
15 
16 _Astra1DDynamic::_Astra1DDynamic(const std::string& filename):
17  _Fieldmap(filename),
18  FourCoefs_m(nullptr) {
19 
20  std::ifstream file;
21  int skippedValues = 0;
22  std::string tmpString;
23  double tmpDouble;
24 
26 
27  // open field map, parse it and disable element on error
28  file.open(Filename_m.c_str());
29  if (file.good()) {
30  bool parsing_passed = true;
31  try {
32  parsing_passed = interpretLine<std::string, int>(file,
33  tmpString,
34  accuracy_m);
35  } catch (GeneralClassicException &e) {
36  parsing_passed = interpretLine<std::string, int, std::string>(file,
37  tmpString,
38  accuracy_m,
39  tmpString);
40 
41  tmpString = Util::toUpper(tmpString);
42  if (tmpString != "TRUE" &&
43  tmpString != "FALSE")
44  throw GeneralClassicException("_Astra1DDynamic::_Astra1DDynamic",
45  "The third string on the first line of 1D field "
46  "maps has to be either TRUE or FALSE");
47 
48  normalize_m = (tmpString == "TRUE");
49  }
50 
51  parsing_passed = parsing_passed &&
52  interpretLine<double>(file, frequency_m);
53  parsing_passed = parsing_passed &&
54  interpretLine<double, double>(file, zbegin_m, tmpDouble);
55 
56  double tmpDouble2 = zbegin_m;
57  while(!file.eof() && parsing_passed) {
58  parsing_passed = interpretLine<double, double>(file, zend_m, tmpDouble, false);
59  if (zend_m - tmpDouble2 > 1e-10) {
60  tmpDouble2 = zend_m;
61  } else if (parsing_passed) {
62  ++ skippedValues;
63  }
64  }
65 
66  num_gridpz_m = lines_read_m - 3 - skippedValues;
67  lines_read_m = 0;
68 
69  if (!parsing_passed && !file.eof()) {
71  zend_m = zbegin_m - 1e-3;
72  throw GeneralClassicException("_Astra1DDynamic::_Astra1DDynamic",
73  "An error occured when reading the fieldmap '" + Filename_m + "'");
74  } else {
75  // conversion from MHz to Hz and from frequency to angular frequency
78  }
79  length_m = 2.0 * num_gridpz_m * (zend_m - zbegin_m) / (num_gridpz_m - 1);
80  file.close();
81  } else {
83  zbegin_m = 0.0;
84  zend_m = -1e-3;
85  }
86 }
87 
89  freeMap();
90 }
91 
92 Astra1DDynamic _Astra1DDynamic::create(const std::string& filename)
93 {
94  return Astra1DDynamic(new _Astra1DDynamic(filename));
95 }
96 
98  if (FourCoefs_m == nullptr) {
99  // declare variables and allocate memory
100  std::ifstream in;
101 
102  bool parsing_passed = true;
103 
104  std::string tmpString;
105 
106  double tmpDouble;
107  double Ez_max = 0.0;
108  double dz = (zend_m - zbegin_m) / (num_gridpz_m - 1);
109 
110  double *RealValues = new double[2 * num_gridpz_m];
111  double *zvals = new double[num_gridpz_m];
112 
113  gsl_spline *Ez_interpolant = gsl_spline_alloc(gsl_interp_cspline, num_gridpz_m);
114  gsl_interp_accel *Ez_accel = gsl_interp_accel_alloc();
115 
116  gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(2 * num_gridpz_m);
117  gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(2 * num_gridpz_m);
118 
119  FourCoefs_m = new double[2 * accuracy_m - 1];
120 
121  // read in and parse field map
122  in.open(Filename_m.c_str());
123  getLine(in, tmpString);
124  getLine(in, tmpString);
125 
126  tmpDouble = zbegin_m - dz;
127  for (int i = 0; i < num_gridpz_m && parsing_passed; /* skip increment of i here */) {
128  parsing_passed = interpretLine<double, double>(in, zvals[i], RealValues[i]);
129  // the sequence of z-position should be strictly increasing
130  // drop sampling points that don't comply to this
131  if (zvals[i] - tmpDouble > 1e-10) {
132  if (std::abs(RealValues[i]) > Ez_max) {
133  Ez_max = std::abs(RealValues[i]);
134  }
135  tmpDouble = zvals[i];
136  ++ i; // increment i only if sampling point is accepted
137  }
138  }
139  in.close();
140 
141  gsl_spline_init(Ez_interpolant, zvals, RealValues, num_gridpz_m);
142 
143  // get equidistant sampling from the, possibly, non-equidistant sampling
144  // using cubic spline for this
145  int ii = num_gridpz_m;
146  for (int i = 0; i < num_gridpz_m - 1; ++ i, ++ ii) {
147  double z = zbegin_m + dz * i;
148  RealValues[ii] = gsl_spline_eval(Ez_interpolant, z, Ez_accel);
149  }
150  RealValues[ii ++] = RealValues[num_gridpz_m - 1];
151  // prepend mirror sampling points such that field values are periodic for sure
152  -- ii; // ii == 2*num_gridpz_m at the moment
153  for (int i = 0; i < num_gridpz_m; ++ i, -- ii) {
154  RealValues[i] = RealValues[ii];
155  }
156 
157  gsl_fft_real_transform(RealValues, 1, 2 * num_gridpz_m, real, work);
158 
159  if (!normalize_m) {
160  Ez_max = 1.0;
161  }
162 
163  FourCoefs_m[0] = RealValues[0] / (Ez_max * Units::Vpm2MVpm* 2. * num_gridpz_m);
164  for (int i = 1; i < 2 * accuracy_m - 1; ++ i) {
165  FourCoefs_m[i] = RealValues[i] / (Ez_max * Units::Vpm2MVpm * num_gridpz_m);
166  }
167 
168  gsl_spline_free(Ez_interpolant);
169  gsl_interp_accel_free(Ez_accel);
170 
171  gsl_fft_real_workspace_free(work);
172  gsl_fft_real_wavetable_free(real);
173 
174  delete[] zvals;
175  delete[] RealValues;
176 
177  INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << endl);
178  }
179 }
180 
182  if (FourCoefs_m != nullptr) {
183  delete[] FourCoefs_m;
184  FourCoefs_m = nullptr;
185  }
186 }
187 
189  // do fourier interpolation in z-direction
190  const double RR2 = R(0) * R(0) + R(1) * R(1);
191 
192  const double kz = Physics::two_pi * (R(2) - zbegin_m) / length_m + Physics::pi;
193 
194  double ez = FourCoefs_m[0];
195  double ezp = 0.0;
196  double ezpp = 0.0;
197  double ezppp = 0.0;
198 
199  int n = 1;
200  for (int l = 1; l < accuracy_m ; ++ l, n += 2) {
201  double somefactor_base = Physics::two_pi / length_m * l; // = \frac{d(kz*l)}{dz}
202  double somefactor = 1.0;
203  double coskzl = cos(kz * l);
204  double sinkzl = sin(kz * l);
205  ez += (FourCoefs_m[n] * coskzl - FourCoefs_m[n + 1] * sinkzl);
206  somefactor *= somefactor_base;
207  ezp += somefactor * (-FourCoefs_m[n] * sinkzl - FourCoefs_m[n + 1] * coskzl);
208  somefactor *= somefactor_base;
209  ezpp += somefactor * (-FourCoefs_m[n] * coskzl + FourCoefs_m[n + 1] * sinkzl);
210  somefactor *= somefactor_base;
211  ezppp += somefactor * (FourCoefs_m[n] * sinkzl + FourCoefs_m[n + 1] * coskzl);
212  }
213  // expand the field off-axis
214  const double f = -(ezpp + ez * xlrep_m * xlrep_m) / 16.;
215  const double fp = -(ezppp + ezp * xlrep_m * xlrep_m) / 16.;
216 
217  const double EfieldR = -(ezp / 2. + fp * RR2);
218  const double BfieldT = (ez / 2. + f * RR2) * xlrep_m / Physics::c;
219 
220  E(0) += EfieldR * R(0);
221  E(1) += EfieldR * R(1);
222  E(2) += ez + 4. * f * RR2;
223  B(0) += -BfieldT * R(1);
224  B(1) += BfieldT * R(0);
225 
226  return false;
227 }
228 
229 bool _Astra1DDynamic::getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &/*B*/, const DiffDirection &/*dir*/) const {
230  const double kz = Physics::two_pi * (R(2) - zbegin_m) / length_m + Physics::pi;
231  double ezp = 0.0;
232 
233  int n = 1;
234  for (int l = 1; l < accuracy_m; ++ l, n += 2)
235  ezp += Physics::two_pi / length_m * l * (-FourCoefs_m[n] * sin(kz * l) - FourCoefs_m[n + 1] * cos(kz * l));
236 
237  E(2) += ezp;
238 
239  return false;
240 }
241 
242 void _Astra1DDynamic::getFieldDimensions(double &zBegin, double &zEnd) const {
243  zBegin = zbegin_m;
244  zEnd = zend_m;
245 }
246 
247 void _Astra1DDynamic::getFieldDimensions(double &/*xIni*/, double &/*xFinal*/, double &/*yIni*/, double &/*yFinal*/, double &/*zIni*/, double &/*zFinal*/) const {}
248 
250 { }
251 
253  (*msg) << Filename_m << " (1D dynamic); zini= " << zbegin_m << " m; zfinal= " << zend_m << " m;" << endl;
254 }
255 
257  return frequency_m;
258 }
259 
260 void _Astra1DDynamic::setFrequency(double freq) {
261  frequency_m = freq;
262 }
263 
264 void _Astra1DDynamic::getOnaxisEz(std::vector<std::pair<double, double> > & F) {
265  double Ez_max = 0.0;
266  double tmpDouble;
267  int tmpInt;
268  std::string tmpString;
269  F.resize(num_gridpz_m);
270 
271  std::ifstream in(Filename_m.c_str());
272  interpretLine<std::string, int>(in, tmpString, tmpInt);
273  interpretLine<double>(in, tmpDouble);
274 
275  for (int i = 0; i < num_gridpz_m; ++ i) {
276  interpretLine<double, double>(in, F[i].first, F[i].second);
277  if (std::abs(F[i].second) > Ez_max) {
278  Ez_max = std::abs(F[i].second);
279  }
280  }
281  in.close();
282 
283  for (int i = 0; i < num_gridpz_m; ++ i) {
284  F[i].second /= Ez_max;
285  }
286 }
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:45
MapType Type
Definition: Fieldmap.h:115
DiffDirection
Definition: Fieldmap.h:55
constexpr double two_pi
The value of .
Definition: Physics.h:33
void noFieldmapWarning()
Definition: Fieldmap.cpp:618
std::shared_ptr< _Astra1DDynamic > Astra1DDynamic
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
_Astra1DDynamic(const std::string &filename)
double *__restrict__ FourCoefs_m
virtual void getOnaxisEz(std::vector< std::pair< double, double > > &F)
int lines_read_m
Definition: Fieldmap.h:119
static std::string typeset_msg(const std::string &msg, const std::string &title)
Definition: Fieldmap.cpp:649
constexpr double pi
The value of .
Definition: Physics.h:30
std::string toUpper(const std::string &str)
Definition: Util.cpp:147
constexpr double Vpm2MVpm
Definition: Units.h:125
#define INFOMSG(msg)
Definition: IpplInfo.h:348
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
virtual double getFrequency() const
virtual void setFrequency(double freq)
static Astra1DDynamic create(const std::string &filename)
virtual void getFieldDimensions(double &zBegin, double &zEnd) const
bool normalize_m
Definition: Fieldmap.h:121
Definition: Inform.h:42
virtual void swap()
virtual void getInfo(Inform *)
virtual void readMap()
void disableFieldmapWarning()
Definition: Fieldmap.cpp:610
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition: TpsMath.h:129
constexpr double MHz2Hz
Definition: Units.h:113
and that you know you can do these things To protect your we need to make restrictions that forbid anyone to deny you these rights or to ask you to surrender the rights These restrictions translate to certain responsibilities for you if you distribute copies of the or if you modify it For if you distribute copies of such a whether gratis or for a you must give the recipients all the rights that you have You must make sure that receive or can get the source code And you must show them these terms so they know their rights We protect your rights with two distribute and or modify the software for each author s protection and we want to make certain that everyone understands that there is no warranty for this free software If the software is modified by someone else and passed we want its recipients to know that what they have is not the so that any problems introduced by others will not reflect on the original authors reputations any free program is threatened constantly by software patents We wish to avoid the danger that redistributors of a free program will individually obtain patent in effect making the program proprietary To prevent we have made it clear that any patent must be licensed for everyone s free use or not licensed at all The precise terms and conditions for distribution and modification follow GNU GENERAL PUBLIC LICENSE TERMS AND CONDITIONS FOR DISTRIBUTION AND MODIFICATION This License applies to any program or other work which contains a notice placed by the copyright holder saying it may be distributed under the terms of this General Public License The refers to any such program or work
Definition: LICENSE:43
void getLine(std::ifstream &in, std::string &buffer)
Definition: Fieldmap.h:122
virtual ~_Astra1DDynamic()
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
constexpr double e
The value of .
Definition: Physics.h:39
std::string Filename_m
Definition: Fieldmap.h:118
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
Tps< T > sin(const Tps< T > &x)
Sine.
Definition: TpsMath.h:111
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
virtual void freeMap()