Repare bug dans la fermeture de l'accelerateur. Ajout d'un test de performance de l'accelerateur.

14 files changed, 333 insertions, 83 deletions

diff --git a/Makefile b/Makefile
index ae25778..da10e0f 100644
--- a/Makefile
+++ b/Makefile
@@ -28,7 +28,7 @@ export BINDIR = $(BASEDIR)/bin
 # CXXFLAGS += -g                      # pour debugger
 # CXXFLAGS += -pg                     # pour profiler
 # LDFLAGS  += -pg                     # pour profiler
-# CXXFLAGS += -O2                     # pour optimiser la vitesse
+ CXXFLAGS += -O2                     # pour optimiser la vitesse
 
 export CXXFLAGS
 
@@ -51,6 +51,11 @@ gui-build: build
 	cd $(SRCDIR)/gui; qmake; make all
 	@echo "Done building GUI."
 
+gui-build-noqmake: build
+	@echo "Building GUI..."
+	cd $(SRCDIR)/gui; make all
+	@echo "Done building GUI."
+
 # Genere la documentation
 doc:
 	@echo "Building doc..."
diff --git a/src/gui/Main.cc b/src/gui/Main.cc
index 79242ca..3925f61 100644
--- a/src/gui/Main.cc
+++ b/src/gui/Main.cc
@@ -23,7 +23,7 @@
 using namespace std;
 using namespace vhc;
 
