Rajoute getHorizontalAt() pour tous les elements, en vue des faisceaux

13 files changed, 38 insertions, 9 deletions

diff --git a/src/main/Accelerator.cc b/src/main/Accelerator.cc
index 37f8dbb..ca5de8a 100644
--- a/src/main/Accelerator.cc
+++ b/src/main/Accelerator.cc
@@ -144,13 +144,15 @@ void Accelerator::step(double dt) {
 	for (ParticleIterator i = particleCollec.begin(); i != particleCollec.end(); ++i) {
 		Particle& particle = **i;
 
+		particle.setForce(Vector3D::Null);
+
 		particle.applyMagneticForce(particle.getElement()->magneticFieldAt(particle.getPosition()), dt);
 
 		Vector3D a = particle.getForce() / (particle.getGamma() * particle.getMass());
 		particle.setVelocity(particle.getVelocity() + a * dt);
 
 		particle.translate(particle.getVelocity() * dt);
-		particle.setForce(Vector3D::Null);
+
 	}
 
 
diff --git a/src/main/CompositeElement.cc b/src/main/CompositeElement.cc
index 3e4a216..48e5757 100644
--- a/src/main/CompositeElement.cc
+++ b/src/main/CompositeElement.cc
@@ -48,6 +48,12 @@ Vector3D CompositeElement::electricFieldAt(const Vector3D& position) const {
 	return e;
 }
 
+Vector3D CompositeElement::getHorizontalAt(const Vector3D& position) const {
+	for (int i(0); i < elements.size(); i++) {
+		if (elements[i]->contains(position)) return elements[i]->getHorizontalAt(position);
+	}
+}
+
 void CompositeElement::accept(const ElementVisitor& v) const {
 	for (int i(0); i < elements.size(); ++i) {
 		elements[i]->accept(v);
diff --git a/src/main/CompositeElement.h b/src/main/CompositeElement.h
index c8d2286..60e0a44 100644
--- a/src/main/CompositeElement.h
+++ b/src/main/CompositeElement.h
@@ -42,6 +42,8 @@ public:
 
 	virtual Vector3D electricFieldAt(const Vector3D& position) const;
 
+	virtual Vector3D getHorizontalAt(const Vector3D& position) const;
+
 	virtual void accept(const ElementVisitor& v) const;
 
 	virtual std::string getType() const;
diff --git a/src/main/CurvedElement.cc b/src/main/CurvedElement.cc
index 4368913..d30bff4 100644
--- a/src/main/CurvedElement.cc
+++ b/src/main/CurvedElement.cc
@@ -59,6 +59,11 @@ bool CurvedElement::isAfter(const Vector3D& position) const {
 	return (position - exitPosition).dot(out) > 0;
 }
 
+Vector3D CurvedElement::getHorizontalAt(const Vector3D& position) const {
+	Vector3D X(position - curvatureCenter);
+	return (X - position.getZ() * Vector3D::k).unit();
+}
+
 std::string CurvedElement::getType() const {return "Curved Element";}
 std::string CurvedElement::toString() const {
 	std::stringstream s;
diff --git a/src/main/CurvedElement.h b/src/main/CurvedElement.h
index ac81746..846f279 100644
--- a/src/main/CurvedElement.h
+++ b/src/main/CurvedElement.h
@@ -54,6 +54,8 @@ public:
 
 	virtual bool isAfter(const Vector3D& position) const;
 
+	virtual Vector3D getHorizontalAt(const Vector3D& position) const;
+
 	/** Retourne la courbure. */
 	double getCurvature() const;
 
diff --git a/src/main/Element.h b/src/main/Element.h
index 1f0afd2..3722be2 100644
--- a/src/main/Element.h
+++ b/src/main/Element.h
@@ -74,6 +74,8 @@ public:
 	bool contains(const Vector3D& position) const;
 	bool contains(const Particle& particle) const;
 
+	virtual Vector3D getHorizontalAt(const Vector3D& position) const = 0;
+
 	/** Retourne le champ magnetique, a l'interieur de cet element a la position donnee. */
 	virtual Vector3D magneticFieldAt(const Vector3D& position) const;
 
diff --git a/src/main/Makefile b/src/main/Makefile
index 9caa7f1..83385fd 100644
--- a/src/main/Makefile
+++ b/src/main/Makefile
@@ -14,7 +14,7 @@ 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 Beam.o
 
 OBJS=$(addprefix $(BINDIR)/$(LOCALDIR)/,$(LOCALOBJS))
 .PHONY = all checkdirs lib
diff --git a/src/main/Particle.cc b/src/main/Particle.cc
index 5784625..6e5eb0b 100644
--- a/src/main/Particle.cc
+++ b/src/main/Particle.cc
@@ -44,7 +44,7 @@ void Particle::applyForce(const Vector3D& f) {force += f;}
 void Particle::applyMagneticForce(const Vector3D& b, double dt) {
 	if (dt != 0 && b != Vector3D::Null) {
 		Vector3D f = charge * velocity.cross(b);
-		applyForce(f.rotate(velocity.cross(f), asin((dt * f.norm()) / (2 * gamma * mass * velocity.norm()))));
+		applyForce(f.rotate(velocity.cross(f), asin(1.0 * (dt * f.norm()) / (2 * gamma * mass * velocity.norm()))));
 	}
 }
 
@@ -70,6 +70,7 @@ std::string Particle::toString() const {
 	s << "\tposition: "	<< position	<< "\n";
 	s << "\tvelocity: "	<< velocity;
 	s << "\tnorm: "		<< velocity.norm() << "\n";
+	s << "\tenergy:   "	<< getEnergy()	<< "\n";
 	s << "\tforce:    "	<< force << "\n";
 	s << "\tgamma: "	<< gamma << "\n";
 	s << "\tmass:     "	<< mass		<< "\n";
diff --git a/src/main/Printable.h b/src/main/Printable.h
index c14aec8..6d2129f 100644
--- a/src/main/Printable.h
+++ b/src/main/Printable.h
@@ -8,6 +8,7 @@
 #ifndef PRINTABLE_H_
 #define PRINTABLE_H_
 
+#include <sstream>
 #include <string>
 
 namespace vhc {
@@ -22,6 +23,8 @@ public:
 	/** Retourne une représentation en chaîne de caractères de cet objet imprimable. */
 	virtual std::string toString() const = 0;
 
+	//virtual std::stringstream toStream(std::stringstream s) const = 0;
+
 };
 
 /** Ajoute la représentation en chaîne de <code>p</code> à <code>output</code>. */
diff --git a/src/main/StraightElement.cc b/src/main/StraightElement.cc
index 420e047..05e88c6 100644
--- a/src/main/StraightElement.cc
+++ b/src/main/StraightElement.cc
@@ -16,8 +16,8 @@ StraightElement::StraightElement(const Vector3D& entry, const Vector3D& exit, do
 StraightElement::~StraightElement() {};
 
 bool StraightElement::isBefore(const Vector3D& position) const {
-	const Vector3D v(position - exitPosition);
-	return (-getDiagonal()).dot(v) > getDiagonal().normSquare();
+	const Vector3D v(position - entryPosition);
+	return (-getDiagonal()).dot(v) > 0;//getDiagonal().normSquare();
 }
 
 bool StraightElement::isBeside(const Vector3D& position) const {
@@ -27,8 +27,12 @@ bool StraightElement::isBeside(const Vector3D& position) const {
 };
 
 bool StraightElement::isAfter(const Vector3D& position) const {
-	const Vector3D v(position - entryPosition);
-	return getDiagonal().dot(v) > getDiagonal().normSquare();
+	const Vector3D v(position - exitPosition);
+	return getDiagonal().dot(v) > 0;
+}
+
+Vector3D StraightElement::getHorizontalAt(const Vector3D& position) const {
+	return Vector3D::k.cross(getDiagonal());
 }
 
 std::string StraightElement::getType() const {return "Straight Element";}
diff --git a/src/main/StraightElement.h b/src/main/StraightElement.h
index 9a3a366..14d21aa 100644
--- a/src/main/StraightElement.h
+++ b/src/main/StraightElement.h
@@ -35,6 +35,8 @@ public:
 
 	virtual void accept(const ElementVisitor& v) const;
 
+	virtual Vector3D getHorizontalAt(const Vector3D& position) const;
+
 	virtual StraightElement* clone() const;
 
 };
diff --git a/src/main/Vector3D.cc b/src/main/Vector3D.cc
index fc7fc2b..9328edf 100644
--- a/src/main/Vector3D.cc
+++ b/src/main/Vector3D.cc
@@ -69,7 +69,7 @@ double Vector3D::tripleProduct(const Vector3D& v, const Vector3D& w) const { ret
 Vector3D Vector3D::rotate(const Vector3D& axis, double t) const {
 	const Vector3D& x = *this;
 	const Vector3D& a = ~axis;
-	return x * cos(t) + a * x.dot(a) * (1-cos(t)) +  a.cross(x) * sin(t);
+	return cos(t) * x +  (1-cos(t)) * x.dot(a) * a + sin(t) * a.cross(x);
 }
 
 bool Vector3D::ae(const Vector3D& u, const Vector3D& v, double epsilon) {return (u - v).norm() <= epsilon;}
diff --git a/src/main/constants.h b/src/main/constants.h
index 5753279..8e0331d 100644
--- a/src/main/constants.h
+++ b/src/main/constants.h
@@ -14,7 +14,7 @@ namespace vhc {
 namespace constants {
 
 /** Vitesse de la lumière dans le vide [m/s]. */
-const double C = 299792458;
+const double C = 299792458.0;
 
 /** Vitesse de la lumière dans le vide au carré [m/s]. */
 const double C2 = C * C;