Ajout des bunchs :)

11 files changed, 316 insertions, 10 deletions

diff --git a/src/gui/Main.cc b/src/gui/Main.cc
index 94f68ee..a5adee9 100644
--- a/src/gui/Main.cc
+++ b/src/gui/Main.cc
@@ -22,6 +22,7 @@
 #include "constants.h"
 #include <vector>
 #include "CircularBeam.h"
+#include "Bunch.h"
 
 using namespace std;
 using namespace vhc;
@@ -143,7 +144,15 @@ Une particule :
 	acc->add(ap1);
 	acc->add(ap2);
 
-	acc->add(CircularBeam(p1, 10, 1));
+	acc->add(CircularBeam(p1, 20, 1));
+
+	double emittance = 5E-6; //m
+	double A_12 = 0.02; //1/m
+	double A_22 = 4;//E-19; // s² m-1 (dépend totalement de l'accélérateur)
+	double length = 300E-12 * constants::C;
+	double stdDev = 0.1;
+	acc->add(Bunch(p1, 500, 1, stdDev, length, emittance, A_12, A_22));
+	acc->add(Bunch(ap1, 500, 1, stdDev, length, emittance, A_12, A_22));
 
 	acc->close();
 
diff --git a/src/gui/Stage.cc b/src/gui/Stage.cc
index 020f31a..4c0328d 100644
--- a/src/gui/Stage.cc
+++ b/src/gui/Stage.cc
@@ -25,7 +25,7 @@ Stage::Stage(QWidget* parent):
 		particleRenderer(),
 		displayMode(FILL),
 		frameTime(0),
