From 64d6cd8a7ed46be4a34b4818705186f8e0998ab1 Mon Sep 17 00:00:00 2001 From: Jakob Odersky Date: Wed, 11 May 2011 14:59:00 +0000 Subject: Rajoute getHorizontalAt() pour tous les elements, en vue des faisceaux --- src/main/Accelerator.cc | 4 +++- src/main/CompositeElement.cc | 6 ++++++ src/main/CompositeElement.h | 2 ++ src/main/CurvedElement.cc | 5 +++++ src/main/CurvedElement.h | 2 ++ src/main/Element.h | 2 ++ src/main/Makefile | 2 +- src/main/Particle.cc | 3 ++- src/main/Printable.h | 3 +++ src/main/StraightElement.cc | 12 ++++++++---- src/main/StraightElement.h | 2 ++ src/main/Vector3D.cc | 2 +- src/main/constants.h | 2 +- 13 files changed, 38 insertions(+), 9 deletions(-) (limited to 'src/main') 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 #include 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 p à output. */ 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; -- cgit v1.2.3