From 4dfe6cc0636809542f525fdd86eaaa2af3181647 Mon Sep 17 00:00:00 2001 From: Jakob Odersky Date: Mon, 23 May 2011 16:19:06 +0000 Subject: Ajout des bunchs :) --- src/gui/Main.cc | 11 ++++- src/gui/Stage.cc | 4 +- src/main/Beam.h | 20 +++++--- src/main/Bunch.cc | 120 ++++++++++++++++++++++++++++++++++++++++++++++++ src/main/Bunch.h | 78 +++++++++++++++++++++++++++++++ src/main/CircularBeam.h | 1 + src/main/Makefile | 3 +- src/main/Particle.cc | 8 ++++ src/main/Particle.h | 6 +++ src/main/random.cc | 49 ++++++++++++++++++++ src/main/random.h | 26 +++++++++++ 11 files changed, 316 insertions(+), 10 deletions(-) create mode 100644 src/main/Bunch.cc create mode 100644 src/main/Bunch.h create mode 100644 src/main/random.cc create mode 100644 src/main/random.h (limited to 'src') 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 #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 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) * ATTENTION: 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 +#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 dt. */ + 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 +#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 +#include + +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_ */ -- cgit v1.2.3