4 #ifdef ENABLE_WAKE_TESTS
14 #include "gsl/gsl_fft_real.h"
15 #include "gsl/gsl_fft_halfcomplex.h"
48 vector<Filter *> filters,
67 direction_m(direction),
68 constLength_m(constLength),
70 filters_m(filters.begin(), filters.end()) {
71 #ifdef ENABLE_WAKE_DEBUG
72 *
gmsg <<
"* ************* W A K E ************************************************************ " <<
endl;
73 *
gmsg <<
"* Entered GreenWakeFunction::GreenWakeFunction " <<
'\n';
74 *
gmsg <<
"* ********************************************************************************** " <<
endl;
104 dist.first = locBunchRange *
Ippl::myNode() + (1 - tmp) * rem;
105 dist.second = dist.first + locBunchRange - 1;
111 #ifdef ENABLE_WAKE_TESTS
122 double spacing, mindist;
123 std::vector<double> OutEnergy(
NBin_m);
132 spacing =
abs(rmax(2) - rmin(2));
135 spacing = rmax(0) * rmax(0) + rmax(1) * rmax(1);
156 std::pair<double, double> meshInfo;
159 #ifdef ENABLE_WAKE_DEBUG
160 *
gmsg <<
"* ************* W A K E ************************************************************ " <<
endl;
162 *
gmsg <<
"* ********************************************************************************** " <<
endl;
166 for(vector<Filter *>::const_iterator fit =
filters_m.begin(); fit !=
filters_m.end(); ++fit) {
182 for(
unsigned int i = 0; i < bunch->
getLocalNum(); i++) {
186 int idx = (int)(
floor((bunch->
R[i](2) - mindist) / spacing));
189 assert(idx >= 0 && idx <
NBin_m);
190 double dE = OutEnergy[idx];
191 bunch->
Ef[i](2) += dE;
197 for(
unsigned int i = 0; i < bunch->
getLocalNum(); i++) {
200 int idx = (int)(
floor((bunch->
R[i](2) - mindist) / spacing));
203 assert(idx >= 0 && idx <
NBin_m);
204 double dE = OutEnergy[idx];
207 double dist =
sqrt(bunch->
R[i](0) * bunch->
R[i](0) + bunch->
R[i](1) * bunch->
R[i](1));
209 bunch->
Ef[i](0) += dE * bunch->
R[i](0) / dist;
210 bunch->
Ef[i](1) += dE * bunch->
R[i](1) / dist;
219 #ifdef ENABLE_WAKE_DUMP
220 ofstream f2(
"OutEnergy.dat");
221 f2 <<
"# Energy of the Wake calculated in Opal\n"
222 <<
"# Z0 = " <<
Z0_m <<
"\n"
223 <<
"# radius = " << radius <<
"\n"
224 <<
"# sigma = " << sigma <<
"\n"
225 <<
"# c = " <<
c <<
"\n"
226 <<
"# acMode = " << acMode <<
"\n"
227 <<
"# tau = " << tau <<
"\n"
228 <<
"# direction = " << direction <<
"\n"
229 <<
"# spacing = " << spacing <<
"\n"
230 <<
"# Lbunch = " <<
NBin_m <<
"\n";
231 for(
int i = 0; i <
NBin_m; i++) {
232 f2 << i + 1 <<
" " << OutEnergy[i] <<
"\n";
238 #endif //ENABLE_WAKE_TESTS
245 #ifdef ENABLE_WAKE_TESTS
248 double charge = 0.8e-9;
249 double K = 0.20536314319923724e-9;
252 std::vector<double> OutEnergy(
NBin_m);
263 ofstream f2(
"OutEnergy.dat");
264 f2 <<
"# Energy of the Wake calculated in Opal\n"
265 <<
"# Z0 = " <<
Z0_m <<
"\n"
266 <<
"# radius = " <<
radius_m <<
"\n"
267 <<
"# sigma = " <<
sigma_m <<
"\n"
268 <<
"# acMode = " <<
acMode_m <<
"\n"
269 <<
"# tau = " <<
tau_m <<
"\n"
271 <<
"# spacing = " << spacing_m <<
"\n"
272 <<
"# Lbunch = " <<
NBin_m <<
"\n";
273 for(
int i = 0; i <
NBin_m; i++) {
274 f2 << i + 1 <<
" " << OutEnergy[i] <<
"\n";
292 const double *lambda,
296 std::vector<double> pLambda(N);
298 gsl_fft_halfcomplex_wavetable *hc;
299 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(N);
300 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(N);
303 for(
int i = 0; i <
NBin_m; i ++) {
304 pLambda[N-i-1] = 0.0;
305 pLambda[i] = lambda[i];
309 gsl_fft_real_transform(pLambda.data(), 1, N,
real, work);
310 gsl_fft_real_wavetable_free(real);
314 for(
int i = 1; i < N; i += 2) {
315 double temp = pLambda[i];
321 hc = gsl_fft_halfcomplex_wavetable_alloc(N);
323 gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc, work);
326 for(
int i = 0; i <
NBin_m; i ++) {
327 OutEnergy[i] = -charge * K * pLambda[i] / (2.0 *
NBin_m) * N;
336 gsl_fft_halfcomplex_wavetable_free(hc);
337 gsl_fft_real_workspace_free(work);
352 vector<double> lambda,
356 std::vector<double> pLambda(N);
358 gsl_fft_halfcomplex_wavetable *hc;
359 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(N);
360 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(N);
363 for(
int i = 0; i <
NBin_m; i ++) {
364 pLambda[N-i-1] = 0.0;
365 pLambda[i] = lambda[i];
369 gsl_fft_real_transform(pLambda.data(), 1, N,
real, work);
370 gsl_fft_real_wavetable_free(real);
375 for(
int i = 1; i < N; i += 2) {
376 double temp = pLambda[i];
382 hc = gsl_fft_halfcomplex_wavetable_alloc(N);
384 gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc, work);
387 for(
int i = 0; i <
NBin_m; i ++) {
388 OutEnergy[i] = -charge * K * pLambda[i] / (2.0 *
NBin_m) * N;
397 gsl_fft_halfcomplex_wavetable_free(hc);
398 gsl_fft_real_workspace_free(work);
409 double a = 1, b = 1000000;
410 unsigned int N = 1000000;
413 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(M);
414 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(M);
417 const int lowIndex = myDist.first;
418 const int hiIndex = myDist.second;
420 #ifdef ENABLE_WAKE_TESTS
424 file.open(
"wake.dat");
425 file <<
"# Wake calculated in Opal" <<
"\n"
426 <<
"# Z0 = " <<
Z0_m <<
"\n"
427 <<
"# radius = " <<
radius_m <<
"\n"
428 <<
"# sigma = " <<
sigma_m <<
"\n"
430 <<
"# tau = " <<
tau_m <<
"\n"
432 <<
"# spacing = " << spacing <<
"\n"
433 <<
"# Lbunch = " <<
NBin_m <<
"\n";
437 for(
int i = 0; i < M; i ++) {
446 for(
int i = lowIndex; i <= hiIndex; i ++) {
464 #ifdef ENABLE_WAKE_TESTS
466 for(
int i = 0; i <
NBin_m; i++) {
474 #ifdef ENABLE_WAKE_TESTS_FFT_OUT
475 std::vector<double> wf(2*
NBin_m-1);
476 for(
int i = 0; i < 2 *
NBin_m - 1; ++ i) {
484 #ifdef ENABLE_WAKE_TESTS_FFT_OUT
485 ofstream f2(
"FFTwake.dat");
486 f2 <<
"# FFT of the Wake calculated in Opal" <<
"\n"
487 <<
"# Z0 = " <<
Z0_m <<
"\n"
488 <<
"# radius = " <<
radius_m <<
"\n"
489 <<
"# sigma = " <<
sigma_m <<
"\n"
491 <<
"# tau = " <<
tau_m <<
"\n"
493 <<
"# spacing = " << spacing <<
"\n"
494 <<
"# Lbunch = " << NBin_m <<
"\n";
496 f2 <<
"0\t" <<
FftWField_m[0] <<
"\t0.0\t" << wf[0] <<
"\n";
497 for(
int i = 1; i < M; i += 2) {
498 f2 << (i + 1) / 2 <<
"\t"
501 << wf[(i+1)/2] <<
"\n";
508 gsl_fft_real_wavetable_free(real);
509 gsl_fft_real_workspace_free(work);
521 gsl_fft_real_wavetable *
real;
522 gsl_fft_real_workspace *work;
529 "Open file operation failed, please check if \""
534 msg <<
" SSDS1 read = " << name <<
endl;
535 if(name.compare(
"SDDS1") != 0) {
537 " No SDDS1 File. A SDDS1 file should start with a SDDS1 String. Check file \""
541 for(
int i = 0; i < 6; i++) {
542 fs.getline(temp, 256);
543 msg <<
"line " << i <<
" : " << temp <<
endl;
547 msg <<
" header read" <<
endl;
550 " The particle number should be bigger than zero! Please check the first line of file \""
554 msg <<
" Np = " << Np <<
endl;
555 std::vector<double> wake(Np);
556 std::vector<double> dist(Np);
559 for(
int i = 0; i < Np; i ++) {
561 fs >> dist[i] >> wake[i] >> dummy;
565 " End of file reached before the whole wakefield is imported, please check file \""
573 for(
int i = 0; i <
NBin_m; i ++) {
575 while(dist[j] < i * spacing) {
580 FftWField_m[i] = wake[j] + ((wake[j+1] - wake[j]) / (dist[j+1] - dist[j]) * (i * spacing - dist[j]));
583 real = gsl_fft_real_wavetable_alloc(NBin_m);
584 work = gsl_fft_real_workspace_alloc(NBin_m);
588 gsl_fft_real_wavetable_free(real);
589 gsl_fft_real_workspace_free(work);
593 return "GreenWakeFunction";
std::vector< double > FftWField_m
FFT of the zero padded wakefield.
void setWakeFromFile(int NBin, double spacing)
reads in the wakefield from file
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
constexpr double e
The value of .
ParticleAttrib< Vector_t > Ef
const unsigned int nBins_m
Interface for basic beam line object.
GreenWakeFunction(const std::string &name, ElementBase *element, std::vector< Filter * > filters, int NBIN, double Z0, double radius, double sigma, int acMode, double tau, int direction, bool constLength, std::string fname)
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
double sigma_m
material constant
std::pair< int, int > distrIndices(int vectLen)
given a vector of length N, distribute the indexes among the available processors ...
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
virtual const std::string getType() const
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 reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
void calcBeamParameters()
const double testLambda[294]
int acMode_m
conductivity either 1="AC" or 2="DC"
constexpr double c
The velocity of light in m/s.
int direction_m
direction either 1="Longitudinal" 2= "Transversal"
void testApply(PartBunchBase< double, 3 > *bunch)
Just a test function.
double getChargePerParticle() const
get the macro particle charge
void apply(PartBunchBase< double, 3 > *bunch)
double tau_m
material constant
void CalcWakeFFT(double spacing)
Calculate the FFT of the Wakefunction.
size_t getLocalNum() const
std::vector< Filter * > filters_m
Tps< T > sqrt(const Tps< T > &x)
Square root.
int NBin_m
divides the particle bunch in NBin slices
bool constLength_m
true if the length of the particle bunch is considered as constant
std::string filename_m
filename of the wakefield
void get_bounds(Vector_t &rmin, Vector_t &rmax)
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
void calcLineDensity(unsigned int nBins, std::vector< double > &lineDensity, std::pair< double, double > &meshInfo)
calculates the 1d line density (not normalized) and append it to a file.
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Inform & endl(Inform &inf)