Modification de `ParticleTest'.

15 files changed, 396 insertions, 24 deletions

diff --git a/src/gui/GLWidget.cc b/src/gui/GLWidget.cc
new file mode 100644
index 0000000..60d411b
--- /dev/null
+++ b/src/gui/GLWidget.cc
@@ -0,0 +1,78 @@
+#include <QtOpenGL>
+#include "GLWidget.h"
+
+GLWidget::GLWidget (QWidget* parent)
+	: QGLWidget (parent)
+{
+	resize (sizeHint ());
+}
+
+GLWidget::~GLWidget () {}
+
+QSize GLWidget::minimumSizeHint () const {
+	return QSize (50, 50);
+}
+
+QSize GLWidget::sizeHint () const {
+	return QSize (800, 600);
+}
+
+void GLWidget::initializeGL () {
+	glClearColor (255, 255, 255, 1.0);
+	glEnable (GL_DEPTH_TEST);
+	gluPerspective(65.0, 4.0/3, 1.0, 1000.0);
+}
+
+void GLWidget::paintGL () {
+	glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+	glLoadIdentity ();
+	//gluLookAt(0,0,0,
+	//	0,0,0,
+	//	0,0,1);
+	//glTranslated (0.0, 0.0, -10.0);
+	//glColor3d (0.0, 0.0, 1.0);
+	//glScaled (300.0, 300.0, 300.0);
+
+	glBegin (GL_POLYGON);
+	glColor3d(1,0,0);
+	glVertex3d (-300, -300, 0.0);
+	glColor3d(0,1,0);
+	glVertex3d (300 + 10, -300, 0.0);
+	glColor3d(0,0,1);
+	glVertex3d (0, 300, 0.0);
+	glColor3d(1,1,1);
+	glVertex3d (0, 0, 100);
+	glColor3d(0,1,0);
+	glVertex3d (300 + 10, -300, 0.0);
+	glEnd ();
+}
+
+void GLWidget::resizeGL (int width, int height) {
+	glViewport (0, 0, width, height);
+	glMatrixMode (GL_PROJECTION);
+	glLoadIdentity ();
+	glOrtho (-width, width, -height, height, -100.0, 100.0);
+	glMatrixMode (GL_MODELVIEW);
+}
+
+void GLWidget::mousePressEvent (QMouseEvent* event) {
+}
+
+void GLWidget::mouseReleaseEvent (QMouseEvent* event) {
+}
+
+void GLWidget::keyPressEvent (QKeyEvent* event) {
+	switch (event->key()) {
+		case Qt::Key_Escape:
+			qApp->quit();
+			break;
+		case Qt::Key_Up:
+			qApp->quit();
+			break;
+		default:
+			break;
+	}
+}
+
+void GLWidget::keyReleaseEvent (QKeyEvent* event) {
+}
diff --git a/src/gui/GLWidget.h b/src/gui/GLWidget.h
new file mode 100644
index 0000000..b2cd02d
--- /dev/null
+++ b/src/gui/GLWidget.h
@@ -0,0 +1,27 @@
+#ifndef GLWIDGET_H
+#define GLWIDGET_H
+
+#include <QGLWidget>
+class GLWidget : public QGLWidget
+{
+	Q_OBJECT
+
+public:
+	GLWidget (QWidget* parent = 0);
+	~GLWidget ();
+
+	QSize minimumSizeHint () const;
+	QSize sizeHint () const;
+
+protected:
+	
+	void initializeGL ();
+	void paintGL ();
+	void resizeGL (int width, int height);
+	void mousePressEvent (QMouseEvent* event);
+	void mouseReleaseEvent (QMouseEvent* event);
+	void keyPressEvent (QKeyEvent* event);
+	void keyReleaseEvent (QKeyEvent* event);
+};
+
+#endif
diff --git a/src/gui/Makefile b/src/gui/Makefile
index 8595ac2..aed836b 100644
--- a/src/gui/Makefile
+++ b/src/gui/Makefile
@@ -1,6 +1,6 @@
 #############################################################################
 # Makefile for building: $(BINDIR)/gui/gui
-# Generated by qmake (2.01a) (Qt 4.7.0) on: Tue Mar 29 19:14:31 2011
+# Generated by qmake (2.01a) (Qt 4.7.0) on: Tue Mar 29 20:46:53 2011
 # Project:  gui.pro
 # Template: app
 # Command: /usr/bin/qmake -o Makefile gui.pro
diff --git a/src/gui/moc_GLWidget.cpp b/src/gui/moc_GLWidget.cpp
new file mode 100644
index 0000000..bef8a88
--- /dev/null
+++ b/src/gui/moc_GLWidget.cpp
@@ -0,0 +1,69 @@
+/****************************************************************************
+** Meta object code from reading C++ file 'GLWidget.h'
+**
+** Created: Mon Mar 28 20:36:46 2011
+**      by: The Qt Meta Object Compiler version 62 (Qt 4.7.0)
+**
+** WARNING! All changes made in this file will be lost!
+*****************************************************************************/
+
+#include "GLWidget.h"
+#if !defined(Q_MOC_OUTPUT_REVISION)
+#error "The header file 'GLWidget.h' doesn't include <QObject>."
+#elif Q_MOC_OUTPUT_REVISION != 62
+#error "This file was generated using the moc from 4.7.0. It"
+#error "cannot be used with the include files from this version of Qt."
+#error "(The moc has changed too much.)"
+#endif
+
+QT_BEGIN_MOC_NAMESPACE
+static const uint qt_meta_data_GLWidget[] = {
+
+ // content:
+       5,       // revision
+       0,       // classname
+       0,    0, // classinfo
+       0,    0, // methods
+       0,    0, // properties
+       0,    0, // enums/sets
+       0,    0, // constructors
+       0,       // flags
+       0,       // signalCount
+
+       0        // eod
+};
+
+static const char qt_meta_stringdata_GLWidget[] = {
+    "GLWidget\0"
+};
+
+const QMetaObject GLWidget::staticMetaObject = {
+    { &QGLWidget::staticMetaObject, qt_meta_stringdata_GLWidget,
+      qt_meta_data_GLWidget, 0 }
+};
+
+#ifdef Q_NO_DATA_RELOCATION
+const QMetaObject &GLWidget::getStaticMetaObject() { return staticMetaObject; }
+#endif //Q_NO_DATA_RELOCATION
+
+const QMetaObject *GLWidget::metaObject() const
+{
+    return QObject::d_ptr->metaObject ? QObject::d_ptr->metaObject : &staticMetaObject;
+}
+
+void *GLWidget::qt_metacast(const char *_clname)
+{
+    if (!_clname) return 0;
+    if (!strcmp(_clname, qt_meta_stringdata_GLWidget))
+        return static_cast<void*>(const_cast< GLWidget*>(this));
+    return QGLWidget::qt_metacast(_clname);
+}
+
+int GLWidget::qt_metacall(QMetaObject::Call _c, int _id, void **_a)
+{
+    _id = QGLWidget::qt_metacall(_c, _id, _a);
+    if (_id < 0)
+        return _id;
+    return _id;
+}
+QT_END_MOC_NAMESPACE
diff --git a/src/main/CurvedElement.cc b/src/main/CurvedElement.cc
index d6ceff6..16816d9 100644
--- a/src/main/CurvedElement.cc
+++ b/src/main/CurvedElement.cc
@@ -7,7 +7,7 @@
 
 #include <assert.h>
 #include <math.h>
-#include "exceptions.h"
+#include "exception.h"
 #include "CurvedElement.h"
 
 using namespace std;
diff --git a/src/main/CurvedElement.h b/src/main/CurvedElement.h
index e65bcc5..7bbb71c 100644
--- a/src/main/CurvedElement.h
+++ b/src/main/CurvedElement.h
@@ -45,7 +45,8 @@ public:
 
 	virtual bool isOutside(const Particle& particle) const {
 		Vector3D x(particle.getPosition() - entryPosition);
-		return (x - Vector3D(x.getX(), x.getY(), 0).unit() / fabs(curvature)).norm() > sectionRadius;
+		if (x == Vector3D::Null) return false;
+		else return (x - Vector3D(x.getX(), x.getY(), 0).unit() / fabs(curvature)).norm() > sectionRadius;
 	}
 
 	virtual bool isPast(const Particle& particle) const {
diff --git a/src/main/FODO.cc b/src/main/FODO.cc
new file mode 100644
index 0000000..c87324c
--- /dev/null
+++ b/src/main/FODO.cc
@@ -0,0 +1,52 @@
+/*
+ * FODO.cc
+ *
+ *  Created on: Mar 29, 2011
+ *      Author: jakob
+ */
+
+#include "FODO.h"
+#include "CompositeElement.h"
+#include "StraightElement.h"
+#include "Quadrupole.h"
+#include "Vector3D.h"
+#include "exception.h"
+
+namespace vhc {
+
+FODO::FODO(const Vector3D& entry, const Vector3D& exit, double sectionRadius, double straightLength, double focalisingCoefficient, Element* next):
+	CompositeElement(entry, exit, sectionRadius, next),
+	focalisingQuadrupole(NULL),
+	defocalisingQuadrupole(NULL),
+	straightElement1(NULL),
+	straightElement2(NULL) {
+
+	double L = straightLength;
+	double l = getDiagonal().norm() / 2 - L;
+	if (l < 0) throw IllegalArgumentException("Length of straight elements must be less than half of FODO length.");
+
+	Vector3D d = getDiagonal().unit();
+	focalisingQuadrupole = new Quadrupole(entry, entry + d * l, sectionRadius, focalisingCoefficient);
+	elements.push_back(focalisingQuadrupole);
+	straightElement1 = new StraightElement(focalisingQuadrupole->getExitPosition(), focalisingQuadrupole->getExitPosition() + d * L, sectionRadius);
+	elements.push_back(straightElement1);
+	defocalisingQuadrupole = new Quadrupole(straightElement1->getExitPosition(), straightElement1->getExitPosition() + d * l, sectionRadius, -focalisingCoefficient);
+	elements.push_back(defocalisingQuadrupole);
+	straightElement2 = new StraightElement(defocalisingQuadrupole->getExitPosition(), defocalisingQuadrupole->getExitPosition() + d * L, sectionRadius);
+	elements.push_back(straightElement2);
+
+	//connection des elements
+	for (int i(0); i < elements.size() - 1; i++) {
+		elements[i]->setNext(elements[i+1]);
+	}
+
+}
+
+FODO::~FODO() {
+	for (int i(0); i < elements.size(); i++) {
+		delete elements[i];
+		elements[i] = NULL;
+	}
+}
+
+}
diff --git a/src/main/FODO.h b/src/main/FODO.h
new file mode 100644
index 0000000..99241a6
--- /dev/null
+++ b/src/main/FODO.h
@@ -0,0 +1,31 @@
+/*
+ * FODO.h
+ *
+ *  Created on: Mar 29, 2011
+ *      Author: jakob
+ */
+
+#ifndef FODO_H_
+#define FODO_H_
+
+#include "Quadrupole.h"
+#include "CompositeElement.h"
+#include "Vector3D.h"
+
+namespace vhc {
+
+class FODO: public CompositeElement {
+	Quadrupole* focalisingQuadrupole;
+	Quadrupole* defocalisingQuadrupole;
+	StraightElement* straightElement1;
+	StraightElement* straightElement2;
+
+public:
+	FODO(const Vector3D& entry, const Vector3D& exit, double sectionRadius, double straightLength, double focalisingCoefficient,  Element* next = NULL);
+
+	virtual ~FODO();
+};
+
+}
+
+#endif /* FODO_H_ */
diff --git a/src/main/Particle.h b/src/main/Particle.h
index d2cc7d8..9c0ba1a 100644
--- a/src/main/Particle.h
+++ b/src/main/Particle.h
@@ -38,7 +38,7 @@ private:
 	/** Charge de cette particule. */
 	double charge;
 
-	double energy;
+	double gamma;
 
 public:
 
@@ -51,15 +51,14 @@ public:
 			charge(charge)
 	{};
 
-	Particle(const Vector3D& position0, double energy, const Vector3D& direction, double mass, double charge):
+	Particle(const Vector3D& position, double mass, double charge, double energy, const Vector3D& direction):
 		position(position),
-		velocity(0, 0, 0),
+		velocity(constants::c * sqrt(1 - (mass * mass) / (energy * energy)) * direction),
 		force(0, 0, 0),
 		mass(mass),
-		charge(charge) {
-
-		velocity = constants::c * sqrt(1 - (mass * mass) / (energy * energy)) * direction;
-	};
+		charge(charge),
+		gamma(energy/(mass*constants::c2))
+	{};
 
 
 	/** Retourne la position de cette particule. */
@@ -77,6 +76,13 @@ public:
 	/** Applique une force sur cette particule. */
 	void applyForce(const Vector3D& f) {force = force + f;}
 
+	void applyMagneticForce(const Vector3D& b, double dt) {
+		if (b != Vector3D::Null) {
+			Vector3D f = charge * velocity.cross(b);
+			force = force + f.rotate(velocity.cross(f), (dt * f.norm()) / (2 * gamma * mass * velocity.norm()));
+		}
+	}
+
 	/** Retourne la masse de cette particule. */
 	double getMass() const {return mass;}
 
@@ -86,10 +92,12 @@ public:
 	/** Retourne la vitesse de cette particule. */
 	Vector3D getVelocity() const {return velocity;}
 
+	void setVelocity(const Vector3D& v) {velocity = v;}
+
 	//GeV
-	double getEnergy() const {return energy;}
+	double getEnergy() const {return gamma * mass * constants::c2;}
 
-	double getGamma() const {return energy / (mass * constants::c2);}
+	double getGamma() const {return gamma;}
 
 	/** Retourne une représentation en chaîne de cette particule. */
 	virtual std::string toString() const {
@@ -99,7 +107,7 @@ public:
 		s << "\tvelocity: "	<< velocity	<< "\n";
 		s << "\tmass:     "	<< mass		<< "\n";
 		s << "\tcharge:   "	<< charge	<< "\n";
-		s << "\tforce:    "	<< force	<< "\n";
+		s << "\tforce:    "	<< force;
 		return s.str();
 	}
 
diff --git a/src/main/Quadrupole.cc b/src/main/Quadrupole.cc
new file mode 100644
index 0000000..7236150
--- /dev/null
+++ b/src/main/Quadrupole.cc
@@ -0,0 +1,14 @@
+/*
+ * Quadrupole.cc
+ *
+ *  Created on: Mar 29, 2011
+ *      Author: jakob
+ */
+
+#include "Quadrupole.h"
+
+namespace vhc {
+
+
+
+}
diff --git a/src/main/Quadrupole.h b/src/main/Quadrupole.h
new file mode 100644
index 0000000..017a543
--- /dev/null
+++ b/src/main/Quadrupole.h
@@ -0,0 +1,51 @@
+/*
+ * Quadrupole.h
+ *
+ *  Created on: Mar 29, 2011
+ *      Author: jakob
+ */
+
+#ifndef QUADRUPOLE_H_
+#define QUADRUPOLE_H_
+
+#include <string>
+#include "StraightElement.h"
+#include "Vector3D.h"
+
+namespace vhc {
+
+class Quadrupole: public StraightElement {
+
+private:
+	double focusingCoefficient;
+
+public:
+	Quadrupole(const Vector3D& entry, const Vector3D& exit, double sectionRadius, double focusingCoefficient, Element* next = NULL):
+		StraightElement(entry, exit, sectionRadius, next),
+		focusingCoefficient(focusingCoefficient)
+	{};
+
+	virtual ~Quadrupole() {};
+
+	virtual Vector3D magneticFieldAt(const Vector3D& position) {
+		Vector3D x = position - getEntryPosition();
+		Vector3D d = getDiagonal().unit();
+		Vector3D y = x - x.dot(d) * d;
+		Vector3D u = Vector3D::k.cross(d);
+		return focusingCoefficient * (y.dot(u) * Vector3D::k + position.getZ() * u);
+	}
+
+	double getFocusingCoefficient() const {
+		return focusingCoefficient;
+	}
+
+	void setFocusingCoefficient(double value) {
+		focusingCoefficient = value;
+	}
+
+	virtual std::string getType() const {return "Quadrupole";}
+};
+
+}
+
+#endif /* QUADRUPOLE_H_ */
diff --git a/src/main/Vector3D.h b/src/main/Vector3D.h
index 6d4dfed..979c9a2 100644
--- a/src/main/Vector3D.h
+++ b/src/main/Vector3D.h
@@ -10,7 +10,7 @@
 
 #include <sstream>
 #include <math.h>
-#include "exceptions.h"
+#include "exception.h"
 #include "Printable.h"
 
 namespace vhc {
@@ -113,7 +113,7 @@ public:
 
 	/** Rotation vectorielle. Retourne le vecteur courant, évalué dans la formule en <code>a</code>,
 	 * le vecteur de l'axe, et en <code>t</code>, l'angle de rotation. */
-	Vector3D rotate(const Vector3D& axis,const double& t) const {
+	Vector3D rotate(const Vector3D& axis, double t) const {
 		const Vector3D& x = *this;
 		const Vector3D& a = ~axis;
 
diff --git a/src/main/exceptions.h b/src/main/exception.h
index 944af4c..720cbe3 100644
--- a/src/main/exceptions.h
+++ b/src/main/exception.h
@@ -1,12 +1,12 @@
 /*
- * exceptions.h
+ * exception.h
  *
  *  Created on: Mar 24, 2011
  *      Author: jakob
  */
 
-#ifndef EXCEPTIONS_H_
-#define EXCEPTIONS_H_
+#ifndef EXCEPTION_H_
+#define EXCEPTION_H_
 
 #include <sstream>
 #include <string>
@@ -74,4 +74,4 @@ public:
 
 }
 
-#endif /* EXCEPTIONS_H_ */
+#endif /* EXCEPTION_H_ */
diff --git a/src/test/ParticleTest.cc b/src/test/ParticleTest.cc
index 73632bc..8b65027 100644
--- a/src/test/ParticleTest.cc
+++ b/src/test/ParticleTest.cc
@@ -5,16 +5,57 @@
  *      Author: jakob
  */
 
-#include "Particle.h"
 #include <iostream>
 #include <string>
+#include <vector>
+#include "StraightElement.h"
+#include "Dipole.h"
+#include "Quadrupole.h"
+#include "Particle.h"
 #include "Vector3D.h"
 
 using namespace vhc;
 using namespace std;
-/** lance le test*/
+
+Particle particle(Vector3D(0,0,0), 0.2, constants::e, 450, Vector3D::i);
+Element* element;
+
+
+void step(double h);
+
 int main() {
-	Particle p(Vector3D::Null, 0, 0);
-	cout << p << endl;
+	//Particle p(Vector3D::Null, 0, 0);
+	//cout << p << endl;
+	element = new Dipole(Vector3D::Null, Vector3D::i-Vector3D::j, 0.2, 1.0, Vector3D(0, 0, 1) * particle.getGamma() * particle.getMass() * 1 * particle.getVelocity().norm() / particle.getCharge());
+	cout << *element << endl;
+
+	double t = 0;
+	double h = 1E-9;
+	char c('0');
+	bool hit = false;
+	do {
+		//if (element->isOutside(particle)) hit = true;
+		cout << particle << endl;
+		step(h);
+		t += h;
+		cin.get(c);
+	} while (c != 'x' && !hit);
+
+	if (hit) {
+		cout << "Particle hit the wall!!!" << endl;
+	}
+
+
+	delete element;
 	return 0;
 }
+
+void step(double h) {
+	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);
+
+}
+
diff --git a/src/test/Vector3DTest.cc b/src/test/Vector3DTest.cc
index 5e1fde9..0065f00 100644
--- a/src/test/Vector3DTest.cc
+++ b/src/test/Vector3DTest.cc
@@ -9,7 +9,7 @@
 #include <assert.h>
 #include <iomanip>
 #include <limits>
-#include "exceptions.h"
+#include "exception.h"
 #include "Vector3D.h"
 
 using namespace std;