Merge pull request #512 from ehlertdo/ReflectiveShell

Implemented ReflectiveShell as optional boundary

Author: JulienDoerner

CHANGELOG.md

@@ -7,12 +7,15 @@
 * Fixed wrong generation of interval ranges in ObserverTimeEvolution
 
 ### New features:
-* Added new backwards-compatible function particleMass that returns particle mass also for non-nuclei
-* Added the new Galactic magnetic field models from Unger&Farrar arXiv:2311.12120 and Korochkin et al. arXiv:2407.02148
-* Added EBL model from Finke et al. 2022
-* Added overridable getTime function to ObserverTimeEvolution which is called instead of detList. 
+
+ * Added new backwards-compatible function particleMass that returns particle mass also for non-nuclei
+ * Added the new Galactic magnetic field models from Unger&Farrar arXiv:2311.12120
+ * Added EBL model from Finke et al. 2022
+ * Added ReflectiveShell boundary condition
+ * Added overridable getTime function to ObserverTimeEvolution which is called instead of detList. 
   Old functionalities are preserved by adding old functions and setter/getter functions.
 
+
 ### Interface changes:
 
 ### Features that are deprecated and will be removed after this release

doc/pages/Simulation-Modules.md

@@ -41,7 +41,7 @@ General interactions/processes
 Conditional modules implement certain conditions for stopping propagation.
 They provide interfaces to act onReject or onAccept of a cosmic ray.
 Boundary modules can be used to limit the simulation volume.
-Periodic- and ReflectiveBox implement boundary conditions for the particles. They are useful for a 3D setup where an initial volume is to be repeated (periodically or reflectively). When using them, a couple of things need to be considered. Observers will shadow the volume behind if they are set to inactivate particles. Also Observers should be placed at a distance to the boundaries that is larger than the maximum step size of the propagator, since step size limitation does not work beyond periodic/reflective boundaries.
+Periodic-/ReflectiveBox and ReflectiveShell implement boundary conditions for the particles. They are useful for a 3D setup where an initial volume is to be repeated (periodically or reflectively). When using them, a couple of things need to be considered. Observers will shadow the volume behind if they are set to inactivate particles. Also Observers should be placed at a distance to the boundaries that is larger than the maximum step size of the propagator, since step size limitation does not work beyond periodic/reflective boundaries.
 
 * **MaximumTrajectoryLength** - Stop after reaching maximum trajectory length
 * **MinimumEnergy** - Stop after reaching a minimum energy
@@ -53,6 +53,7 @@ Periodic- and ReflectiveBox implement boundary conditions for the particles. The
 * **CylindricalBoundary** - Cylindric simulation volume
 * **PeriodicBox** - Periodic boundary conditions for the particle: If a particle leaves the box it will enter from the opposite side and the initial position will be changed as if it had come from that side.
 * **ReflectiveBox** - Reflective boundary conditions for the particle: If a particle leaves the box it will be reflected (mirrored) and the initial position will be changed as if it had come from that side.
+* **ReflectiveShell** - Reflective spherical boundary condition for the particle; If a particle crosses the shell during the current step it will be reflected back at the boundary according to the law of reflection. Particles coming from the outside are only reflected if the current step ends inside the shell.
 * **DetectionLength** - Detects the candidate at a given trajectory length.
 
 ### Observers

include/crpropa/module/Boundary.h

@@ -37,6 +37,32 @@ class PeriodicBox: public Module {
 	std::string getDescription() const;
 };
 
