From a91ecfe471c4659e65be573e0217c8b72d5fa958 Mon Sep 17 00:00:00 2001 From: Jakob Odersky Date: Sun, 24 Apr 2011 13:31:11 +0000 Subject: *Ajout d'un pointeur sur les element precedents. *Implementation de la methode 'close()' dans l'accelerateur. *Restructuration des methodes geometriques des elements. -Ajout d'une methode isBefore(): determine si un point est contenu dans l'espace avant le plan de la face d'entree. -Ajout d'une methode isBeside(): determine si un point est plus proche de la trajectoire ideal que le rayon de section. -Ajout d'une methode isAfter(): determine si un point est contenu dans l'espace apres le plan de la face de sortie. -Ajout d'une methode contains(): determine si un point est contenu dans l'element (contains = !(isBefore || isBeside || isAfter)) -La methode isPast() (passe_au_suivant selon projet) a ete supprimee, elle est a remplacer par isAfter(). -La methode hitWall() (heurte_bord selon projet) a ete supprimee, elle est a remplacer par isBeside(). *Resolution d'un bug dans la methode isBeside() de StraightElement, avant une particule ne pouvait jamais heuter le bord! *Ajout des cibles gui et gui-build dans eclipse. --- src/gui/ElementRenderer.cc | 2 +- src/gui/Main.cc | 12 +++++++++--- src/gui/ParticleRenderer.cc | 3 +-- src/gui/Stage.cc | 2 +- src/gui/util.cc | 2 +- src/main/Accelerator.h | 43 +++++++++++++++++++++++++++++++++---------- src/main/CompositeElement.cc | 13 ++++++++----- src/main/CompositeElement.h | 8 ++++---- src/main/CurvedElement.cc | 16 ++++++++++------ src/main/CurvedElement.h | 11 +++++------ src/main/Element.cc | 24 +++++++++++++++++++++--- src/main/Element.h | 27 ++++++++++++++++++++++----- src/main/Quadrupole.cc | 6 ++---- src/main/StraightElement.cc | 18 +++++++++++------- src/main/StraightElement.h | 6 ++++-- src/main/Vector3D.cc | 5 +++++ src/main/Vector3D.h | 10 ++++++++-- 17 files changed, 146 insertions(+), 62 deletions(-) (limited to 'src') diff --git a/src/gui/ElementRenderer.cc b/src/gui/ElementRenderer.cc index 5f99491..a547851 100644 --- a/src/gui/ElementRenderer.cc +++ b/src/gui/ElementRenderer.cc @@ -60,7 +60,7 @@ void ElementRenderer::visit(const Quadrupole* quadrupole) const { void ElementRenderer::visit(const Dipole* dipole) const { - glColor4d(0.2, 0, 0.4, 0.9); + glColor4d(0.8, 0.4, 0, 0.9); glPushMatrix(); glTranslated(dipole->getCurvatureCenter().getX(), dipole->getCurvatureCenter().getY(), dipole->getCurvatureCenter().getZ()); Vector3D d = dipole->getExitPosition() - dipole->getCurvatureCenter(); diff --git a/src/gui/Main.cc b/src/gui/Main.cc index 03ba472..09465c2 100644 --- a/src/gui/Main.cc +++ b/src/gui/Main.cc @@ -9,6 +9,7 @@ #include #include #include +#include "exceptions.h" #include "Stage.h" #include "Accelerator.h" #include "StraightElement.h" @@ -112,14 +113,16 @@ Accelerator* singleDipole() { int main(int argc, char *argv[]) { + try { + QApplication app(argc, argv); vhc::Stage window; - Accelerator* acc = singleDipole(); + Accelerator* acc = linear(); window.accelerator = acc; - window.showFullScreen(); + //window.showFullScreen(); - //window.resize(QSize(500, 500)); + window.resize(QSize(500, 500)); window.setWindowTitle("Virtual Hadron Collider"); window.show(); @@ -127,4 +130,7 @@ int main(int argc, char *argv[]) delete acc; acc = NULL; + } catch (Exception& ex){ + std::cerr << ex.toString() << "\n"; + } } diff --git a/src/gui/ParticleRenderer.cc b/src/gui/ParticleRenderer.cc index 38736c7..3c8fe86 100644 --- a/src/gui/ParticleRenderer.cc +++ b/src/gui/ParticleRenderer.cc @@ -12,12 +12,11 @@ namespace vhc { ParticleRenderer::ParticleRenderer() { - //Auto-generated constructor stub } ParticleRenderer::~ParticleRenderer() { - //Auto-generated destructor stub + } void ParticleRenderer::render(const Particle& particle) const { diff --git a/src/gui/Stage.cc b/src/gui/Stage.cc index 6928b29..6dec091 100644 --- a/src/gui/Stage.cc +++ b/src/gui/Stage.cc @@ -137,7 +137,7 @@ void Stage::mouseMoveEvent(QMouseEvent* event) { int dpitch = -event->y() + center.y(); camera.addHeading(1.0 * dheading * frameTime / 4000); camera.addPitch(1.0 * dpitch * frameTime / 4000); - QCursor::setPos(center); + QCursor::setPos(QWidget::mapToGlobal(center)); update(); } diff --git a/src/gui/util.cc b/src/gui/util.cc index e460e36..5f026bf 100644 --- a/src/gui/util.cc +++ b/src/gui/util.cc @@ -24,7 +24,7 @@ void torus(double R, double r, double fraction, int slices, int stacks) { x = (R+r*cos(s*twopi/slices))*cos(t*twopi/stacks); y = (R+r*cos(s*twopi/slices))*sin(t*twopi/stacks); z = r * sin(s * twopi / slices); - glColor3d(0, 1 - 1.0 * j / (stacks * fraction), 1.0 * j / (stacks * fraction)); + // glColor3d(0, 1 - 1.0 * j / (stacks * fraction), 1.0 * j / (stacks * fraction)); glVertex3d(x, y, z); } } diff --git a/src/main/Accelerator.h b/src/main/Accelerator.h index a8d1562..3275a1e 100644 --- a/src/main/Accelerator.h +++ b/src/main/Accelerator.h @@ -9,6 +9,7 @@ #define ACCELERATOR_H_ #include #include +#include #include "exceptions.h" #include "Vector3D.h" #include "Particle.h" @@ -29,9 +30,13 @@ private : protected: - /** Attributs d'un accélérateur : une collection d'éléments, et une de particules. */ - std::vector< Element* > elementCollec; - std::vector< Particle* > particleCollec; + /*Attributs d'un accélérateur : une collection d'éléments, et une de particules. */ + + /** Collection d'elements contenus dans cet accelerateur. */ + std::vector< Element * > elementCollec; + + /** Collection de particules contenus danc cet accelerateur. */ + std::vector< Particle *> particleCollec; public: @@ -44,28 +49,45 @@ public: virtual ~Accelerator() {clear();} /** Retourne un pointeur sur un élément de l'accélérateur.*/ - Element* getElement(int rank) const { return elementCollec[rank]; } + Element *const getElement(int rank) const { return elementCollec[rank];} /* Retourne un pointeur sur une particule de l'accélérateur.*/ - Particle* getParticle(int rank) const { return particleCollec[rank];} + Particle *const getParticle(int rank) const { return particleCollec[rank];} /** Retourne une représentation en chaîne de caractères de cet accélérateur. */ virtual std::string toString() const; /** Copie un élément dans l'accélérateur. */ - Element* add(const Element& element) { + Element *const add(const Element& element) { Element* e = element.clone(); elementCollec.push_back(e); return e; } /** Copie une particule dans l'accélérateur. */ - Particle* add(const Particle& particle){ + Particle *const add(const Particle& particle){ Particle* p = particle.clone(); particleCollec.push_back(p); return p; } + void close() { + for (int i = 0; i < elementCollec.size(); ++i) { + //element i + Element* current = elementCollec[i]; + + //element i+1 + Element* next = elementCollec[(i+1) % elementCollec.size()]; + + // est-ce que les elements se suivent? + 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)"); + } + + } + /** Efface tous les éléments et les particules. */ void clear() { for (unsigned int i = 0; i < particleCollec.size(); ++i) { @@ -119,8 +141,8 @@ public: particle.setPosition(particle.getPosition() + particle.getVelocity() * dt); particle.setForce(Vector3D::Null); - if (particle.getElement()->hasHit(particle)) std::cout << "Particle hit wall!" << std::endl; - if (particle.getElement()->isPast(particle)) { // si la particule est passee, qui sait si elle est dans l'element suivant? + if (particle.getElement()->isBeside(particle)) std::cout << "Particle hit wall!" << std::endl; + if (particle.getElement()->isAfter(particle)) { // si la particule est passee, qui sait si elle est dans l'element suivant? if (!particle.getElement()->isConnected()) throw Exception("Element in accelerator not connected."); else particle.setElement(particle.getElement()->getNext()); } @@ -128,12 +150,13 @@ public: } + /* void* stepMonteCarlo(); void* stepDeterministBackUp(); void* stepDeterminist(); - +*/ }; } diff --git a/src/main/CompositeElement.cc b/src/main/CompositeElement.cc index 4c4170e..81f77aa 100644 --- a/src/main/CompositeElement.cc +++ b/src/main/CompositeElement.cc @@ -15,16 +15,19 @@ CompositeElement::CompositeElement(const Vector3D& entry, const Vector3D& exit, CompositeElement::~CompositeElement() {}; -bool CompositeElement::hasHit(const Particle& particle) const { +bool CompositeElement::isBefore(const Vector3D& position) const { + return elements.front()->isBefore(position); +} + +bool CompositeElement::isBeside(const Vector3D& position) const { for (int i(0); i < elements.size(); ++i) { - if (elements[i]->hasHit(particle)) return true; + if (elements[i]->isBeside(position)) return true; } return false; } -bool CompositeElement::isPast(const Particle& particle) const { - if (elements[elements.size() - 1]->isPast(particle)) return true; - else return false; +bool CompositeElement::isAfter(const Vector3D& position) const { + return elements.back()->isAfter(position); } Vector3D CompositeElement::magneticFieldAt(const Vector3D& position) const { diff --git a/src/main/CompositeElement.h b/src/main/CompositeElement.h index 4cf55a8..c8d2286 100644 --- a/src/main/CompositeElement.h +++ b/src/main/CompositeElement.h @@ -32,11 +32,11 @@ public: CompositeElement(const Vector3D& entry, const Vector3D& exit, double sectionRadius, Element* next = NULL); virtual ~CompositeElement(); - /** Determine si une particule a heurte le bord de cet element, donc d'un des elements composants. */ - virtual bool hasHit(const Particle& particle) const; + virtual bool isBefore(const Vector3D& position) const; - /** Determine si une particule a passe cet element, donc passe le dernier element composant. */ - virtual bool isPast(const Particle& particle) const; + virtual bool isBeside(const Vector3D& position) const; + + virtual bool isAfter(const Vector3D& position) const; virtual Vector3D magneticFieldAt(const Vector3D& position) const; diff --git a/src/main/CurvedElement.cc b/src/main/CurvedElement.cc index ca1a4bd..30efebb 100644 --- a/src/main/CurvedElement.cc +++ b/src/main/CurvedElement.cc @@ -13,7 +13,6 @@ using namespace std; namespace vhc { - /** Cf. CurvedElement.h */ CurvedElement::CurvedElement(const Vector3D& entry, const Vector3D& exit, double sectionRadius, double curvature, Element* next): Element(entry, exit, sectionRadius, next), curvature(curvature), @@ -44,15 +43,20 @@ double CurvedElement::getAngle() const { return acos((entryPosition - curvatureCenter).unit().dot((exitPosition - curvatureCenter).unit())); } -bool CurvedElement::hasHit(const Particle& particle) const { - Vector3D X(particle.getPosition() - curvatureCenter); - Vector3D u = (X - particle.getPosition().getZ() * Vector3D::k).unit(); +bool CurvedElement::isBefore(const Vector3D& position) const { + Vector3D out = (exitPosition - curvatureCenter).cross(entryPosition - curvatureCenter).cross(entryPosition - curvatureCenter); + return (position - entryPosition).dot(out) > 0; +} + +bool CurvedElement::isBeside(const Vector3D& position) const { + Vector3D X(position - curvatureCenter); + Vector3D u = (X - position.getZ() * Vector3D::k).unit(); return (X - fabs(1.0 / curvature) * u).norm() > sectionRadius; } -bool CurvedElement::isPast(const Particle& particle) const { +bool CurvedElement::isAfter(const Vector3D& position) const { Vector3D out = (entryPosition - curvatureCenter).cross(exitPosition - curvatureCenter).cross(exitPosition - curvatureCenter); - return (particle.getPosition() - exitPosition).dot(out) > 0; + return (position - exitPosition).dot(out) > 0; } std::string CurvedElement::getType() const {return "Curved Element";} diff --git a/src/main/CurvedElement.h b/src/main/CurvedElement.h index f7cc902..ac81746 100644 --- a/src/main/CurvedElement.h +++ b/src/main/CurvedElement.h @@ -45,15 +45,14 @@ public: * @param next pointeur sur l'element suivant */ CurvedElement(const Vector3D& entry, const Vector3D& exit, double sectionRadius, double curvature, Element* next = NULL); - + /** Destructeur virtuel. */ virtual ~CurvedElement(); - //virtual CurvedElement* clone() const {return new CurvedElement(*this);} - //TODO !!! erreurs dans l'algorithme - virtual bool hasHit(const Particle& particle) const; + virtual bool isBefore(const Vector3D& position) const; + + virtual bool isBeside(const Vector3D& position) const; - //TODO !!! erreurs dans l'algorithme - virtual bool isPast(const Particle& particle) const; + virtual bool isAfter(const Vector3D& position) const; /** Retourne la courbure. */ double getCurvature() const; diff --git a/src/main/Element.cc b/src/main/Element.cc index ed66d94..25efa8c 100644 --- a/src/main/Element.cc +++ b/src/main/Element.cc @@ -15,11 +15,25 @@ Element::Element(const Vector3D& entry, const Vector3D& exit, double sectionRadi entryPosition(entry), exitPosition(exit), sectionRadius(sectionRadius), - next(next) { + previous(NULL), + next(next) + {} +Element::~Element() {}; + +bool Element::isAfter(const Particle& particle) const {return isAfter(particle.getPosition());} + +bool Element::isBefore(const Particle& particle) const {return isBefore(particle.getPosition());} + +bool Element::isBeside(const Particle& particle) const {return isBeside(particle.getPosition());} + +bool Element::contains(const Vector3D& position) const { + return !(isAfter(position) || isBefore(position) || isBeside(position)); } -Element::~Element() {}; +bool Element::contains(const Particle& particle) const { + return contains(particle.getPosition()); +} Vector3D Element::magneticFieldAt(const Vector3D& position) const {return Vector3D::Null;} @@ -33,7 +47,11 @@ Vector3D Element::getExitPosition() const {return exitPosition;} double Element::getSectionRadius() const {return sectionRadius;} -Element* Element::getNext() const {return next;} +Element* const Element::getPrevious() const {return previous;} + +void Element::setPrevious(Element* p) {previous = p;} + +Element* const Element::getNext() const {return next;} void Element::setNext(Element* n) {next = n;} diff --git a/src/main/Element.h b/src/main/Element.h index a34200c..28fc9ed 100644 --- a/src/main/Element.h +++ b/src/main/Element.h @@ -38,6 +38,10 @@ protected: /** Rayon de la chambre à vide. */ double sectionRadius; + /** Pointeur sur l'élément precedent. + * NULL si aucun element n'existe. */ + Element *previous; + /** Pointeur sur l'élément suivant. * NULL si aucun element n'existe. */ Element *next; @@ -58,11 +62,17 @@ public: * ATTENTION: La delocation de memoire est sous la responsabilite de l'appelant. */ virtual Element* clone() const = 0; - /** Determine si la particule donnee a heurte le bord de cet element. */ - virtual bool hasHit(const Particle& particle) const = 0; + virtual bool isBefore(const Vector3D& position) const = 0; + bool isBefore(const Particle& particle) const; + + virtual bool isBeside(const Vector3D& position) const = 0; + bool isBeside(const Particle& particle) const; + + virtual bool isAfter(const Vector3D& position) const = 0; + bool isAfter(const Particle& particle) const; - /** Determine si la particule donnee a passe cet element. */ - virtual bool isPast(const Particle& particle) const = 0; + bool contains(const Vector3D& position) const; + bool contains(const Particle& particle) const; /** Retourne le champ magnetique, a l'interieur de cet element a la position donnee. */ virtual Vector3D magneticFieldAt(const Vector3D& position) const; @@ -93,7 +103,13 @@ public: //void setSectionRadius(double r) {sectionRadius = r;} /** Retourne un pointeur l'element suivant. NULL s'il n'existe pas. */ - Element* getNext() const; + Element* const getPrevious() const; + + /** Assigne un pointeur sur l'element suivant. */ + void setPrevious(Element* p); + + /** Retourne un pointeur sur l'element suivant. NULL s'il n'existe pas. */ + Element* const getNext() const; /** Assigne un pointeur sur l'element suivant. */ void setNext(Element* n); @@ -111,6 +127,7 @@ public: //TODO expl. virtual void accept(const ElementVisitor& v) const = 0; + }; } diff --git a/src/main/Quadrupole.cc b/src/main/Quadrupole.cc index 911ae18..ff0e952 100644 --- a/src/main/Quadrupole.cc +++ b/src/main/Quadrupole.cc @@ -19,14 +19,12 @@ Quadrupole::Quadrupole(const Vector3D& entry, const Vector3D& exit, double secti Quadrupole::~Quadrupole() {}; Vector3D Quadrupole::magneticFieldAt(const Vector3D& position) const { + if (!Quadrupole::contains(position)) return Vector3D::Null; Vector3D x = position - getEntryPosition(); Vector3D d = getDiagonal().unit(); Vector3D y = x - x.dot(d) * d; Vector3D u = Vector3D::k.cross(d); - if ((x.dot(d) >= 0) && ((position - getExitPosition()).dot(-d) >= 0)) - return focalizingCoefficient * (y.dot(u) * Vector3D::k + position.getZ() * u); - else - return Vector3D::Null; + return focalizingCoefficient * (y.dot(u) * Vector3D::k + position.getZ() * u); } double Quadrupole::getFocalizingCoefficient() const { diff --git a/src/main/StraightElement.cc b/src/main/StraightElement.cc index 1739486..24d63aa 100644 --- a/src/main/StraightElement.cc +++ b/src/main/StraightElement.cc @@ -15,16 +15,20 @@ 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(); +} -bool StraightElement::hasHit(const Particle& particle) const { - Vector3D a(particle.getPosition() - entryPosition); - const Vector3D b = (particle.getPosition() - entryPosition); - return (a.cross(b)).norm() / getDiagonal().norm() > sectionRadius; +bool StraightElement::isBeside(const Vector3D& position) const { + Vector3D X = Vector3D(position - entryPosition); + Vector3D d = getDiagonal().unit(); + return (X - X.dot(d) * d).norm() > sectionRadius; }; -bool StraightElement::isPast(const Particle& particle) const { - const Vector3D v(particle.getPosition() - entryPosition); - return getDiagonal().dot(v) > getDiagonal().dot(getDiagonal()); +bool StraightElement::isAfter(const Vector3D& position) const { + const Vector3D v(position - entryPosition); + return getDiagonal().dot(v) > getDiagonal().normSquare(); } std::string StraightElement::getType() const {return "Straight Element";} diff --git a/src/main/StraightElement.h b/src/main/StraightElement.h index 6de596a..9a3a366 100644 --- a/src/main/StraightElement.h +++ b/src/main/StraightElement.h @@ -24,9 +24,11 @@ public: virtual ~StraightElement(); - virtual bool hasHit(const Particle& particle) const; + virtual bool isBefore(const Vector3D& position) const; - virtual bool isPast(const Particle& particle) const; + virtual bool isBeside(const Vector3D& position) const; + + virtual bool isAfter(const Vector3D& position) const; virtual std::string getType() const; virtual std::string toString() const; diff --git a/src/main/Vector3D.cc b/src/main/Vector3D.cc index 44ccd47..94952c3 100644 --- a/src/main/Vector3D.cc +++ b/src/main/Vector3D.cc @@ -5,6 +5,7 @@ * Author: Jakob Odersky */ +#include #include #include "Vector3D.h" #include "exceptions.h" @@ -51,6 +52,8 @@ double Vector3D::norm() const { return normCache; } +double Vector3D::minNorm() const {return min(x, min(y, z));} + double Vector3D::normSquare() const {return dot(*this);} std::string Vector3D::toString() const { @@ -69,6 +72,8 @@ Vector3D Vector3D::rotate(const Vector3D& axis, double t) const { return x * cos(t) + a * x.dot(a) * (1-cos(t)) + a.cross(x) * sin(t); } +bool Vector3D::ae(const Vector3D& u, const Vector3D& v, double epsilon) {return (u - v).norm() <= epsilon;} + //c.f. header Vector3D const Vector3D::Null = Vector3D(0.0, 0.0, 0.0); diff --git a/src/main/Vector3D.h b/src/main/Vector3D.h index 3c41147..ee7ffcc 100644 --- a/src/main/Vector3D.h +++ b/src/main/Vector3D.h @@ -8,7 +8,7 @@ #ifndef VECTOR3D_H_ #define VECTOR3D_H_ -#include +#include #include "Printable.h" namespace vhc { @@ -93,9 +93,12 @@ public: /** Retourne le vecteur unitaire */ Vector3D unit() const; - /** Retourne la norme du vecteur. */ + /** Retourne la norme euclidienne du vecteur. */ double norm() const; + /** Retourne la norme de ce vecteur correspondant au minimum des composantes. */ + double minNorm() const; + /** Retourne la norme du vecteur au carre. */ double normSquare() const; @@ -112,6 +115,9 @@ public: * le vecteur de l'axe, et en t, l'angle de rotation. */ Vector3D rotate(const Vector3D& axis, double t) const; + /** Verifie si les vecteurs u et v sont approximativement égaux, i.e. que la norme de leur difference est plus petit qu'un epsilon. */ + static bool ae(const Vector3D& u, const Vector3D& v, double epsilon = std::numeric_limits::epsilon()); + /** Vecteur nul. (0,0,0) */ static const Vector3D Null; -- cgit v1.2.3