20 filters_m(filters.begin(), filters.end()),
37 const double sPos = bunch->
get_sPos();
38 std::pair<double, double> meshInfo;
40 const double &meshSpacing = meshInfo.second;
41 const double &meshOrigin = meshInfo.first + 0.5 * meshSpacing;
51 if (sPos + smax(2) < FieldBegin_m)
return;
58 double angleOfSlice = 0.0;
59 double pathLengthOfSlice = minPathLength + i * meshSpacing;
60 if(pathLengthOfSlice > 0.0)
70 for(
unsigned int i = 0; i < bunch->
getLocalNum(); ++i) {
72 double distanceToOrigin = (
R(2) - meshOrigin) / meshSpacing;
74 unsigned int indexz = (
unsigned int)
floor(distanceToOrigin);
75 double leverz = distanceToOrigin - indexz;
78 bunch->
Ef[i](2) += (1. - leverz) *
Ez_m[indexz] + leverz *
Ez_m[indexz + 1];
82 static std::string oldBendName;
83 static unsigned long counter = 0;
88 bool print_criterion = (counter + 1) % every == 0;
90 static unsigned int file_number = 0;
91 if(counter == 0) file_number = 0;
93 std::stringstream filename_str;
94 filename_str <<
"data/" <<
bendName_m <<
"-CSRWake" << std::setw(5) << std::setfill(
'0') << file_number <<
".txt";
96 std::ofstream csr(filename_str.str().c_str());
97 csr << std::setprecision(8);
98 csr <<
"# " << sPos + smin(2) - FieldBegin_m <<
"\t" << sPos + smax(2) - FieldBegin_m <<
std::endl;
100 csr << i *meshSpacing <<
"\t"
106 msg <<
"** wrote " << filename_str.str() <<
endl;
134 std::vector<Filter *>::const_iterator fit;
146 const double meshSpacingsup =
pow(meshSpacing, -1. / 3.);
148 double relativeSlippageLength = SlippageLength / meshSpacing;
149 if(relativeSlippageLength > sliceNumber) {
156 double dx1 =
pow(sliceNumber, 2. / 3.);
157 double dx2 =
pow(sliceNumber, 5. / 3.);
158 double dx3 =
pow(sliceNumber - 1., 5. / 3.);
159 Ez_m[sliceNumber] += 0.3 * meshSpacingsup *
dlineDensitydz_m[0] * (5. * dx1 - 3. * dx2 + 3. * dx3);
160 for(
unsigned int j = 1; j < sliceNumber; ++ j) {
163 dx3 =
pow(sliceNumber - j - 1., 5. / 3.);
168 }
else if(relativeSlippageLength < 1) {
171 if(4.0 * relativeSlippageLength <= 1) {
176 if(4.0 * relativeSlippageLength < sliceNumber) {
178 int j = sliceNumber -
static_cast<int>(
floor(4.0 * relativeSlippageLength));
179 double frac = 4.0 * relativeSlippageLength - (sliceNumber - j);
184 Ez_m[sliceNumber] += (relativeSlippageLength *
lineDensity_m[sliceNumber - 1] + (1. - relativeSlippageLength) *
lineDensity_m[sliceNumber]) /
pow(SlippageLength, 1. / 3.);
193 if(4. * relativeSlippageLength < sliceNumber) {
195 int j = sliceNumber -
static_cast<int>(
floor(4. * relativeSlippageLength));
196 double frac = 4. * relativeSlippageLength - (sliceNumber - j);
201 int j = sliceNumber -
static_cast<int>(
floor(SlippageLength / meshSpacing));
202 double frac = relativeSlippageLength - (sliceNumber - j);
205 double dx1 =
pow(sliceNumber - j + frac, 2. / 3.);
206 double dx2 =
pow(sliceNumber - j, 2. / 3.);
207 double dx3 =
pow(sliceNumber - j + frac, 5. / 3.);
208 double dx4 =
pow(sliceNumber - j, 5. / 3.);
210 Ez_m[sliceNumber] += 1.5 * meshSpacingsup *
dlineDensitydz_m[j - 1] * (dx1 - dx2);
215 dx3 =
pow(sliceNumber - j - 1., 5. / 3.);
216 Ez_m[sliceNumber] += 0.3 * meshSpacingsup *
dlineDensitydz_m[j] * (5.*dx1 - 3.*dx2 + 3.*dx3);
217 for(
unsigned int k = j + 1; k < sliceNumber; ++ k) {
220 dx3 =
pow(sliceNumber - k - 1., 5. / 3.);
221 Ez_m[sliceNumber] += 0.9 * meshSpacingsup *
dlineDensitydz_m[k] * (dx1 - 2.*dx2 + dx3);
226 Ez_m[sliceNumber] *= prefactor;
239 if(Ds_max2 / meshSpacing < sliceNumber) {
240 j = sliceNumber -
static_cast<int>(
floor(Ds_max2 / meshSpacing));
241 frac = Ds_max2 / meshSpacing - (sliceNumber - j);
247 if(Ds_max / meshSpacing < sliceNumber) {
248 j = sliceNumber -
static_cast<int>(
floor(Ds_max / meshSpacing));
249 frac = Ds_max / meshSpacing - (sliceNumber - j);
258 if(Ds_max / meshSpacing < sliceNumber) {
261 if(
Psi_m[k] > 0 &&
Psi_m[k] < totalBendAngle_m)
265 if(
Psi_m[0] > 0 &&
Psi_m[0] < totalBendAngle_m)
270 for(
unsigned int l = sliceNumber - k + 1; l < sliceNumber; ++ l) {
272 if(
Psi_m[l] > 0 &&
Psi_m[l] < totalBendAngle_m)
278 Psi_m[sliceNumber] =
calcPsi(
Psi_m[sliceNumber], angleOverlap, meshSpacing / 4.0);
279 if(
Psi_m[sliceNumber] > 0 &&
Psi_m[sliceNumber] < totalBendAngle_m)
283 Ez_m[sliceNumber] *= prefactor;
294 const int Nmax = 100;
295 const double eps = 1
e-10;
297 double residual = 0.0;
299 if(psiInitial != 0.0) psi = psiInitial;
301 for(
int i = 0; i < Nmax; ++i) {
302 residual =
bendRadius_m * psi * psi * psi * (psi + 4. * x) - 24. * Ds * psi - 24. * Ds * x;
314 ERRORMSG(
"In CSRWakeFunction::calcPsi(): exceed maximum number of iterations!" <<
endl);
319 return "CSRWakeFunction";
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
void calculateContributionInside(size_t sliceNumber, double angleOfSlice, double meshSpacing)
constexpr double e
The value of .
virtual void getDimensions(double &sBegin, double &sEnd) const
ParticleAttrib< Vector_t > Ef
const unsigned int nBins_m
Interface for basic beam line object.
CSRWakeFunction(const std::string &name, ElementBase *element, std::vector< Filter * > filters, const unsigned int &N)
double getBendRadius() const
virtual const std::string & getName() const
Get element name.
std::shared_ptr< Filter > defaultFilter_m
void initialize(const ElementBase *ref)
virtual ElementType getType() const =0
Get element type std::string.
double getEffectiveLength() const
std::vector< Filter * > filters_m
double getBendAngle() const
std::vector< double > Psi_m
constexpr double pi
The value of .
void calculateContributionAfter(size_t sliceNumber, double angleOfSlice, double meshSpacing)
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
std::vector< double > Ez_m
double calcPsi(const double &psiInitial, const double &x, const double &Ds) const
size_t getLocalNum() const
LineDensity dlineDensitydz_m
double searchRoot(const double &tol)
void calculateLineDensity(PartBunchBase< double, 3 > *bunch, std::pair< double, double > &meshInfo)
bool hasPositiveRealRoots()
constexpr double epsilon_0
The permittivity of vacuum in As/Vm.
double getEffectiveCenter() const
void get_bounds(Vector_t &rmin, Vector_t &rmax)
LineDensity lineDensity_m
virtual const std::string getType() const
virtual void calc_derivative(std::vector< double > &histogram, const double &h)=0
void apply(PartBunchBase< double, 3 > *bunch)
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)