Repository: PlasmaPy/PlasmaPy
Pull request #1035 (2021-02-23 14:51:13)

PRad Bugfix
I found a few minor bugs in sections of `proton_radiography.SyntheticProtonRadiograph` that I used rarely: 

- Accessing the `max_deflection` property raised an exception if particles had been deflected away from the detector (since `v` and `v_init` were then un-equal in length). I fixed this by removing deflected particles from `v_init` as well in `_coast_to_detector`.

- `_adaptive_dt` depended on the velocity of the particles `v0`, which was calculated in `create_particles` but not defined if the user provided particles in `load_particles`. Furthermore, user-provided particles don't necessarily have the same energy. I solved that by using `max(v)` in `_adaptive_dt`, which I'm calculating once at the beginning of `run()`.

- `_adaptive_dt` now calculates the cyclotron period using the provided particle charge and mass (previously assumed protons).

- [x] I have added a changelog entry for this pull request.

codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035?src=pr&el=h1) Report
> Merging [#1035](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035?src=pr&el=desc) (baba2db) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/d4371b59a7e1906fd9e5a18568a629b0eee3e39d?el=desc) (d4371b5) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035/graphs/tree.svg?width=650&height=150&src=pr&token=LzZ5mnXVad)](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master    #1035   +/-   ##
=======================================
  Coverage   96.82%   96.82%           
=======================================
  Files          69       69           
  Lines        6737     6739    +2     
=======================================
+ Hits         6523     6525    +2     
  Misses        214      214           
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/diagnostics/proton\_radiography.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035/diff?src=pr&el=tree#diff-cGxhc21hcHkvZGlhZ25vc3RpY3MvcHJvdG9uX3JhZGlvZ3JhcGh5LnB5) | `99.65% <100.00%> (+<0.01%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035?src=pr&el=footer). Last update [d4371b5...baba2db](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1035?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/changelog/1035.bugfix.rst b/changelog/1035.bugfix.rst
new file mode 100644
index 00000000..58ddb5ec
--- /dev/null
+++ b/changelog/1035.bugfix.rst
@@ -0,0 +1,4 @@
+Fixed a few bugs in `proton_radiography.SyntheticProtonRadiograph`:
+- The `max_deflection` property now works properly even when particles are deflected from the grid.
+- The adaptive timestep function now works when particles are provided using `load_particles`
+- Adaptive dt now calculates the cyclotron period using the provided particle charge and mass (previously assumed protons).
diff --git a/plasmapy/diagnostics/proton_radiography.py b/plasmapy/diagnostics/proton_radiography.py
index ff165bf1..52e21358 100644
--- a/plasmapy/diagnostics/proton_radiography.py
+++ b/plasmapy/diagnostics/proton_radiography.py
@@ -304,9 +304,7 @@ def _log(self, msg):
             print(msg)
 
     # Define some constants so they don't get constantly re-evaluated
-    _e = const.e.si.value
     _c = const.c.si.value
-    _m_p = const.m_p.si.value
 
     # *************************************************************************
     # Particle creation methods
@@ -436,7 +434,7 @@ def create_particles(
 
         # Calculate the velocity corresponding to the proton energy
         ER = self.proton_energy * 1.6e-19 / (self.m * self._c ** 2)
-        self.v0 = self._c * np.sqrt(1 - 1 / (ER + 1) ** 2)
+        v0 = self._c * np.sqrt(1 - 1 / (ER + 1) ** 2)
 
         if distribution == "monte-carlo":
             theta, phi = self._angles_monte_carlo()
@@ -449,9 +447,9 @@ def create_particles(
 
         # Construct the velocity distribution around the z-axis
         self.v = np.zeros([self.nparticles, 3])
-        self.v[:, 0] = self.v0 * np.sin(theta) * np.cos(phi)
-        self.v[:, 1] = self.v0 * np.sin(theta) * np.sin(phi)
-        self.v[:, 2] = self.v0 * np.cos(theta)
+        self.v[:, 0] = v0 * np.sin(theta) * np.cos(phi)
+        self.v[:, 1] = v0 * np.sin(theta) * np.sin(phi)
+        self.v[:, 2] = v0 * np.cos(theta)
 
         # Calculate the rotation matrix that rotates the z-axis
         # onto the source-detector axis
@@ -543,7 +541,7 @@ def _adaptive_dt(self, Ex, Ey, Ez, Bx, By, Bz):
 
         # Compute the timestep indicated by the grid resolution
         ds = self.grid.grid_resolution.to(u.m).value
-        gridstep = 0.5 * (np.min(ds) / self.v0)
+        gridstep = 0.5 * (np.min(ds) / self.vmax)
 
         # If not, compute a number of possible timesteps
         # Compute the cyclotron gyroperiod
@@ -553,7 +551,7 @@ def _adaptive_dt(self, Ex, Ey, Ez, Bx, By, Bz):
         if Bmag == 0:
             gyroperiod = np.inf
         else:
-            gyroperiod = 2 * np.pi * self._m_p / (self._e * np.max(Bmag))
+            gyroperiod = 2 * np.pi * self.m / (self.q * np.max(Bmag))
 
         # TODO: introduce a minimum timestep based on electric fields too!
 
@@ -632,6 +630,7 @@ def _coast_to_detector(self):
         # Drop the other particles
         self.x = self.x[ind, :]
         self.v = self.v[ind, :]
+        self.v_init = self.v_init[ind, :]
         self.nparticles_grid = self.x.shape[0]
         t = t[ind]
 
@@ -807,6 +806,10 @@ def run(
         # This will be used later to calculate the maximum deflection
         self.v_init = np.copy(self.v)
 
+        # Calculate the maximum velocity
+        # Used for determining the grid crossing maximum timestep
+        self.vmax = np.max(np.linalg.norm(self.v, axis=-1))
+
         # Determine which particles should be tracked
         # This array holds the indices of all particles that WILL hit the grid
         # Only these particles will actually be pushed through the fields


Test patch:
diff --git a/plasmapy/diagnostics/tests/test_proton_radiography.py b/plasmapy/diagnostics/tests/test_proton_radiography.py
index c4116070..214f24a3 100644
--- a/plasmapy/diagnostics/tests/test_proton_radiography.py
+++ b/plasmapy/diagnostics/tests/test_proton_radiography.py
@@ -365,13 +365,18 @@ def test_run_options():
     # Test extreme deflections -> warns user
     # This requires instatiating a whole new example field with a really
     # big B-field
-    grid = _test_grid("constant_bz", num=50, B0=500 * u.T)
+    grid = _test_grid("constant_bz", num=50, B0=250 * u.T)
     source = (0 * u.mm, -10 * u.mm, 0 * u.mm)
     detector = (0 * u.mm, 200 * u.mm, 0 * u.mm)
     sim = prad.SyntheticProtonRadiograph(grid, source, detector, verbose=False)
-    sim.create_particles(1e4, 3 * u.MeV, max_theta=1 * u.deg)
+    sim.create_particles(1e4, 3 * u.MeV, max_theta=0.1 * u.deg)
     with pytest.warns(RuntimeWarning):
         sim.run(field_weighting="nearest neighbor", dt=1e-12 * u.s)
+    # Calc max deflection: should be between 0 and pi/2
+    # Note: that's only true because max_theta is very small
+    # More generally, max_deflection can be a bit bigger than pi/2 for
+    # particles that begin at an angle then deflect all the way around.
+    assert 0 < sim.max_deflection.to(u.rad).value < np.pi / 2
 
 
 def test_synthetic_radiograph():