+/**
+ @class ReflectiveShell
+ @brief Shell with reflective boundary
+
+ If a particle crosses the boundary it is reflected back inside or outside (position and velocity) depending where it came from.
+ Particles can overshoot (be outside of the box during the step) since the step size is not limited by this module.
+ */
+class ReflectiveShell: public Module {
+private:
+	Vector3d center;
+	double radius;
+
+public:
+	/** Constructor
+	 @param center	vector corresponding to the center of the sphere
+	 @param r		value corresponding to the radius of the shell
+	 */
+	ReflectiveShell(Vector3d center, double r);
+	double distance(const Vector3d &point) const;
+	Vector3d normal(const Vector3d &point) const;
+	void process(Candidate *candidate) const;
+	void setCenter(Vector3d center);
+	void setRadius(double r);
+	std::string getDescription() const;
+};
+
 /**
  @class ReflectiveBox
  @brief Rectangular box with reflective boundaries.

src/module/Boundary.cpp

@@ -40,6 +40,57 @@ std::string PeriodicBox::getDescription() const {
 	return s.str();
 }
 
+
+ReflectiveShell::ReflectiveShell(Vector3d center, double r) :
+		center(center), radius(r) {
+		}
+
+double ReflectiveShell::distance(const Vector3d &point) const {
+	Vector3d dR = point - center;
+	return dR.getR() - radius;
+}
+
+Vector3d ReflectiveShell::normal(const Vector3d& point) const {
+	Vector3d d = point - center;
+	return d.getUnitVector();
+}
+
+void ReflectiveShell::process(Candidate *c) const {
+	double currentDistance = distance(c->current.getPosition());
+	double previousDistance = distance(c->previous.getPosition());
+	// check if cosmic ray crossed boundary in last step
+	if (currentDistance * previousDistance < 0){
+		Vector3d currentDirection = c->current.getDirection();
+		Vector3d previousPosition = c->previous.getPosition();
+		// get point where trajectory intersects the shell boundary
+		double p_half = previousPosition.dot(currentDirection) / currentDirection.getR2();
+		double q = (previousPosition.getR2() - radius * radius) / currentDirection.getR2();
+		double k1 = - p_half + sqrt(p_half * p_half - q);
+		Vector3d intersectPoint = previousPosition + k1 * currentDirection;
+		// flip component of velocity normal to surface
+		Vector3d surfaceNormal = normal(intersectPoint);
+		Vector3d newDirection = currentDirection - 2 * currentDirection.dot(surfaceNormal) * surfaceNormal;
+		c->current.setDirection(newDirection.getUnitVector());
+		// also update position
+		c->current.setPosition(intersectPoint + newDirection * (c->getCurrentStep() - (k1 * currentDirection).getR()) / newDirection.getR());
+	}
+}
+
+void ReflectiveShell::setCenter(Vector3d c) {
+	center = c;
+}
+void ReflectiveShell::setRadius(double r) {
+	radius = r;
+}
+
+std::string ReflectiveShell::getDescription() const {
+	std::stringstream s;
+	s << "Reflective sphere: center: " << center / Mpc << " Mpc, ";
+	s << "radius: " << radius / Mpc << " Mpc";
+	return s.str();
+};
+
+
 ReflectiveBox::ReflectiveBox() :
 		origin(Vector3d(0, 0, 0)), size(Vector3d(0, 0, 0)) {
 }

test/testBreakCondition.cpp

@@ -413,6 +413,34 @@ TEST(PeriodicBox, low) {
 	EXPECT_DOUBLE_EQ(3, c.created.getPosition().z);
 }
 
+TEST(ReflectiveShell, inside) {
+	// Tests if the reflective boundaries place the particle back inside the shell
+	Vector3d center(0, 0, 0);
+	double radius = 100;
+	ReflectiveShell shell(center, radius);
+
+	Candidate c;
+	c.setCurrentStep(20);
+	c.previous.setPosition(Vector3d(80, 20, 30));
+	c.previous.setDirection(Vector3d(10, -1, -1));
+	// un-reflected new position (outside shell) after full step
+	Vector3d currentPosition = c.previous.getPosition() + c.previous.getDirection() * c.getCurrentStep() / c.previous.getDirection().getR();
+	c.current.setPosition(currentPosition);
+	c.current.setDirection(Vector3d(10, -1, -1));
+	// process reflection
+	shell.process(&c);
+
+	// expected position & direction after reflection
+	// calculated by hand with the same algorithm as implemented in Boundary.cpp
+	EXPECT_NEAR(90.06356247, c.current.getPosition().x, 1e-7);
+	EXPECT_NEAR(16.09256317, c.current.getPosition().y, 1e-7);
+	EXPECT_NEAR(25.05646296, c.current.getPosition().z, 1e-7);
+
+	EXPECT_NEAR(-0.67179589, c.current.getDirection().x, 1e-7);
+	EXPECT_NEAR(-0.42786503, c.current.getDirection().y, 1e-7);
+	EXPECT_NEAR(-0.60466668, c.current.getDirection().z, 1e-7);
+}
+
 TEST(ReflectiveBox, high) {
 	// Tests if the reflective boundaries place the particle back inside the box and translate the initial position accordingly.
 	// Also the initial and final directions are to be reflected