OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
Astra1D_fast.cpp
Go to the documentation of this file.
1 
2 #include "Fields/Astra1D_fast.h"
3 #include "Fields/Fieldmap.hpp"
4 #include "Physics/Physics.h"
5 #include "gsl/gsl_fft_real.h"
6 
7 #include <fstream>
8 #include <ios>
9 
10 _Astra1D_fast::_Astra1D_fast(const std::string& filename):
11  _Fieldmap(filename) {
12  normalize_m = true;
13 }
14 
16  freeMap();
17 }
18 
20 { }
21 
23  if(onAxisField_m != nullptr) {
24  for(int i = 0; i < 4; ++i) {
25  gsl_spline_free(onAxisInterpolants_m[i]);
26  gsl_interp_accel_free(onAxisAccel_m[i]);
27  }
28 
29  delete[] onAxisField_m;
30  onAxisField_m = nullptr;
31  delete[] zvals_m;
32  zvals_m = nullptr;
33  }
34 }
35 
36 void _Astra1D_fast::getOnaxisEz(std::vector<std::pair<double, double> > & /*F*/)
37 { }
38 
39 bool _Astra1D_fast::determineNumSamplingPoints(std::ifstream &file) {
40  double tmpDouble, tmpDouble2;
41  unsigned int skippedValues = 0;
42  bool passed = interpretLine<double, double>(file, zbegin_m, tmpDouble);
43 
44  tmpDouble2 = zbegin_m;
45  while (!file.eof() && passed) {
46  passed = interpretLine<double, double>(file, zend_m, tmpDouble, false);
47  if (zend_m - tmpDouble2 > 1e-10) {
48  tmpDouble2 = zend_m;
49  } else if (passed) {
50  ++ skippedValues;
51  }
52  }
53 
54  num_gridpz_m = lines_read_m - numHeaderLines_m - skippedValues;
55  lines_read_m = 0;
56 
57  return passed || file.eof();
58 }
59 
60 double _Astra1D_fast::readFieldData(std::ifstream &file) {
61  double tmpDouble = zbegin_m - hz_m;
62  double maxValue = 0;
63  int nsp;
64 
65  for (nsp = 0; nsp < num_gridpz_m; ) {
66  interpretLine<double, double>(file, zvals_m[nsp], onAxisField_m[nsp]);
67 
68  if(zvals_m[nsp] - tmpDouble > 1e-10) {
69  if(std::abs(onAxisField_m[nsp]) > maxValue) {
70  maxValue = std::abs(onAxisField_m[nsp]);
71  }
72  tmpDouble = zvals_m[nsp];
73  ++ nsp;
74  }
75  }
76  num_gridpz_m = nsp;
77  hz_m = (zend_m - zbegin_m) / (num_gridpz_m - 1);
78 
79  if (!normalize_m) return 1.0;
80 
81  return maxValue;
82 }
83 
84 void _Astra1D_fast::normalizeFieldData(double maxValue) {
85  int ii = num_gridpz_m - 1;
86  for(int i = 0; i < num_gridpz_m; ++ i, -- ii) {
87  onAxisField_m[i] /= maxValue;
88  zvals_m[ii] -= zvals_m[0];
89  }
90 }
91 
93  std::vector<double> zvals(num_gridpz_m, 0.0);
94 
95  for(int i = 1; i < num_gridpz_m - 1; ++ i) {
96  zvals[i] = zvals[i - 1] + hz_m;
97  }
98  zvals[num_gridpz_m - 1] = zvals_m[num_gridpz_m - 1];
99 
100  return zvals;
101 }
102 
103 std::vector<double> _Astra1D_fast::interpolateFieldData(std::vector<double> &samplingPoints) {
104  std::vector<double> realValues(num_gridpz_m);
105 
106  onAxisInterpolants_m[0] = gsl_spline_alloc(gsl_interp_cspline, num_gridpz_m);
107  onAxisAccel_m[0] = gsl_interp_accel_alloc();
108 
110 
111  for (int i = 0; i < num_gridpz_m - 1; ++ i) {
112  realValues[i] = gsl_spline_eval(onAxisInterpolants_m[0], samplingPoints[i], onAxisAccel_m[0]);
113  }
114  realValues[num_gridpz_m - 1] = onAxisField_m[num_gridpz_m - 1];
115 
116  return realValues;
117 }
118 
119 std::vector<double> _Astra1D_fast::computeFourierCoefficients(int accuracy,
120  std::vector<double> &evenFieldSampling) {
121  std::vector<double> fourierCoefficients(2 * accuracy - 1);
122  std::vector<double> RealValues(2 * num_gridpz_m, 0.0);
123  int ii = num_gridpz_m, iii = num_gridpz_m - 1;
124 
125  for(int i = 0; i < num_gridpz_m; ++ i, ++ ii, -- iii) {
126  RealValues[ii] = evenFieldSampling[i];
127  RealValues[iii] = evenFieldSampling[i];
128  }
129 
130  gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(2 * num_gridpz_m);
131  gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(2 * num_gridpz_m);
132 
133  gsl_fft_real_transform(&RealValues[0], 1, 2 * num_gridpz_m, real, work);
134 
135  gsl_fft_real_workspace_free(work);
136  gsl_fft_real_wavetable_free(real);
137 
138  // normalize to Ez_max = 1 MV/m
139  fourierCoefficients[0] = RealValues[0] / (2 * num_gridpz_m);
140  for(int i = 1; i < 2 * accuracy - 1; i++) {
141  fourierCoefficients[i] = RealValues[i] / num_gridpz_m;
142  }
143 
144  return fourierCoefficients;
145 }
146 
147 void _Astra1D_fast::computeFieldDerivatives(std::vector<double> & fourierComponents, int accuracy) {
148  double interiorDerivative, base;
149  double coskzl, sinkzl, z = 0.0;
150  std::vector<double> zvals(num_gridpz_m);
151  std::vector<double> higherDerivatives[3];
152  higherDerivatives[0].resize(num_gridpz_m, 0.0);
153  higherDerivatives[1].resize(num_gridpz_m, 0.0);
154  higherDerivatives[2].resize(num_gridpz_m, 0.0);
155 
156  for(int i = 0; i < num_gridpz_m; ++ i, z += hz_m) {
157  zvals[i] = z;
158  const double kz = Physics::two_pi * (zvals[i] / length_m + 0.5);
159 
160  for(int l = 1; l < accuracy; ++l) {
161  int coefIndex = 2 * l - 1;
162  base = Physics::two_pi / length_m * l;
163  interiorDerivative = base;
164  coskzl = cos(kz * l);
165  sinkzl = sin(kz * l);
166 
167  higherDerivatives[0][i] += interiorDerivative * (-fourierComponents[coefIndex] * sinkzl
168  - fourierComponents[coefIndex + 1] * coskzl);
169  interiorDerivative *= base;
170  higherDerivatives[1][i] += interiorDerivative * (-fourierComponents[coefIndex] * coskzl
171  + fourierComponents[coefIndex + 1] * sinkzl);
172  interiorDerivative *= base;
173  higherDerivatives[2][i] += interiorDerivative * (fourierComponents[coefIndex] * sinkzl
174  + fourierComponents[coefIndex + 1] * coskzl);
175  }
176  }
177 
178  for(int i = 1; i < 4; ++i) {
179  onAxisInterpolants_m[i] = gsl_spline_alloc(gsl_interp_cspline, num_gridpz_m);
180  onAxisAccel_m[i] = gsl_interp_accel_alloc();
181  gsl_spline_init(onAxisInterpolants_m[i], &zvals[0], &higherDerivatives[i - 1][0], num_gridpz_m);
182  }
183 }
virtual ~_Astra1D_fast()
int numHeaderLines_m
Definition: Astra1D_fast.h:49
constexpr double two_pi
The value of .
Definition: Physics.h:33
double length_m
Definition: Astra1D_fast.h:46
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
double zbegin_m
Definition: Astra1D_fast.h:44
int lines_read_m
Definition: Fieldmap.h:119
gsl_spline * onAxisInterpolants_m[4]
Definition: Astra1D_fast.h:39
std::vector< double > interpolateFieldData(std::vector< double > &samplingPoints)
double * onAxisField_m
Definition: Astra1D_fast.h:37
virtual void freeMap()
void normalizeFieldData(double maxEz)
virtual void getOnaxisEz(std::vector< std::pair< double, double > > &F)
double readFieldData(std::ifstream &file)
bool normalize_m
Definition: Fieldmap.h:121
std::vector< double > getEvenlyDistributedSamplingPoints()
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition: TpsMath.h:129
_Astra1D_fast(const std::string &filename)
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
gsl_interp_accel * onAxisAccel_m[4]
Definition: Astra1D_fast.h:40
constexpr double e
The value of .
Definition: Physics.h:39
void computeFieldDerivatives(std::vector< double > &fourierComponents, int accuracy)
std::vector< double > computeFourierCoefficients(int accuracy, std::vector< double > &evenSampling)
Tps< T > sin(const Tps< T > &x)
Sine.
Definition: TpsMath.h:111
virtual void readMap()
double * zvals_m
Definition: Astra1D_fast.h:38
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
bool determineNumSamplingPoints(std::ifstream &file)