-		h(1E-11),
+		h(1E-12),
 		paused(true) {
 
 	timer = new QTimer(this);
@@ -145,7 +145,7 @@ void Stage::paintGL() {
 		camera.move(mv);
 
 
-	if (!paused) for (int i = 0; i < 10; ++i) accelerator->step(h);
+	if (!paused) for (int i = 0; i < 30; ++i) accelerator->step(h);
 
 	glColor3d(1,1,0);
 	util::crosshair();
diff --git a/src/main/Beam.h b/src/main/Beam.h
index 70f9f86..7572908 100644
--- a/src/main/Beam.h
+++ b/src/main/Beam.h
@@ -14,6 +14,7 @@
 
 namespace vhc {
 
+/** Classe de base pour les faisceaux. */
 class Beam: public Cloneable {
 
 public:
@@ -21,6 +22,9 @@ public:
 	typedef std::list<Particle*> ParticleCollection;
 
 	/** Cree un nouveau faisceaux.
+	 *  @param refererenceParticle particule de reference (son element ne joue pas de role)
+	 *  @param quantity la quantite de particules a generer
+	 *  @param lambda le facteur de macroparticules (il y aura lambda fois moins de particules dans le faisceau mais a charges et masses lambdas fois plus grands)
 	 *  <b>ATTENTION:</b> un `Beam' est abstrait et n'initialise pas les particules a partir de la reference.
 	 *  C'est au client de faire cela! */
 	Beam(const Particle& referenceParticle, int quantity, int lambda);
@@ -40,8 +44,11 @@ public:
 	 *  peut survenir si un element est trop petit ou si la simulation est faite avec des pas de temps trop grands. */
 	void updateParticles();
 
+	/** Initialise les particules de ce faisceau.
+	 *  Cette methode est appelee par la classe Accelerator. */
 	virtual void initializeParticles() = 0;
 
+	/** Simule ce faisceau pendant un pas de temps dt. */
 	void step(double dt);
 
 	/** Retourne la quantite de particules contenus dans ce faisceau. */
@@ -65,6 +72,7 @@ public:
 	/** Retourne l'emmitance horizontale de ce faisceau. */
 	double getHorizontalEmittance() const;
 
+	/** Vide ce faisceau en supprimant tous les particules. */
 	void clear();
 
 	virtual Beam* clone() const = 0;
@@ -75,22 +83,22 @@ public:
 	//-------------------------------------------------------------------
 
 	/** Retourne coefficient des ellipses de phases vertical. */
-	double getVerticalA11() const;
+	virtual double getVerticalA11() const;
 
 	/** Retourne coefficient des ellipses de phases vertical. */
-	double getVerticalA12() const;
+	virtual double getVerticalA12() const;
 
 	/** Retourne coefficient des ellipses de phases vertical. */
-	double getVerticalA22() const;
+	virtual double getVerticalA22() const;
 
 	/** Retourne coefficient des ellipses de phases horizontal. */
-	double getHorizontalA11() const;
+	virtual double getHorizontalA11() const;
 
 	/** Retourne coefficient des ellipses de phases horizontal. */
-	double getHorizontalA12() const;
+	virtual double getHorizontalA12() const;
 
 	/** Retourne coefficient des ellipses de phases horizontal. */
-	double getHorizontalA22() const;
+	virtual double getHorizontalA22() const;
 
 protected:
 
diff --git a/src/main/Bunch.cc b/src/main/Bunch.cc
new file mode 100644
index 0000000..e32de97
--- /dev/null
+++ b/src/main/Bunch.cc
@@ -0,0 +1,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;
+}
+
+}
diff --git a/src/main/Bunch.h b/src/main/Bunch.h
new file mode 100644
index 0000000..7794312
--- /dev/null
+++ b/src/main/Bunch.h
@@ -0,0 +1,78 @@
+/*
+ * Bunch.h
+ *
+ *  Created on: May 23, 2011
+ *      Author: jakob
+ */
+
+#ifndef BUNCH_H_
+#define BUNCH_H_
+
+#include "Beam.h"
+
+namespace vhc {
+
+/** Un bunch est un faisceau genere par une source de particules. */
+class Bunch: public Beam {
+protected:
+
+	/** Ecart-type voulu. */
+	double standardDeviation;
+
+	/** Longueur de ce bunch. */
+	double length;
+
+	/** Emittance voulue. */
+	double targetEmittance;
+
+	/** Coefficient A12 de l'ellipse de phase. */
+	double A12;
+
+	/** Coefficient A22 de l'ellipse de phase. */
+	double A22;
+
+public:
+
+	/** Cree un nouveau bunch.
+	 *  @param referenceParticle particule de reference
+	 *  @param quantity quantite de particules a generer
+	 *  @param lambda facteur lambda
+	 *  @param standardDeviation ecart-type du bunch
+	 *  @param length longueur du bunch [m]
+	 *  @param targetEmittance emittance voulue
+	 *  @param A12 coefficient A12 de l'ellipse de phase
+	 *  @param A22 coefficient A22 de l'ellipse de phase */
+	Bunch(const Particle& referenceParticle, int quantity, int lambda, double standardDeviation, double length, double targetEmittance, double A12, double A22);
+
+	virtual ~Bunch();
+
+	virtual void initializeParticles();
+
+	virtual Bunch* clone() const;
+
+	/** Cree un bunch genere pendant <code>dt</code>. */
+	void create(double dt);
+
+	/** Retourne le coefficient A12 de l'ellipse de phase. */
+	double getA12() const;
+
+	/** Retourne le coefficient A22 de l'ellipse de phase. */
+	double getA22() const;
+
+	/** Retourne le coefficient A11 de l'ellipse de phase. */
+	double getA11() const;
+
+	/** Retourne l'ecart-type de ce bunch. */
+	double getStandardDeviation() const;
+
+	/** Retourne la longueur de ce faisceau. */
+	double getLength() const;
+
+	/** Retourne l'emittance voulue de ce faisceau. */
+	double getTargetEmittance() const;
+
+};
+
+}
+
+#endif /* BUNCH_H_ */
diff --git a/src/main/CircularBeam.h b/src/main/CircularBeam.h
index 853774d..5112fd9 100644
--- a/src/main/CircularBeam.h
+++ b/src/main/CircularBeam.h
@@ -12,6 +12,7 @@
 
 namespace vhc {
 
+/** Represente un faisceau dont les particules sont unifomement repartis dans une suite d'elements. */
 class CircularBeam: public Beam {
 
 
diff --git a/src/main/Makefile b/src/main/Makefile
index 96f4d98..09ca21e 100644
--- a/src/main/Makefile
+++ b/src/main/Makefile
@@ -15,7 +15,8 @@ LOCALDIR = main
 # 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 Beam.o SingleBeam.o CircularBeam.o
+			Accelerator.o exceptions.o Stepper.o Beam.o SingleBeam.o CircularBeam.o \
+			random.o Bunch.o
 
 OBJS=$(addprefix $(BINDIR)/$(LOCALDIR)/,$(LOCALOBJS))
 .PHONY = all checkdirs lib
diff --git a/src/main/Particle.cc b/src/main/Particle.cc
index 6e5eb0b..79f8672 100644
--- a/src/main/Particle.cc
+++ b/src/main/Particle.cc
@@ -50,8 +50,16 @@ void Particle::applyMagneticForce(const Vector3D& b, double dt) {
 
 double Particle::getMass() const {return mass;}
 
+void Particle::setMass(double m) {
+	mass = m;
+}
+
 double Particle::getCharge() const {return charge;}
 
+void Particle::setCharge(double q) {
+	charge = q;
+}
+
 
 Vector3D Particle::getVelocity() const {return velocity;}
 
diff --git a/src/main/Particle.h b/src/main/Particle.h
index 7c16176..43e059d 100644
--- a/src/main/Particle.h
+++ b/src/main/Particle.h
@@ -88,9 +88,15 @@ public:
 	/** Retourne la masse de cette particule. [kg] */
 	double getMass() const;
 
+	/** Affecte la masse. */
+	void setMass(double m);
+
 	/** Retourne la charge de cette particule. [C] */
 	double getCharge() const;
 
+	/** Affecte la charge */
+	void setCharge(double q);
+
 	/** Retourne la vitesse de cette particule. [m/s]*/
 	Vector3D getVelocity() const;
 
diff --git a/src/main/random.cc b/src/main/random.cc
new file mode 100644
index 0000000..58987fd
--- /dev/null
+++ b/src/main/random.cc
@@ -0,0 +1,49 @@
+/*
+ * random.cc
+ *
+ *  Created on: May 23, 2011
+ *      Author: jakob
+ */
+#include <math.h>
+#include "random.h"
+
+namespace vhc {
+
+namespace random {
+
+void seed(unsigned int s) {
+    if (s == 0) {
+    	srand(time(0));
+    } else {
+    	srand(s);
+    }
+}
+
+double uniform(double a, double b) {
+    return a + (rand() / double(RAND_MAX)) * (b - a);
+}
+
+double gaussian(double average, double standardDeviation) {
+  double v1;
+  double v2;
+  double s ;
+
+  do {
+    v1 = uniform(-1.0, 1.0);
+    v2 = uniform(-1.0, 1.0);
+    s = v1*v1 + v2*v2;
+  } while ((s >= 1.0) or (s == 0.0));
+
+  double x(sqrt(-2.0 * log(s) / s));
+  if (uniform(0.0, 1.0) > 0.5)
+    x *= v1;
+  else
+    x *= v2;
+
+  return average + standardDeviation * x;
+}
+
+} //random
+
+} //vhc
+
diff --git a/src/main/random.h b/src/main/random.h
new file mode 100644
index 0000000..6ba6264
--- /dev/null
+++ b/src/main/random.h
@@ -0,0 +1,26 @@
+/*
+ * random.h
+ *
+ *  Created on: May 23, 2011
+ *      Author: jakob
+ */
+
+#ifndef RANDOM_H_
+#define RANDOM_H_
+#include <cstdlib>
+#include <time.h>
+
+namespace vhc {
+
+namespace random {
+
+void seed(unsigned int s);
+
+double uniform(double a, double b);
+
+double gaussian(double average, double standardDeviation);
+
+} //random
+
+} //vhc
+#endif /* RANDOM_H_ */