1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
/*
* Bunch.cc
*
* Created on: May 23, 2011
* Author: jakob
*/
#include <math.h>
#include "Bunch.h"
#include "random.h"
#include "Vector3D.h"
#include "Element.h"
namespace vhc {
Bunch::Bunch(const Particle& referenceParticle, int quantity, int lambda, double standardDeviation, double length, double targetEmittance, double A12, double A22):
Beam(referenceParticle, quantity, lambda),
standardDeviation(standardDeviation),
length(length),
targetEmittance(targetEmittance),
A12(A12),
A22(A22) {}
Bunch::~Bunch() {
}
void Bunch::initializeParticles() {
create(getLength());
}
void Bunch::create(double dt) {
random::seed(0);
double debit = quantity / length;
double fraction(debit * dt); // debit "vrai", mais a priori non entier
int number(fraction); // partie entière
fraction -= number; // partie fractionnaire
// on ajoute 1 au hasard, proportionnellement à la partie fractionnaire :
if ( random::uniform(0.0, 1.0) < fraction ) ++number;
//calculation des proprietes des particules
for (int i = 0; i < quantity / lambda; ++i) {
Particle* particle = getReferenceParticle().clone();
double sa[2]; // = [r, z]
double va[2]; // = [vr, vz]
//calculation de r, vr, z et vz
for (int i = 0; i < 2; ++i) {
double A11 = getA11();
double A12 = getA12();
double A22 = getA22();
double theta = atan(2 * A12 / (A11 - A22));
double a = A11 * cos(theta) * cos(theta) + 2 * A12 * cos(theta) * sin(theta) + A22 * sin(theta) * sin(theta);
double b = A11 * sin(theta) * sin(theta) + 2 * A12 * cos(theta) * sin(theta) + A22 * cos(theta) * cos(theta);
double x = random::gaussian(0, sqrt(getTargetEmittance() / a));
double y = random::gaussian(0, sqrt(getTargetEmittance() / b));
sa[i] = cos(theta) * x + sin(theta) * y;
va[i] = -sin(theta) * x + cos(theta) * y;
}
double r = sa[0];
double vr = va[0];
double z = sa[1];
double vz = va[1];
double norm = random::gaussian(getReferenceParticle().getVelocity().norm(), getStandardDeviation());
Vector3D u = getReferenceParticle().getElement()->getHorizontalAt(getReferenceParticle().getPosition());
Vector3D v = vr * u + vz * Vector3D::k + sqrt(norm * norm - vr * vr - vz * vz) * u.cross(Vector3D::k);
//composante s de position
double s = random::gaussian(0, getLength());
//ecart de position
Vector3D dp = r * u + z * Vector3D::k + s * u.cross(Vector3D::k);
particle->setMass(getReferenceParticle().getMass() * getLambda());
particle->setCharge(getReferenceParticle().getCharge() * getLambda());
particle->translate(dp);
particles.push_back(particle);
}
}
Bunch* Bunch::clone() const {
return new Bunch(*this);
}
double Bunch::getA11() const {
return (1 + getA12() * getA12() ) / getA22();
}
double Bunch::getA12() const {
return A12;
}
double Bunch::getA22() const {
return A22;
}
double Bunch::getStandardDeviation() const {
return standardDeviation;
}
double Bunch::getLength() const {
return length;
}
double Bunch::getTargetEmittance() const {
return targetEmittance;
}
}
|
