22 #include "gsl/gsl_fft_real.h"
23 #include "gsl/gsl_fft_halfcomplex.h"
57 std::vector<Filter *> filters,
76 direction_m(direction),
77 constLength_m(constLength),
79 filters_m(filters.
begin(), filters.
end()) {
80 #ifdef ENABLE_WAKE_DEBUG
81 *
gmsg <<
"* ************* W A K E ************************************************************ " <<
endl;
82 *
gmsg <<
"* Entered GreenWakeFunction::GreenWakeFunction " <<
'\n';
83 *
gmsg <<
"* ********************************************************************************** " <<
endl;
106 std::pair<int, int> dist;
113 dist.first = locBunchRange *
Ippl::myNode() + (1 - tmp) * rem;
114 dist.second = dist.first + locBunchRange - 1;
127 double spacing, mindist;
128 std::vector<double> OutEnergy(
NBin_m);
137 spacing =
std::abs(rmax(2) - rmin(2));
141 spacing = rmax(0) * rmax(0) + rmax(1) * rmax(1);
146 "Invalid direction specified");
165 std::pair<double, double> meshInfo;
168 #ifdef ENABLE_WAKE_DEBUG
169 *
gmsg <<
"* ************* W A K E ************************************************************ " <<
endl;
171 *
gmsg <<
"* ********************************************************************************** " <<
endl;
175 for (std::vector<Filter *>::const_iterator fit =
filters_m.begin(); fit !=
filters_m.end(); ++fit) {
191 for (
unsigned int i = 0; i < bunch->
getLocalNum(); i++) {
195 int idx = (int)(
floor((bunch->
R[i](2) - mindist) / spacing));
199 double dE = OutEnergy[idx];
200 bunch->
Ef[i](2) += dE;
206 for (
unsigned int i = 0; i < bunch->
getLocalNum(); i++) {
209 int idx = (int)(
floor((bunch->
R[i](2) - mindist) / spacing));
213 double dE = OutEnergy[idx];
216 double dist =
sqrt(bunch->
R[i](0) * bunch->
R[i](0) + bunch->
R[i](1) * bunch->
R[i](1));
218 bunch->
Ef[i](0) += dE * bunch->
R[i](0) / dist;
219 bunch->
Ef[i](1) += dE * bunch->
R[i](1) / dist;
226 "Invalid direction specified");
230 #ifdef ENABLE_WAKE_DUMP
231 ofstream f2(
"OutEnergy.dat");
232 f2 <<
"# Energy of the Wake calculated in Opal\n"
233 <<
"# Z0 = " <<
Z0_m <<
"\n"
234 <<
"# radius = " << radius <<
"\n"
235 <<
"# sigma = " << sigma <<
"\n"
236 <<
"# c = " <<
c <<
"\n"
237 <<
"# acMode = " << acMode <<
"\n"
238 <<
"# tau = " << tau <<
"\n"
239 <<
"# direction = " << direction <<
"\n"
240 <<
"# spacing = " << spacing <<
"\n"
241 <<
"# Lbunch = " <<
NBin_m <<
"\n";
242 for (
int i = 0; i <
NBin_m; i++) {
243 f2 << i + 1 <<
" " << OutEnergy[i] <<
"\n";
261 const double* lambda,
265 std::vector<double> pLambda(N);
267 gsl_fft_halfcomplex_wavetable *hc;
268 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(N);
269 gsl_fft_real_workspace *
work = gsl_fft_real_workspace_alloc(N);
272 for (
int i = 0; i <
NBin_m; i ++) {
273 pLambda[N-i-1] = 0.0;
274 pLambda[i] = lambda[i];
278 gsl_fft_real_transform(pLambda.data(), 1, N,
real,
work);
279 gsl_fft_real_wavetable_free(real);
283 for (
int i = 1; i < N; i += 2) {
284 double temp = pLambda[i];
290 hc = gsl_fft_halfcomplex_wavetable_alloc(N);
292 gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc,
work);
295 for (
int i = 0; i <
NBin_m; i ++) {
296 OutEnergy[i] = -charge * K * pLambda[i] / (2.0 *
NBin_m) * N;
305 gsl_fft_halfcomplex_wavetable_free(hc);
306 gsl_fft_real_workspace_free(work);
321 std::vector<double> lambda,
325 std::vector<double> pLambda(N);
327 gsl_fft_halfcomplex_wavetable *hc;
328 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(N);
329 gsl_fft_real_workspace *
work = gsl_fft_real_workspace_alloc(N);
332 for (
int i = 0; i <
NBin_m; i ++) {
333 pLambda[N-i-1] = 0.0;
334 pLambda[i] = lambda[i];
338 gsl_fft_real_transform(pLambda.data(), 1, N,
real,
work);
339 gsl_fft_real_wavetable_free(real);
344 for (
int i = 1; i < N; i += 2) {
345 double temp = pLambda[i];
351 hc = gsl_fft_halfcomplex_wavetable_alloc(N);
353 gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc,
work);
356 for (
int i = 0; i <
NBin_m; i ++) {
357 OutEnergy[i] = -charge * K * pLambda[i] / (2.0 *
NBin_m) * N;
366 gsl_fft_halfcomplex_wavetable_free(hc);
367 gsl_fft_real_workspace_free(work);
378 double a = 1, b = 1000000;
379 unsigned int N = 1000000;
382 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(M);
383 gsl_fft_real_workspace *
work = gsl_fft_real_workspace_alloc(M);
386 const int lowIndex = myDist.first;
387 const int hiIndex = myDist.second;
389 for (
int i = 0; i < M; i ++) {
398 for (
int i = lowIndex; i <= hiIndex; i ++) {
415 #ifdef ENABLE_WAKE_TESTS_FFT_OUT
416 std::vector<double> wf(2*
NBin_m-1);
417 for (
int i = 0; i < 2 *
NBin_m - 1; ++ i) {
425 #ifdef ENABLE_WAKE_TESTS_FFT_OUT
426 ofstream f2(
"FFTwake.dat");
427 f2 <<
"# FFT of the Wake calculated in Opal" <<
"\n"
428 <<
"# Z0 = " <<
Z0_m <<
"\n"
429 <<
"# radius = " <<
radius_m <<
"\n"
430 <<
"# sigma = " <<
sigma_m <<
"\n"
432 <<
"# tau = " <<
tau_m <<
"\n"
434 <<
"# spacing = " << spacing <<
"\n"
435 <<
"# Lbunch = " << NBin_m <<
"\n";
437 f2 <<
"0\t" <<
FftWField_m[0] <<
"\t0.0\t" << wf[0] <<
"\n";
438 for (
int i = 1; i < M; i += 2) {
439 f2 << (i + 1) / 2 <<
"\t"
442 << wf[(i+1)/2] <<
"\n";
449 gsl_fft_real_wavetable_free(real);
450 gsl_fft_real_workspace_free(work);
461 gsl_fft_real_wavetable *
real;
462 gsl_fft_real_workspace *
work;
469 "Open file operation failed, please check if '"
474 msg <<
" SSDS1 read = " << name <<
endl;
475 if (name.compare(
"SDDS1") != 0) {
477 " No SDDS1 File. A SDDS1 file should start with a SDDS1 String. Check file '"
481 for (
int i = 0; i < 6; i++) {
482 fs.getline(temp, 256);
483 msg <<
"line " << i <<
" : " << temp <<
endl;
487 msg <<
" header read" <<
endl;
490 " The particle number should be bigger than zero! Please check the first line of file '"
494 msg <<
" Np = " << Np <<
endl;
495 std::vector<double> wake(Np);
496 std::vector<double> dist(Np);
499 for (
int i = 0; i < Np; i ++) {
501 fs >> dist[i] >> wake[i];
505 " End of file reached before the whole wakefield is imported, please check file '"
513 for (
int i = 0; i <
NBin_m; i ++) {
515 while(dist[j] < i * spacing) {
520 FftWField_m[i] = wake[j] + ((wake[j+1] - wake[j]) / (dist[j+1] - dist[j]) * (i * spacing - dist[j]));
523 real = gsl_fft_real_wavetable_alloc(NBin_m);
524 work = gsl_fft_real_workspace_alloc(NBin_m);
528 gsl_fft_real_wavetable_free(real);
529 gsl_fft_real_workspace_free(work);
const unsigned int nBins_m
void setWakeFromFile(int NBin, double spacing)
reads in the wakefield from file
Tps< T > sqrt(const Tps< T > &x)
Square root.
constexpr double c
The velocity of light in m/s.
std::string filename_m
filename of the wakefield
bool constLength_m
true if the length of the particle bunch is considered as constant
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
double simpson(F &f, double a, double b, unsigned int N)
Simpson-Integration from the function f from a to b with N steps.
static std::string getWakeDirectionString(const WakeDirection &direction)
ParticleAttrib< Vector_t > Ef
void apply(PartBunchBase< double, 3 > *bunch) override
void get_bounds(Vector_t &rmin, Vector_t &rmax) const
void calcBeamParameters()
GreenWakeFunction(const std::string &name, std::vector< Filter * > filters, int NBIN, double Z0, double radius, double sigma, int acMode, double tau, WakeDirection direction, bool constLength, std::string fname)
clearpage the user may choose between constant or variable radius This model includes fringe fields begin
Inform & endl(Inform &inf)
size_t getLocalNum() const
int NBin_m
divides the particle bunch in NBin slices
int acMode_m
conductivity either 1="AC" or 2="DC"
std::vector< double > lineDensity_m
save the line Density of the particle bunch
static const std::map< WakeDirection, std::string > wakeDirectiontoString_s
double getChargePerParticle() const
get the macro particle charge
std::vector< double > FftWField_m
FFT of the zero padded wakefield.
std::vector< Filter * > filters_m
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
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.
double sigma_m
material constant
double tau_m
material constant
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
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
WakeDirection direction_m
direction either "Longitudinal" - "Transversal"
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
void CalcWakeFFT(double spacing)
Calculate the FFT of the Wakefunction.
virtual WakeType getType() const override
std::pair< int, int > distrIndices(int vectLen)
given a vector of length N, distribute the indexes among the available processors ...
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)