-std::vector< Particle > createParticles(const Vector3D& position, int n, double mass = constants::ELECTRON_MASS, double charge = constants::E, double energy = 1E9 * constants::E, Vector3D direction = Vector3D::i) {
+std::vector< Particle > createParticles(const Vector3D& position, int n, double mass = constants::PROTON_MASS, double charge = constants::E, double energy = 2E9 * constants::E, Vector3D direction = -Vector3D::j) {
 	std::vector< Particle > v;
 
 	double r = 0.1;
@@ -105,11 +105,28 @@ Une particule :
 	acc->add(e7);
 	acc->add(e8);
 
+	acc->close();
+
+	//proton
 	Particle p1 = Particle(Vector3D(3.00, 0, 0), constants::PROTON_MASS, constants::E, 2 * constants::GeV, -Vector3D::j);
 	Particle p2 = Particle(Vector3D(2.99, 0, 0), constants::PROTON_MASS, constants::E, 2 * constants::GeV, -Vector3D::j);
 	acc->add(p1);
 	acc->add(p2);
-	acc->close();
+
+	//anti-proton
+	Particle ap1 = Particle(Vector3D(3.00, 0, 0), constants::PROTON_MASS, -constants::E, 2 * constants::GeV, Vector3D::j);
+	Particle ap2 = Particle(Vector3D(2.99, 0, 0), constants::PROTON_MASS, -constants::E, 2 * constants::GeV, Vector3D::j);
+	acc->add(ap1);
+	acc->add(ap2);
+
+
+
+	/*std::vector< Particle > ps = createParticles(e1.getEntryPosition(), 1000);
+
+	for (int i = 0; i < ps.size(); ++i) {
+		acc->add(ps[i]);
+	}*/
+
 
 	return acc;
 }
@@ -201,8 +218,6 @@ int main(int argc, char *argv[])
 
     vhc::Stage window;
     Accelerator* acc = standard();
-    cout << (*acc) << endl;
-    cout.flush();
     window.accelerator = acc;
 
     window.setWindowTitle("Virtual Hadron Collider");
diff --git a/src/gui/ParticleRenderer.cc b/src/gui/ParticleRenderer.cc
index 3c8fe86..05f1b63 100644
--- a/src/gui/ParticleRenderer.cc
+++ b/src/gui/ParticleRenderer.cc
@@ -9,9 +9,11 @@
 #include "ParticleRenderer.h"
 #include "util.h"
 
+using namespace vhc::util;
+
 namespace vhc {
 
-ParticleRenderer::ParticleRenderer() {
+ParticleRenderer::ParticleRenderer(): _drawSpheres(false) {
 
 }
 
@@ -22,11 +24,20 @@ ParticleRenderer::~ParticleRenderer() {
 void ParticleRenderer::render(const Particle& particle) const {
 	glPushMatrix();
 	glTranslated(particle.getPosition().getX(), particle.getPosition().getY(), particle.getPosition().getZ());
-	glBegin(GL_POINTS);
-	//glVertex3d(0, 0, 0);
-	util::sphere(0.01);
-	glEnd();
+
+	if (_drawSpheres) {
+		glBegin(GL_POINTS);
+		glVertex3d(0, 0, 0);
+		glEnd();
+	}
+	else {
+		util::sphere(0.01);
+	}
 	glPopMatrix();
 }
 
+void ParticleRenderer::drawSpheres(bool value) {
+	_drawSpheres = value;
+}
+
 }
diff --git a/src/gui/ParticleRenderer.h b/src/gui/ParticleRenderer.h
index c1c71b0..ffa6707 100644
--- a/src/gui/ParticleRenderer.h
+++ b/src/gui/ParticleRenderer.h
@@ -10,14 +10,24 @@
 
 #include "Renderer.h"
 #include "Particle.h"
+#include "util.h"
 
 namespace vhc {
 
 class ParticleRenderer: public Renderer<Particle> {
+
+private:
+	bool _drawSpheres;
+
 public:
+
 	ParticleRenderer();
 	virtual ~ParticleRenderer();
 	void render(const Particle& particle) const;
+
+	void drawSpheres(bool value);
+
+
 };
 
 }
diff --git a/src/gui/Stage.cc b/src/gui/Stage.cc
index 13c3411..b3f054a 100644
--- a/src/gui/Stage.cc
+++ b/src/gui/Stage.cc
@@ -25,7 +25,7 @@ Stage::Stage(QWidget* parent):
 		displayMode(FILL),
 		keys(0),
 		frameTime(0),
-		h(1E-12),
+		h(1E-11),
 		paused(true) {
 
 	timer = new QTimer(this);
@@ -58,16 +58,31 @@ void Stage::resizeGL (int width, int height) {
 	glMatrixMode (GL_MODELVIEW);
 }
 
+void Stage::displayText(QString text[], int size) {
+	for (int i = 0; i < size; ++i) {
+		renderText(0, (i + 1) * 12, text[i]);
+	}
+}
+
 void Stage::paintGL() {
 	//time.start();
 	glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 	glLoadIdentity();
 
 	camera.setView();
-	renderText(0,12,QString("fps: ") + QString::number(1E3 / frameTime));
-	renderText(0,24,QString("camera coordinates: ") + camera.getPosition().toString().c_str());
-	renderText(0,36,QString("heading: ") + QString::number(camera.getHeading()));
-	renderText(0,48,QString("pitch: ") + QString::number(camera.getPitch()));
+	QString text[] = {
+			QString("-----simulation-----"),
+			QString("fps: ") + QString::number(1E3 / frameTime),
+			QString("status: ") + (paused ? QString("paused") : QString("running")),
+			QString("-----camera-----"),
+			QString("coordinates: ") + camera.getPosition().toString().c_str(),
+			QString("heading: ") + QString::number(camera.getHeading()),
+			QString("pitch: ") + QString::number(camera.getPitch()),
+			QString("-----accelerator-----"),
+			QString("Elements: ") + QString::number(accelerator->getElements().size()),
+			QString("Particles: ") + QString::number(accelerator->getParticles().size())
+	};
+	displayText(text, 10);
 	//renderText(0,60,QString("") + accelerator->getParticle(0)->getElement()->magneticFieldAt(accelerator->getParticle(0)->getPosition()).toString().c_str());
 	//renderText(0,72,QString("") + accelerator->getParticle(0)->toString().c_str());
 	axes();
@@ -128,7 +143,7 @@ void Stage::paintGL() {
 		camera.move(mv);
 	}
 
-	if (!paused) for (int i = 0; i < 100; ++i) accelerator->step(h);
+	if (!paused) accelerator->step(h);
 
 	glColor3d(1,1,0);
 	util::crosshair();
diff --git a/src/gui/Stage.h b/src/gui/Stage.h
index 445ce63..7eba62a 100644
--- a/src/gui/Stage.h
+++ b/src/gui/Stage.h
@@ -16,16 +16,11 @@
 #include "ParticleRenderer.h"
 #include "Element.h"
 #include "Accelerator.h"
+#include "util.h"
 
 
 namespace vhc {
 
-enum DisplayMode {
-	FILL = 0,
-	WIREFRAME = 1,
-	POINTS = 2
-};
-
 class Stage: public QGLWidget {
 
 	Q_OBJECT
@@ -40,6 +35,8 @@ public:
 
 protected:
 
+	void displayText(QString text[], int size);
+
 	void initializeGL();
 	void resizeGL (int width, int height);
 	void paintGL();
@@ -64,7 +61,7 @@ private:
 	ElementRenderer elementRenderer;
 	ParticleRenderer particleRenderer;
 
-	DisplayMode displayMode;
+	util::DisplayMode displayMode;
 
 	int keys;
 	QPoint center;
diff --git a/src/gui/gui.pro b/src/gui/gui.pro
index 3b68f40..0a73f23 100644
--- a/src/gui/gui.pro
+++ b/src/gui/gui.pro
@@ -9,6 +9,11 @@ INCLUDEPATH += $(SRCDIR)/main/
 LIBS += -L$(BINDIR)/main -lvhc
 QT += opengl
 
+#CONFIG += DEBUG
+#QMAKE_CXXFLAGS_DEBUG += -pg
+#QMAKE_LFLAGS_DEBUG += -pg
+
+
 # Input
 HEADERS += Camera.h ElementRenderer.h util.h Stage.h ParticleRenderer.h Renderer.h
 SOURCES += Camera.cc ElementRenderer.cc util.cc Stage.cc ParticleRenderer.cc Main.cc 
diff --git a/src/gui/util.h b/src/gui/util.h
index 866ef96..9cfded6 100644
--- a/src/gui/util.h
+++ b/src/gui/util.h
@@ -13,6 +13,12 @@ namespace vhc {
 /** Contient des fonctions utilitaires. */
 namespace util {
 
+enum DisplayMode {
+	FILL = 0,
+	WIREFRAME = 1,
+	POINTS = 2
+};
+
 /** Dessine un tore autour de l'origine, sur le plan xy.
  *  @param R rayon du tore
  *  @param r rayon du "cylindre" du tore
diff --git a/src/main/Accelerator.cc b/src/main/Accelerator.cc
index 124cec1..37f8dbb 100644
--- a/src/main/Accelerator.cc
+++ b/src/main/Accelerator.cc
@@ -13,96 +13,135 @@ using namespace std;
 
 namespace vhc {
 
+Accelerator::Accelerator():
+		elementCollec(0),
+		particleCollec(0),
+		allowLinear(false),
+		closed(false)
+	{};
+
+Accelerator::~Accelerator() {
+	clear();
+}
+
+void Accelerator::initializeParticles() {
+	//rajouter les particules dans leurs elements respectifs
+	for (ParticleIterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
+
+		for (ElementIterator j = elementCollec.begin(); j != elementCollec.end(); ++j) {
+			if ((**j).contains(**i)) {
+				(**i).setElement(*j);
+				break;
+			}
+		}
+
+		//si une particule n'est pas contenue dans un element elle est supprimee
+		if ((**i).getElement() == NULL) {
+			delete *i;
+			i = particleCollec.erase(i);
+			--i;
+			//std::cout << "Particle hit wall. Removed from simulation" << std::endl;
+		}
+	}
+}
+
 Element& Accelerator::add(const Element& element) {
 	Element* e = element.clone();
 	elementCollec.push_back(e);
+	closed = false;
 	return *e;
 }
 
 Particle& Accelerator::add(const Particle& particle) {
 	Particle* p = particle.clone();
 	particleCollec.push_back(p);
+	closed = false;
 	return *p;
 }
 
-const list<Element*>& Accelerator::getElements() const { return elementCollec;}
+const Accelerator::ElementCollection& Accelerator::getElements() const { return elementCollec;}
 
-const list<Particle*>& Accelerator::getParticles() const { return particleCollec;}
+const Accelerator::ParticleCollection& Accelerator::getParticles() const { return particleCollec;}
 
 void Accelerator::close() {
-	for (list<Element*>::iterator current = elementCollec.begin(); current != elementCollec.end(); ++current) {
+	for (ElementIterator current = elementCollec.begin(); current != elementCollec.end(); ++current) {
 
-		list<Element*>::iterator next = current;
+		ElementIterator next = current;
 		++next;
 		if (next == elementCollec.end()) next = elementCollec.begin();
 
 
-		// est-ce que les elements se suivent?
+		// est-ce que les elements se suivent (sont connectes)?
 		if (Vector3D::ae((*current)->getExitPosition(), (*next)->getEntryPosition())) {
 			(**current).setNext(*next);
 			(**next).setPrevious(*current);
-		} else
-			throw UnsupportedOperationException("Cannot close accelerator. Two succeeding elements are not physically connected! (not close enough)");
-	}
 
-	//rajouter les particules dans leurs elements respectifs
-	for (list<Particle*>::iterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
+		//sinon est-ce qu'il s'agit du dernier element?
+		} else if (next == elementCollec.begin()) {
+			if (!allowLinear) throw UnsupportedOperationException("Cannot close accelerator. Linear Accelerators are not allowed.");
 
-		for (list<Element*>::iterator j = elementCollec.begin(); j != elementCollec.end(); ++j) {
-			if ((**j).contains(**i)) {
-				(**i).setElement(*j);
-				break;
-			}
-		}
+		//sinon
+		} else throw UnsupportedOperationException("Cannot close accelerator. Two succeeding elements are not physically connected. (not close enough)");
+	}
 
-		//si une particule n'est pas contenue dans un element elle est supprimee
-		if ((**i).getElement() == NULL) {
-			delete *i;
-			i = particleCollec.erase(i);
-			--i;
-			std::cout << "Particle hit wall. Removed from simulation" << std::endl;
-		}
+	initializeParticles();
 
-	}
+	closed = true;
 
 }
 
-
 void Accelerator::updateParticles() {
-	for (list<Particle*>::iterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
+	for (ParticleIterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
 		Particle& particle = **i;
-		if (particle.getElement()->isBeside(particle)) {
-			std::cout << "Particle hit wall. Removed from simulation" << std::endl;
-			delete *i;
-			i = particleCollec.erase(i);
-			--i;
-		} else if (particle.getElement()->isAfter(particle)) {
-			if (particle.getElement()->getNext() == NULL) throw Exception("Element in accelerator not connected to next.");
+		if (particle.getElement()->isAfter(particle)) {
+			if (particle.getElement()->getNext() == NULL)
+				if (allowLinear) {
+					delete *i;
+					i = particleCollec.erase(i);
+					--i;
+					//cout << "Particle reached end of accelerator. Removed from simulation" << std::endl;
+				}
+				else throw Exception("Element in accelerator not connected to next.");
 			else particle.setElement(particle.getElement()->getNext());
 		} else if (particle.getElement()->isBefore(particle)) {
-			if (particle.getElement()->getPrevious() == NULL) throw Exception("Element in accelerator not connected to previous.");
+			if (particle.getElement()->getPrevious() == NULL)
+				if (allowLinear) {
+					delete *i;
+					i = particleCollec.erase(i);
+					--i;
+					//cout << "Particle reached beginning of accelerator. Removed from simulation" << std::endl;
+				}
+				else throw Exception("Element in accelerator not connected to previous.");
 			else particle.setElement(particle.getElement()->getPrevious());
+		} else if (particle.getElement()->isBeside(particle)) {
+			//std::cout << "Particle hit wall. Removed from simulation" << std::endl;
+			delete *i;
+			i = particleCollec.erase(i);
+			--i;
 		}
 	}
 }
 
 void Accelerator::clear() {
-	for (list<Particle*>::iterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
+	for (ParticleIterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
 		delete *i;
 		*i = NULL;
 	}
 	particleCollec.clear();
 
-	for (list<Element*>::iterator i = elementCollec.begin(); i != elementCollec.end(); ++i) {
+	for (ElementIterator i = elementCollec.begin(); i != elementCollec.end(); ++i) {
 		delete *i;
 		*i = NULL;
 	}
 	elementCollec.clear();
+
+	closed = false;
 }
 
 void Accelerator::step(double dt) {
+	if (!closed) close();
 
-	for (list<Particle*>::iterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
+	for (ParticleIterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
 		Particle& particle = **i;
 
 		particle.applyMagneticForce(particle.getElement()->magneticFieldAt(particle.getPosition()), dt);
@@ -118,6 +157,10 @@ void Accelerator::step(double dt) {
 	updateParticles();
 }
 
+void Accelerator::enableLinear(bool value) {
+	allowLinear = value;
+}
+
 /** Cf. Accelerator.h */
 std::string Accelerator::toString() const {
 	std::stringstream s;
@@ -129,7 +172,7 @@ std::string Accelerator::toString() const {
 		s << elementCollec.front()->toString()<<"\n";
 	} else {
 		s << "This accelerator is made of the " << elementCollec.size() << " following elements :" << "\n";
-		for (list<Element*>::const_iterator i = elementCollec.begin(); i != elementCollec.end(); ++i) {
+		for (ElementCollection::const_iterator i = elementCollec.begin(); i != elementCollec.end(); ++i) {
 			s << (*i)->toString() << "\n";
 		}
 	}
diff --git a/src/main/Accelerator.h b/src/main/Accelerator.h
index d3c26e0..310c26d 100644
--- a/src/main/Accelerator.h
+++ b/src/main/Accelerator.h
@@ -17,7 +17,7 @@ namespace vhc {
 /** Représente un accélérateur. Cette classe contient en particulier une méthode
  *  qui fait évoluer les particules qu'elle contient. */
 class Accelerator: public Printable {
-private :
+private:
 
 	/** Constructeur de copie ne faisant rien. */
 	Accelerator (Accelerator const& other) {}
@@ -27,42 +27,78 @@ private :
 
 protected:
 
+	typedef std::list<Particle*> ParticleCollection;
+	typedef std::list<Element*>  ElementCollection;
+
+	typedef ParticleCollection::iterator ParticleIterator;
+	typedef ElementCollection::iterator ElementIterator;
+
+
 	/** Collection d'elements contenus dans cet accelerateur. */
-	std::list< Element * > elementCollec;
+	ElementCollection elementCollec;
 
 	/** Collection de particules contenus danc cet accelerateur. */
-	std::list< Particle *> particleCollec;
-
+	ParticleCollection particleCollec;
+
+	/** Autorise les accelerateurs lineaires.
+	 *  @see enableLinear */
+	bool allowLinear;
+
+	/** Indique l'etat ferme de l'accelarateur.
+	 *  @see close */
+	bool closed;
+
+	/** Initialize les particules en leur attribuant l'element dans lequel ils sont contenus.
+	 *  Les particules non-contenus sont supprimes de l'accelerateur. */
+	void initializeParticles();
+
+	/** Met a jour les particules en leur attribuant l'element dans lequel ils sont contenus.
+	 *  Contrairement a <code>initializeParticles()</code>, les elements consideres sont:
+	 *  - l'element actuel de la particule
+	 *  - l'element precedent
+	 *  - l'element suivant
+	 *  Si la particule se situe a cote de son element, elle est supprimee de l'accelerateur.
+	 *  Attention: si la particule saute un element, elle est tout de meme consideree comme etant dans l'element suivant (ou precedent)! Ceci
+	 *  peut survenir si un element est trop petit ou si la simulation est faite avec des pas de temps trop grands. */
+	void updateParticles();
 
 public:
 
-	/** Constructeur d'un accélérateur. */
-	Accelerator ():
-		elementCollec(0),
-		particleCollec(0)
-	{};
+	/** Cree un nouveau accelerateur vide.
+	 *  */
+	Accelerator ();
 
-	virtual ~Accelerator() {clear();}
+	virtual ~Accelerator();
 
-	/** Copie un élément dans l'accélérateur. */
+	/** Copie un élément dans l'accélérateur.
+	 *  L'accelerateur est ouvert en ajoutant un element. */
 	Element& add(const Element& element);
 
 	/** Copie une particule dans l'accélérateur. */
 	Particle& add(const Particle& particle);
 
-	/** Retourne la liste d'elements contenus dans cet accelerateur. */
-	const std::list<Element*> & getElements() const;
+	/** Retourne la liste d'elements contenus dans cet accelerateur.
+	 *  <b>ATTENTION:</b> les elements peuvent etre supprimes sans preavis par l'accelerateur! */
+	const ElementCollection & getElements() const;
 
-	/** Retourne la liste des particules contenus dans cet accelerateur. */
-	const std::list<Particle*> & getParticles() const;
+	/** Retourne la liste des particules contenus dans cet accelerateur.
+	 *  <b>ATTENTION:</b> les particules peuvent etre supprimes sans preavis par l'accelerateur! */
+	const ParticleCollection & getParticles() const;
 
-	/** Ferme l'accelerateur. */
+	/** Ferme l'accelerateur.
+	 *  En invoquant cette methode, la continuite des elements est verifiee et les particules sont attribues leurs elements respectifs.
+	 *  @throws UnsupportedOperationException si les elements sauf le dernier ne sont pas continus
+	 *  @throws UnsupportedOperationException si le dernier element n'est pas continu et que les accelerateurs
+	 *  lineaires ne sont pas autorises. */
 	void close();
 
-	/** Efface tous les éléments et les particules. */
+	/** Efface tous les éléments et les particules.
+	 *  Ouvre l'accelerateur. */
 	void clear();
 
-	void updateParticles();
+	/** Autorise ou interdit la linearite de cet accelerateur.
+	 *  Un accelerateur lineaire peut ne pas etre continu en son dernier element. C'est-a-dire que ce n'est pas une boucle fermee. */
+	void enableLinear(bool value);
 
 	/** Fait évoluer l'accélérateur d'un laps de temps dt. */
 	void step(double dt);
diff --git a/src/main/FODO.cc b/src/main/FODO.cc
index f8ed892..d2555d9 100644
--- a/src/main/FODO.cc
+++ b/src/main/FODO.cc
@@ -45,6 +45,9 @@ FODO::FODO(const Vector3D& entry, const Vector3D& exit, double sectionRadius, do
 	for (int i(0); i < elements.size() - 1; i++) {
 		elements[i]->setNext(elements[i+1]);
 	}
+	for (int i(elements.size() - 1); i > 0; i--) {
+		elements[i]->setPrevious(elements[i-1]);
+	}
 
 }
 
diff --git a/src/main/Makefile b/src/main/Makefile
index 236a85f..9caa7f1 100644
--- a/src/main/Makefile
+++ b/src/main/Makefile
@@ -14,8 +14,8 @@ LOCALDIR = main
 # Si un objet nécessite une compilation non-standard (i.e. pas de règle du style Foo.o : Foo.cc Foo.h), rajouter
 # cette règle.
 LOCALOBJS = Vector3D.o Particle.o Printable.o Element.o CurvedElement.o StraightElement.o \
-			CompositeElement.o Dipole.o Quadrupole.o FODO.o ElementVisitor.o Cloneable.o Accelerator.o exceptions.o \
-			Stepper.o
+			CompositeElement.o Dipole.o Quadrupole.o FODO.o ElementVisitor.o Cloneable.o Accelerator.o exceptions.o Stepper.o
+
 OBJS=$(addprefix $(BINDIR)/$(LOCALDIR)/,$(LOCALOBJS))
 .PHONY = all checkdirs lib
 
@@ -36,4 +36,4 @@ lib: $(OBJS)
 
 # Règle implicite qui modifie le répertoire de destination des fichiers .o
 $(BINDIR)/$(LOCALDIR)/%.o: %.cc
-	$(CXX) $(CXXFLAGS) -c $< -o $@
\ No newline at end of file
+	$(CXX) $(CXXFLAGS) -c $< -o $@
diff --git a/src/test/AcceleratorBenchmarkTest.cc b/src/test/AcceleratorBenchmarkTest.cc
new file mode 100644
index 0000000..671cf03
--- /dev/null
+++ b/src/test/AcceleratorBenchmarkTest.cc
@@ -0,0 +1,101 @@
+/*
+ * AcceleratorBenchmarkTest.cc
+ *
+ *  Created on: May 2, 2011
+ *      Author: jakob
+ */
+
+#include <iostream>
+#include <cmath>
+#include <vector>
+#include <stdlib.h>
+#include "exceptions.h"
+#include "Accelerator.h"
+#include "StraightElement.h"
+#include "Dipole.h"
+#include "Particle.h"
+#include "FODO.h"
+#include "Vector3D.h"
+#include "constants.h"
+
+using namespace vhc;
+using namespace std;
+
+std::vector< Particle > createParticles(const Vector3D& position, int n, double mass = constants::ELECTRON_MASS, double charge = constants::E, double energy = 1E9 * constants::E, Vector3D direction = -Vector3D::j) {
+	std::vector< Particle > v;
+
+	double r = 0.1;
+
+	for (int i = 0; i < n; ++i) {
+		double x = (rand() % 1000) / 1000.0 * r;
+		double y = (rand() % 1000) / 1000.0 * sqrt(r * r - x * x);
+		double z = (rand() % 1000) / 1000.0 * sqrt(r * r - y * y - x * x);
+		v.push_back(Particle(position + Vector3D(x, y, z), mass, charge, energy, direction));
+	}
+
+	return v;
+}
+
+Accelerator* standard() {
+	double B = 5.8915820038873;
+	double b = 1.2;
+	FODO e1 = FODO(Vector3D(3, 2, 0), Vector3D(3, -2, 0), 0.1, 1.0, b);
+	Dipole e2 = Dipole(e1.getExitPosition(), Vector3D(2, -3, 0), 0.1, 1, Vector3D(0, 0, B));
+	FODO e3 = FODO(e2.getExitPosition(), Vector3D(-2, -3, 0), 0.1, 1, b);
+	Dipole e4 = Dipole(e3.getExitPosition(), Vector3D(-3, -2, 0), 0.1, 1, Vector3D(0, 0, B));
+	FODO e5 = FODO(e4.getExitPosition(), Vector3D(-3, 2, 0), 0.1, 1.0, b);
+	Dipole e6 = Dipole(e5.getExitPosition(), Vector3D(-2, 3, 0), 0.1, 1, Vector3D(0, 0, B));
+	FODO e7 = FODO(e6.getExitPosition(), Vector3D(2, 3, 0), 0.1, 1.0, b);
+	Dipole e8 = Dipole(e7.getExitPosition(), e1.getEntryPosition(), 0.1, 1, Vector3D(0, 0, B));
+	Accelerator* acc = new Accelerator();
+	acc->add(e1);
+	acc->add(e2);
+	acc->add(e3);
+	acc->add(e4);
+	acc->add(e5);
+	acc->add(e6);
+	acc->add(e7);
+	acc->add(e8);
+
+	acc->close();
+
+	Particle p1 = Particle(Vector3D(3.00, 0, 0), constants::PROTON_MASS, constants::E, 2 * constants::GeV, -Vector3D::j);
+	Particle p2 = Particle(Vector3D(2.99, 0, 0), constants::PROTON_MASS, constants::E, 2 * constants::GeV, -Vector3D::j);
+	acc->add(p1);
+	acc->add(p2);
+
+	std::vector< Particle > ps = createParticles(e1.getEntryPosition(), 10000);
+
+	for (int i = 0; i < ps.size(); ++i) {
+		acc->add(ps[i]);
+	}
+
+
+	return acc;
+}
+
+Accelerator* accelerator;
+
+int main() {
+	accelerator = standard();
+
+	int steps = 1000;
+	double dt = 1E-11;
+
+	cout << "Simulating " << steps << " steps with " << accelerator->getParticles().size() << " particles...";
+	cout.flush();
+	int t0 = clock();
+	for (int i = 0; i < steps; ++i) {
+		accelerator->step(dt);
+	}
+	int t1 = clock() - t0;
+	cout << "DONE" << endl;
+
+	cout << "Time taken: " << t1 << " ticks @ " << CLOCKS_PER_SEC << " ticks/s ~ " << 1.0 * t1 / CLOCKS_PER_SEC << "s" << endl;
+	cout << "Average: " << 1.0 * t1 / CLOCKS_PER_SEC / steps << " s/step" << endl;
+	cout << "Average: " << 1.0 * t1 / CLOCKS_PER_SEC / steps / accelerator->getParticles().size() << " s/step/particle" << endl;
+
+	return 0;
+}
+
+
diff --git a/src/test/P10ExerciceTest.cc b/src/test/P10ExerciceTest.cc
index cc18ec3..dd54f6d 100644
--- a/src/test/P10ExerciceTest.cc
+++ b/src/test/P10ExerciceTest.cc
@@ -9,6 +9,7 @@
 #include "StraightElement.h"
 #include "Quadrupole.h"
 #include "FODO.h"
+#include "exceptions.h"
 #include <iostream>
 #include <string>
 #include <vector>
@@ -69,10 +70,12 @@ int main() {
 	a.add(*p1);
 	a.add(*p2);
 
-	//TODO a.close();
+	a.enableLinear(true);
+	a.close();
 
 	cout << a << endl;
 	a.clear();
 
+
 	return 0;
 }