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
|
/*
* ParticleTest.cc
*
* Created on: Mar 16, 2011
* Author: jakob
*/
#include <iostream>
#include <string>
#include <vector>
#include "constants.h"
#include "StraightElement.h"
#include "Dipole.h"
#include "Quadrupole.h"
#include "Particle.h"
#include "Vector3D.h"
using namespace vhc;
using namespace std;
Element* element = NULL;
Particle* particle = NULL;
void step(double h);
double toKg(double mGeV) {
return constants::E / constants::C2 * 1E9 * mGeV;
}
int main() {
Vector3D entry = Vector3D(0, 10, 0);
Vector3D exit = Vector3D(10, 0, 0);
double sectionRadius = 0.2;
double curvature = 0.1;
Vector3D direction = entry.cross(Vector3D::k);
double mass = 9.11E-31;
double charge = constants::E;
double energy = 1 * 1E9 * constants::E;
particle = new Particle(entry, mass, charge, energy, direction);
double Bz = particle->getGamma() * particle->getMass() * curvature * particle->getVelocity().norm() / particle->getCharge();
cout << Bz << endl;
element = new Dipole(entry, exit, sectionRadius, curvature, Vector3D::k * Bz);
cout << *element << endl;
double t = 0;
double h = 1E-12;
char c('0');
bool hit = false;
do {
if (element->hasHit(*particle)) hit = true;
cout << "t = " << t << endl;
step(h);
t += h;
cin.get(c);
} while (c != 'c' && !hit);
if (hit) {
cout << "Particle hit the wall!!!" << endl;
}
delete particle;
delete element;
return 0;
}
void step(double h) {
cout << *particle << endl;
particle->applyMagneticForce(element->magneticFieldAt(particle->getPosition()), h);
Vector3D a = particle->getForce() / (particle->getGamma() * particle->getMass());
particle->setVelocity(particle->getVelocity() + a * h);
particle->setPosition(particle->getPosition() + particle->getVelocity() * h);
particle->setForce(Vector3D::Null);
}
|
