OPAL (Object Oriented Parallel Accelerator Library)  2021.1.99
OPAL
GreenWakeFunction.h
Go to the documentation of this file.
1 //
2 // Class GreenWakeFunction
3 //
4 // Copyright (c) 2008 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
5 // All rights reserved
6 //
7 // This file is part of OPAL.
8 //
9 // OPAL is free software: you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation, either version 3 of the License, or
12 // (at your option) any later version.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with OPAL. If not, see <https://www.gnu.org/licenses/>.
16 //
17 #ifndef GREENWAKEFUNCTION_HH
18 #define GREENWAKEFUNCTION_HH
19 
20 #include "Filters/Filter.h"
21 #include "Solvers/WakeFunction.h"
22 #include "Physics/Physics.h"
23 #include "Utility/IpplInfo.h"
24 #include "Utility/PAssert.h"
25 
26 #include <vector>
27 #include <map>
28 #include <string>
29 #include <complex>
30 
31 #ifdef WITH_UNIT_TESTS
32 #include <gtest/gtest_prod.h>
33 #endif
34 
36 typedef std::map<std::string, int> FilterOptions;
37 
39 public:
41  //IFF: changed direction to int (was double)
42  //IFF: changed acMode to int (was double)
43  GreenWakeFunction(const std::string &name,
44  std::vector<Filter *> filters,
45  int NBIN,
46  double Z0,
47  double radius,
48  double sigma,
49  int acMode,
50  double tau,
51  int direction,
52  bool constLength,
53  std::string fname);
54 
55  std::pair<int, int> distrIndices(int vectLen);
56 
57  void apply(PartBunchBase<double, 3> *bunch);
58  void setWakeFromFile(int NBin, double spacing);
59  virtual const std::string getType() const;
60 
61 private:
62 #ifdef WITH_UNIT_TESTS
63  FRIEND_TEST(GreenWakeFunctionTest, TestApply);
64 #endif
65 
66  class Wake {
67 
68  public:
69 
70  Wake(double s, double Z0, double a, double sigma, int acMode, double tau, int direction)
71  : Z0_(Z0), a_(a), sigma_(sigma), s_(s), acMode_(acMode), tau_(tau), direction_(direction)
72  {}
73 
81  double operator()(double k) {
82 
83  std::complex <double> i(0, 1);
84  std::complex <double> Z(0, 0);
85  double signK = (k > 0 ? 1 : -1);
86 
87  //1 == AC
88  //2 == DC
89  switch(acMode_) {
90  case 1:
91  Z = (Z0_ / (2 * Physics::pi * a_)) * 1.0 / (sqrt(Z0_ * std::abs(k) / 2) * sqrt(sigma_ / (1.0 - i * Physics::c * k * tau_)) * (i + signK) / k - (i * k * a_) / 2.0);
92  break;
93  case 2:
94  Z = (Z0_ / (2 * Physics::pi * a_)) * 1.0 / (sqrt(sigma_ * Z0_ * std::abs(k) / 2) * (i + signK) / k - (i * k * a_) / 2.0);
95  break;
96  }
97  switch(direction_) {
98  case LONGITUDINAL:
99  return real(Z) * cos(k * s_) * 2.0 * Physics::c / Physics::pi;
100  break;
101  case TRANSVERSAL:
102  return real(Z) * Physics::c / k * cos(k * s_) * 2.0 * Physics::c / Physics::pi;
103  break;
104  }
105  ERRORMSG("We should not be here: " << __FILE__ << " L" << __LINE__ << endl);
106 
107  return 0.0;
108  }
109 
110  private:
111 
113  double Z0_;
115  double a_;
117  double sigma_;
119  double s_;
121  int acMode_;
123  double tau_;
126 
127  };
128 
140  template<class F> double simpson(F &f, double a, double b, unsigned int N) {
141  PAssert(b > a);
142  PAssert(N > 0);
143 
144  double result = 0;
145  double h = (b - a) / N;
146 
147  // boundary values
148  result += (f(a) + 4 * f(a + h / 2) + f(b)) / 2.0;
149 
150  // values between boundaries
151  for(unsigned int i = 1; i < N; ++ i) {
152  result += f(a + i * h) + 2 * f(a + (i + 0.5) * h);
153  }
154 
155  result *= h / 3.0;
156 
157  return result;
158 
159  }
161  std::vector<double> lineDensity_m;
163  std::vector<double> FftWField_m;
164 
166  int NBin_m;
168  double Z0_m;
170  double radius_m;
172  double sigma_m;
174  int acMode_m;
176  double tau_m;
182  std::string filename_m;
183 
184  std::vector<Filter *> filters_m;
185 
186  void compEnergy(const double K, const double charge, const double *lambda, double *OutEnergy);
187  void compEnergy(const double K, const double charge, std::vector<double> lambda, double *OutEnergy);
188  void CalcWakeFFT(double spacing);
189 };
190 #endif //GREENWAKEFUNCTION_HH
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition: TpsMath.h:129
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
std::map< std::string, int > FilterOptions
@ LONGITUDINAL
@ TRANSVERSAL
std::complex< double > a
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
#define ERRORMSG(msg)
Definition: IpplInfo.h:350
#define PAssert(c)
Definition: PAssert.h:102
const std::string name
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:51
constexpr double pi
The value of.
Definition: Physics.h:30
std::vector< double > FftWField_m
FFT of the zero padded wakefield.
void setWakeFromFile(int NBin, double spacing)
reads in the wakefield from file
int direction_m
direction either 1="Longitudinal" 2= "Transversal"
int NBin_m
divides the particle bunch in NBin slices
std::string filename_m
filename of the wakefield
void compEnergy(const double K, const double charge, const double *lambda, double *OutEnergy)
just a Testfunction! Calculate the energy of the Wakefunction with the lambda
std::vector< double > lineDensity_m
save the line Density of the particle bunch
GreenWakeFunction(const std::string &name, std::vector< Filter * > filters, int NBIN, double Z0, double radius, double sigma, int acMode, double tau, int direction, bool constLength, std::string fname)
double simpson(F &f, double a, double b, unsigned int N)
Simpson-Integration from the function f from a to b with N steps.
bool constLength_m
true if the length of the particle bunch is considered as constant
void apply(PartBunchBase< double, 3 > *bunch)
int acMode_m
conductivity either 1="AC" or 2="DC"
double sigma_m
material constant
void CalcWakeFFT(double spacing)
Calculate the FFT of the Wakefunction.
virtual const std::string getType() const
double Z0_m
impedance
std::pair< int, int > distrIndices(int vectLen)
given a vector of length N, distribute the indexes among the available processors
double tau_m
material constant
std::vector< Filter * > filters_m
double tau_
material constant
double operator()(double k)
Used to integrate the function.
int direction_
direction either 1="Longitudinal" 0= "Transversal"
int acMode_
conductivity either 1="AC" or 2="DC"
Wake(double s, double Z0, double a, double sigma, int acMode, double tau, int direction)
double s_
distance from the particle
double sigma_
material constant