PlasmaPy · pheuer · Mar 13, 2025 · Mar 13, 2025 · Mar 13, 2025 · Mar 13, 2025
diff --git a/changelog/2968.feature.rst b/changelog/2968.feature.rst
@@ -0,0 +1,4 @@
+Added the ``detector_vdir`` keyword to `~plasmapy.diagnostics.charged_particle_radiography.synthetic_radiography.Tracker`
+to explicitly define the detector's vertical surface vector. Added the ``source_vdir`` keyword to the ``create_particles``
+method of `~plasmapy.diagnostics.charged_particle_radiography.synthetic_radiography.Tracker` to explicitly define the
+orientation of the mean velocity of the source particles.
@@ -175,13 +175,17 @@ class Tracker(ParticleTracker):
         A unit vector (in Cartesian coordinates) defining the horizontal
         direction on the detector plane. By default, the horizontal axis in the
         detector plane is defined to be perpendicular to both the
-        source-to-detector vector and the z-axis (unless the source-to-detector
-        axis is parallel to the z axis, in which case the horizontal axis is
-        the x-axis).
+        origin-to-detector vector  (such that the detector is 'looking at' the origin)
+        and the z-axis (unless the origin-to-detector axis is parallel to the z axis,
+        in which case the horizontal axis is the x-axis).
+
+    detector_vdir : `numpy.ndarray`, shape (3), optional
+        A unit vector (in Cartesian coordinates) defining the vertical
+        direction on the detector plane. By default, the vertical axis in the
+        detector plane is defined to be perpendicular to both the
+        origin-to-detector vector (such that the detector is 'looking at' the origin)
+        and the detector horizontal axis.
 
-        The detector vertical axis is then defined
-        to be orthogonal to both the source-to-detector vector and the
-        detector horizontal axis.
 
     output_directory : `~pathlib.Path`, optional
         Directory for objects that are saved to disk. If a directory is not
@@ -211,6 +215,7 @@ def __init__(
             "volume averaged", "nearest neighbor"
         ] = "volume averaged",
         detector_hdir=None,
+        detector_vdir=None,
         output_directory: Path | None = None,
         output_basename: str = "output",
         fraction_exited_threshold: float = 0.999,
@@ -293,8 +298,11 @@ def __init__(
             self.det_hdir = self._default_detector_hdir()
 
         # Calculate the detector vdir
-        ny = np.cross(self.det_hdir, self.det_n)
-        self.det_vdir = -ny / np.linalg.norm(ny)
+        if detector_vdir is not None:
+            self.det_vdir = detector_vdir / np.linalg.norm(detector_vdir)
+        else:
+            ny = np.cross(self.det_hdir, self.det_n)
+            self.det_vdir = -ny / np.linalg.norm(ny)
 
     def _default_detector_hdir(self):
         """
@@ -609,6 +617,7 @@ def create_particles(
         max_theta=None,
         particle: Particle = Particle("p+"),  # noqa: B008
         distribution: Literal["monte-carlo", "uniform"] = "monte-carlo",
+        source_vdir=None,
         random_seed=None,
     ) -> None:
         r"""
@@ -661,6 +670,11 @@ def create_particles(
             on the image, but will well-sample the field grid with a
             smaller number of particles. The default is ``'monte-carlo'``.
 
+        source_vdir : (3,) |array_like|, default: None
+            A unit vector (in Cartesian coordinates) defining the orientation
+            of the mean of the particle velocities.  By default, the particle
+            velocities will be distributed around the source-detector axis.
+
         random_seed : int, optional
             A random seed to be used when generating random particle
             distributions, e.g. with the ``monte-carlo`` distribution.
@@ -670,6 +684,9 @@ def create_particles(
         # Raise an error if the run method has already been called.
         self._enforce_order()
 
+        if source_vdir is None:
+            source_vdir = self.src_det / np.linalg.norm(self.src_det)
+
         # Load inputs
         num_particles = int(num_particles)
 
@@ -705,10 +722,8 @@ def create_particles(
         v[:, 2] = v0 * np.cos(theta)
 
         # Calculate the rotation matrix that rotates the z-axis
-        # onto the source-detector axis
-        a = np.array([0, 0, 1])
-        b = self.detector - self.source
-        rot = rot_a_to_b(a, b)
+        # onto the vdir vector
+        rot = rot_a_to_b(np.array([0, 0, 1]), source_vdir)
 
         # Apply rotation matrix to calculated velocity distribution
         v = np.matmul(v, rot)

@@ -392,8 +392,11 @@ def test_init() -> None:
     sim = cpr.Tracker(grid, source, detector, verbose=False)
 
     # Test manually setting hdir and vdir
-    hdir = np.array([1, 0, 1])
-    sim = cpr.Tracker(grid, source, detector, verbose=False, detector_hdir=hdir)
+    hdir = np.array([1, 0, 0])
+    vdir = np.array([0, 0, 1])
+    sim = cpr.Tracker(
+        grid, source, detector, verbose=False, detector_hdir=hdir, detector_vdir=vdir
+    )
 
     # Test special case hdir == [0,0,1]
     source = (0 * u.mm, 0 * u.mm, -10 * u.mm)
@@ -433,6 +436,17 @@ def test_create_particles() -> None:
     # Test specifying particle
     sim.create_particles(1e3, 15 * u.MeV, particle="e-", random_seed=42)
 
+    # Test specifying direction
+    src_vdir = np.array([0.1, 1, 0])
+    src_vdir /= np.linalg.norm(src_vdir)
+    sim.create_particles(
+        1e3, 15 * u.MeV, particle="p+", random_seed=42, source_vdir=src_vdir
+    )
+    # Assert particle velocities are actually in that direction
+    vdir = np.mean(sim.v, axis=0)
+    vdir /= np.linalg.norm(vdir)
+    assert np.allclose(vdir, src_vdir, atol=0.05)
+
 
 @pytest.mark.slow
 def test_load_particles() -> None: