Skip to content

Commit 56500b9

Browse files
azukovazukov
committed
Organizing extensions in /ext.
The namespace will be orbit.ext
1 parent 285f689 commit 56500b9

16 files changed

+740
-14
lines changed

‎ext/danilov_envelope/DanilovEnvelopeSolver20.cc‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
#include "DanilovEnvelopeSolver20.hh"
2+
3+
DanilovEnvelopeSolver20::DanilovEnvelopeSolver20(double perveance, double eps_x, double eps_y) : CppPyWrapper(NULL) {
4+
_perveance = perveance;
5+
_eps_x = eps_x;
6+
_eps_y = eps_y;
7+
}
8+
9+
void DanilovEnvelopeSolver20::setPerveance(double perveance) {
10+
_perveance = perveance;
11+
}
12+
13+
void DanilovEnvelopeSolver20::setEmittanceX(double eps_x) {
14+
_eps_x = eps_x;
15+
}
16+
17+
void DanilovEnvelopeSolver20::setEmittanceY(double eps_y) {
18+
_eps_y = eps_y;
19+
}
20+
21+
double DanilovEnvelopeSolver20::getPerveance() {
22+
return _perveance;
23+
}
24+
25+
double DanilovEnvelopeSolver20::getEmittanceX() {
26+
return _eps_x;
27+
}
28+
29+
double DanilovEnvelopeSolver20::getEmittanceY() {
30+
return _eps_y;
31+
}
32+
33+
void DanilovEnvelopeSolver20::trackBunch(Bunch *bunch, double length) {
34+
// Track envelope
35+
double cx = bunch->x(0);
36+
double cy = bunch->y(0);
37+
double sc_term = 2.0 * _perveance / (cx + cy);
38+
double emit_term_x = (_eps_x * _eps_x) / (cx * cx * cx);
39+
double emit_term_y = (_eps_y * _eps_y) / (cy * cy * cy);
40+
bunch->xp(0) += length * (sc_term + emit_term_x);
41+
bunch->yp(0) += length * (sc_term + emit_term_y);
42+
43+
// Track test particles
44+
double cx2 = cx * cx;
45+
double cy2 = cy * cy;
46+
47+
double x;
48+
double y;
49+
double x2;
50+
double y2;
51+
52+
double B;
53+
double C;
54+
double Dx;
55+
double Dy;
56+
double t1;
57+
58+
double delta_xp;
59+
double delta_yp;
60+
bool in_ellipse;
61+
62+
for (int i = 1; i < bunch->getSize(); i++) {
63+
x = bunch->x(i);
64+
y = bunch->y(i);
65+
66+
x2 = x * x;
67+
y2 = y * y;
68+
69+
in_ellipse = ((x2 / cx2) + (y2 / cy2)) <= 1.0;
70+
71+
if (in_ellipse) {
72+
delta_xp = sc_term * x / cx;
73+
delta_yp = sc_term * y / cy;
74+
} else {
75+
// https://arxiv.org/abs/physics/0108040
76+
B = x2 + y2 - cx2 - cy2;
77+
C = x2 * cy2 + y2 * cx2 - cx2 * cy2;
78+
t1 = pow(0.25 * B * B + C, 0.5) + 0.5 * B;
79+
Dx = pow(cx2 + t1, 0.5);
80+
Dy = pow(cy2 + t1, 0.5);
81+
delta_xp = 2.0 * _perveance * x / (Dx * (Dx + Dy));
82+
delta_yp = 2.0 * _perveance * y / (Dy * (Dx + Dy));
83+
}
84+
bunch->xp(i) += delta_xp;
85+
bunch->yp(i) += delta_yp;
86+
}
87+
}

‎ext/danilov_envelope/DanilovEnvelopeSolver20.hh‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
#ifndef DANILOV_ENVELOPE_SOLVER_20_H
2+
#define DANILOV_ENVELOPE_SOLVER_20_H
3+
4+
#include "Bunch.hh"
5+
#include "CppPyWrapper.hh"
6+
7+
using namespace std;
8+
9+
/** Envelope solver for the {2, 0} Danilov distribution (KV distribution).
10+
11+
The ellipse in the x-y plane can be parameterized as:
12+
x = cx * cos(psi),
13+
y = cy * sin(psi),
14+
where 0 <= psi <= 2pi. The first particle in the bunch is used to track the
15+
envelope parameters {a, b}, which are used to apply space charge kicks to the
16+
other particles in the bunch.
17+
*/
18+
class DanilovEnvelopeSolver20 : public OrbitUtils::CppPyWrapper {
19+
public:
20+
DanilovEnvelopeSolver20(double perveance, double eps_x, double eps_y);
21+
void trackBunch(Bunch *bunch, double length);
22+
void setPerveance(double perveance);
23+
void setEmittanceX(double eps_x);
24+
void setEmittanceY(double eps_y);
25+
double getPerveance();
26+
double getEmittanceX();
27+
double getEmittanceY();
28+
29+
private:
30+
double _perveance; // beam perveance
31+
double _eps_x; // (4 * sqrt(<xx><x'x'> - <xx'><xx'>))
32+
double _eps_y; // (4 * sqrt(<yy><y'y'> - <yy'><yy'>))
33+
};
34+
35+
#endif

‎ext/danilov_envelope/DanilovEnvelopeSolver22.cc‎

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
#include "DanilovEnvelopeSolver22.hh"
2+
3+
DanilovEnvelopeSolver22::DanilovEnvelopeSolver22(double perveance) : CppPyWrapper(NULL) {
4+
Q = perveance;
5+
}
6+
7+
void DanilovEnvelopeSolver22::setPerveance(double perveance) {
8+
Q = perveance;
9+
}
10+
11+
double DanilovEnvelopeSolver22::getPerveance() {
12+
return Q;
13+
}
14+
15+
void DanilovEnvelopeSolver22::trackBunch(Bunch *bunch, double length) {
16+
// Compute ellipse size and orientation.
17+
double a = bunch->x(0);
18+
double b = bunch->x(1);
19+
double e = bunch->y(0);
20+
double f = bunch->y(1);
21+
22+
double cov_xx = a * a + b * b; // 4 * <x^2>
23+
double cov_yy = e * e + f * f; // 4 * <y^2>
24+
double cov_xy = a * e + b * f; // 4 * <xy>
25+
26+
double phi = -0.5 * atan2(2.0 * cov_xy, cov_xx - cov_yy);
27+
28+
double _cos = cos(phi);
29+
double _sin = sin(phi);
30+
double cos2 = _cos * _cos;
31+
double sin2 = _sin * _sin;
32+
double sin_cos = _sin * _cos;
33+
34+
double cxn = sqrt(abs(cov_xx * cos2 + cov_yy * sin2 - 2.0 * cov_xy * sin_cos));
35+
double cyn = sqrt(abs(cov_xx * sin2 + cov_yy * cos2 + 2.0 * cov_xy * sin_cos));
36+
double factor = length * (2.0 * Q / (cxn + cyn));
37+
38+
// Track envelope
39+
if (cxn > 0.0) {
40+
bunch->xp(0) += factor * (a * cos2 - e * sin_cos) / cxn;
41+
bunch->xp(1) += factor * (b * cos2 - f * sin_cos) / cxn;
42+
bunch->yp(0) += factor * (e * sin2 - a * sin_cos) / cxn;
43+
bunch->yp(1) += factor * (f * sin2 - b * sin_cos) / cxn;
44+
}
45+
if (cyn > 0.0) {
46+
bunch->xp(0) += factor * (a * sin2 + e * sin_cos) / cyn;
47+
bunch->xp(1) += factor * (b * sin2 + f * sin_cos) / cyn;
48+
bunch->yp(0) += factor * (e * cos2 + a * sin_cos) / cyn;
49+
bunch->yp(1) += factor * (f * cos2 + b * sin_cos) / cyn;
50+
}
51+
52+
// Track test particles
53+
double cxn2 = cxn * cxn;
54+
double cyn2 = cyn * cyn;
55+
double x;
56+
double y;
57+
double x2;
58+
double y2;
59+
60+
double xn;
61+
double yn;
62+
double xn2;
63+
double yn2;
64+
65+
double t1;
66+
double B;
67+
double C;
68+
double Dx;
69+
double Dy;
70+
71+
double delta_xp;
72+
double delta_yp;
73+
double delta_xpn;
74+
double delta_ypn;
75+
bool in_ellipse;
76+
77+
for (int i = 2; i < bunch->getSize(); i++) {
78+
x = bunch->x(i);
79+
y = bunch->y(i);
80+
81+
x2 = x * x;
82+
y2 = y * y;
83+
84+
xn = x * _cos - y * _sin;
85+
yn = x * _sin + y * _cos;
86+
87+
xn2 = xn * xn;
88+
yn2 = yn * yn;
89+
90+
in_ellipse = ((xn2 / cxn2) + (yn2 / cyn2)) <= 1.0;
91+
92+
if (in_ellipse) {
93+
if (cxn > 0.0) {
94+
delta_xpn = factor * xn / cxn;
95+
}
96+
if (cyn > 0.0) {
97+
delta_ypn = factor * yn / cyn;
98+
}
99+
} else {
100+
// https://arxiv.org/abs/physics/0108040
101+
B = xn2 + yn2 - cxn2 - cyn2;
102+
C = xn2 * cyn2 + yn2 * cxn2 - cxn2 * cyn2;
103+
t1 = pow(0.25 * B * B + C, 0.5) + 0.5 * B;
104+
Dx = pow(cxn2 + t1, 0.5);
105+
Dy = pow(cyn2 + t1, 0.5);
106+
delta_xpn = 2.0 * Q * xn / (Dx * (Dx + Dy));
107+
delta_ypn = 2.0 * Q * yn / (Dy * (Dx + Dy));
108+
}
109+
delta_xp = +delta_xpn * _cos + delta_ypn * _sin;
110+
delta_yp = -delta_xpn * _sin + delta_ypn * _cos;
111+
bunch->xp(i) += delta_xp;
112+
bunch->yp(i) += delta_yp;
113+
}
114+
}

‎ext/danilov_envelope/DanilovEnvelopeSolver22.hh‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
#ifndef DANILOV_ENVELOPE_SOLVER_22_H
2+
#define DANILOV_ENVELOPE_SOLVER_22_H
3+
4+
#include "Bunch.hh"
5+
#include "CppPyWrapper.hh"
6+
7+
using namespace std;
8+
9+
/** Envelope solver for the {2, 2} Danilov distribution (KV distribution with
10+
zero emittance in one plane).
11+
12+
The ellipse in the x-y plane can be parameterized as:
13+
x = a * cos(psi) + b * sin(psi),
14+
y = e * cos(psi) + f * sin(psi),
15+
where 0 <= psi <= 2pi. The first two particles in the bunch are used to track
16+
the envelope parameters {a, b, e, f}, which are used to apply space charge
17+
kicks to the other particles in the bunch.
18+
19+
Reference:
20+
V. Danilov, S. Cousineau, S. Henderson, and J. Holmes, "Self-consistent time
21+
dependent two dimensional and three dimensional space charge distributions with
22+
linear force", PPRAB 6, 74–85 (2003).
23+
*/
24+
class DanilovEnvelopeSolver22 : public OrbitUtils::CppPyWrapper {
25+
public:
26+
DanilovEnvelopeSolver22(double perveanceQ);
27+
void trackBunch(Bunch *bunch, double length);
28+
void setPerveance(double perveance);
29+
double getPerveance();
30+
31+
private:
32+
double Q; // beam perveance
33+
};
34+
35+
#endif

‎ext/danilov_envelope/danilov_envelope.cc‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
#include <Python.h>
2+
#include "wrap_orbit_mpi.hh"
3+
#include "wrap_danilov_envelope.hh"
4+
5+
PyMODINIT_FUNC PyInit_danilov_envelope(void) {
6+
return wrap_danilov_envelope::initdanilovenvelope();
7+
}

‎ext/danilov_envelope/meson.build‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
ext_sources = files([
2+
'DanilovEnvelopeSolver20.cc',
3+
'DanilovEnvelopeSolver22.cc',
4+
'wrap_DanilovEnvelopeSolver20.cc',
5+
'wrap_DanilovEnvelopeSolver22.cc',
6+
'wrap_danilov_envelope.cc',
7+
'danilov_envelope.cc',
8+
])
9+
10+
11+
python.extension_module('danilov_envelope',
12+
sources: ext_sources,
13+
include_directories: inc,
14+
cpp_args: ['-fPIC', '-std=c++11'],
15+
dependencies: [core_dep],
16+
install: true,
17+
subdir: 'orbit/ext',
18+
)

0 commit comments

Comments
 (0)