7 #include "gsl/gsl_fft_real.h"
21 if (fieldFile.good()) {
57 std::vector<double> fourierCoefs
65 onAxisFieldPP, onAxisFieldPPP);
71 delete [] onAxisFieldP;
72 delete [] onAxisFieldPP;
73 delete [] onAxisFieldPPP;
102 std::vector<double> fieldComponents;
123 double &rBegin,
double &rEnd)
const {
130 double &yIni,
double &yFinal,
131 double &zIni,
double &zFinal)
const {}
138 <<
" (1D electrotostatic); zini= "
151 bool parsingPassed) {
155 parsingPassed = parsingPassed
156 && interpreteLine<double>(fieldFile, tempDouble);
163 double onAxisFieldP[],
164 double onAxisFieldPP[],
165 double onAxisFieldPPP[]) {
171 onAxisFieldP[zStepIndex] = 0.0;
172 onAxisFieldPP[zStepIndex] = 0.0;
173 onAxisFieldPPP[zStepIndex] = 0.0;
179 double coskzn =
cos(kz *
n);
180 double sinkzn =
sin(kz * n);
182 onAxisFieldP[zStepIndex] += kn * (-fourierCoefs.at(coefIndex) * sinkzn
183 - fourierCoefs.at(coefIndex + 1) * coskzn);
185 double derivCoeff =
pow(kn, 2.0);
186 onAxisFieldPP[zStepIndex] += derivCoeff * (-fourierCoefs.at(coefIndex) * coskzn
187 + fourierCoefs.at(coefIndex + 1) * sinkzn);
189 onAxisFieldPPP[zStepIndex] += derivCoeff * (fourierCoefs.at(coefIndex) * sinkzn
190 + fourierCoefs.at(coefIndex + 1) * coskzn);
198 std::vector<double> fieldComponents)
const {
200 double radiusSq =
pow(
R(0), 2.0) +
pow(
R(1), 2.0);
201 double transverseEFactor = -fieldComponents.at(1) / 2.0
202 + radiusSq * fieldComponents.at(3) / 16.0;
204 E(0) +=
R(0) * transverseEFactor;
205 E(1) +=
R(1) * transverseEFactor;
206 E(2) += fieldComponents.at(0) - fieldComponents.at(2) * radiusSq / 4.0;
211 std::vector<double> &fieldComponents)
const {
226 gsl_fft_real_wavetable *waveTable = gsl_fft_real_wavetable_alloc(totalSize);
227 gsl_fft_real_workspace *workSpace = gsl_fft_real_workspace_alloc(totalSize);
230 double *fieldDataReflected =
new double[totalSize];
232 fieldDataReflected[numberOfGridPoints_m - 1 + dataIndex] =
233 fieldData[dataIndex];
234 fieldDataReflected[numberOfGridPoints_m - 1 - dataIndex] =
235 fieldData[dataIndex];
238 gsl_fft_real_transform(fieldDataReflected, 1,
240 waveTable, workSpace);
242 std::vector<double> fourierCoefs;
243 fourierCoefs.push_back(fieldDataReflected[0] / totalSize);
244 for (
unsigned int coefIndex = 1; coefIndex + 1 < 2 *
accuracy_m; ++ coefIndex)
245 fourierCoefs.push_back(2.0 * fieldDataReflected[coefIndex] / totalSize);
247 delete [] fieldDataReflected;
248 gsl_fft_real_workspace_free(workSpace);
249 gsl_fft_real_wavetable_free(waveTable);
256 double onAxisFieldPP[],
257 double onAxisFieldPPP[]) {
270 z[zStepIndex] =
deltaZ_m * zStepIndex;
275 onAxisFieldP, numberOfGridPoints_m);
277 onAxisFieldPP, numberOfGridPoints_m);
279 onAxisFieldPPP, numberOfGridPoints_m);
304 for (
auto fourierIt = fourierCoefs.begin(); fourierIt < fourierCoefs.end(); ++ fourierIt)
305 *fourierIt *= 1.0e6 / maxEz;
310 double fieldData[]) {
314 interpreteLine<double>(fieldFile, fieldData[dataIndex]);
315 if (
std::abs(fieldData[dataIndex]) > maxEz)
316 maxEz =
std::abs(fieldData[dataIndex]);
327 std::string tempString;
330 bool parsingPassed =
true;
332 parsingPassed = interpreteLine<std::string, unsigned int>(fieldFile,
336 parsingPassed = interpreteLine<std::string, unsigned int, std::string>(fieldFile,
342 if (tempString !=
"TRUE" &&
343 tempString !=
"FALSE")
345 "The third string on the first line of 1D field "
346 "maps has to be either TRUE or FALSE");
350 parsingPassed = parsingPassed &&
351 interpreteLine<double, double, unsigned int>(fieldFile,
355 parsingPassed = parsingPassed &&
356 interpreteLine<double, double, int>(fieldFile,
365 return parsingPassed;
370 std::string tempString;
372 getLine(fieldFile, tempString);
373 getLine(fieldFile, tempString);
374 getLine(fieldFile, tempString);
380 zSampling[zStepIndex] =
deltaZ_m * zStepIndex;
void stripFileHeader(std::ifstream &fieldFile)
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
virtual void getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const
gsl_interp_accel * onAxisFieldPPAccel_m
constexpr double e
The value of .
Tps< T > sin(const Tps< T > &x)
Sine.
constexpr double two_pi
The value of .
gsl_spline * onAxisFieldInterpolants_m
On axis field data.
virtual void setFrequency(double freq)
std::string toUpper(const std::string &str)
std::vector< double > computeFourierCoefficients(double fieldData[])
void disableFieldmapWarning()
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
static std::string typeset_msg(const std::string &msg, const std::string &title)
void computeInterpolationVectors(double onAxisFieldP[], double onAxisFieldPP[], double onAxisFieldPPP[])
double length_m
Longitudinal end of field.
void computeFieldDerivatives(std::vector< double > fourierCoefs, double onAxisFieldP[], double onAxisFieldPP[], double onAxisFieldPPP[])
FM1DElectroStatic_fast(std::string aFilename)
~FM1DElectroStatic_fast()
constexpr double pi
The value of .
unsigned int numberOfGridPoints_m
Field length.
void computeFieldOffAxis(const Vector_t &R, Vector_t &E, Vector_t &B, std::vector< double > fieldComponents) const
double zBegin_m
Maximum radius of field.
void computeFieldOnAxis(double z, std::vector< double > &fieldComponents) const
bool interpreteEOF(std::ifstream &in)
bool readFileHeader(std::ifstream &fieldFile)
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
gsl_spline * onAxisFieldPPInterpolants_m
On axis field first derivative interpolation structure.
gsl_interp_accel * onAxisFieldPAccel_m
gsl_spline * onAxisFieldPPPInterpolants_m
On axis field second derivative interpolation structure.
double rEnd_m
Minimum radius of field.
Tps< T > cos(const Tps< T > &x)
Cosine.
gsl_interp_accel * onAxisFieldAccel_m
On axis field third derivative interpolation structure.
void getLine(std::ifstream &in, std::string &buffer)
double deltaZ_m
Number of grid points in field input file.
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
void normalizeField(double maxEz, std::vector< double > &fourierCoefs)
Inform & level3(Inform &inf)
gsl_spline * onAxisFieldPInterpolants_m
On axis field interpolation structure.
void prepareForMapCheck(std::vector< double > &fourierCoefs)
void checkMap(unsigned int accuracy, std::pair< double, double > fieldDimensions, double deltaZ, const std::vector< double > &fourierCoefficients, gsl_spline *splineCoefficients, gsl_interp_accel *splineAccelerator)
bool checkFileData(std::ifstream &fieldFile, bool parsingPassed)
virtual double getFrequency() const
double readFileData(std::ifstream &fieldFile, double fieldData[])
unsigned int accuracy_m
Field grid point spacing.
virtual void getInfo(Inform *)
Inform & endl(Inform &inf)
double zEnd_m
Longitudinal start of field.
gsl_interp_accel * onAxisFieldPPPAccel_m