From 40f89caa223f6f84bc3a9d84681909c4704bb614 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 8 Jul 2019 23:24:19 -0500
Subject: [PATCH 1/2] Fix cylinder_from_points function

---
 openmc/model/funcs.py            |  6 +++---
 tests/unit_tests/test_surface.py | 35 +++++++++++++++++++++++---------
 2 files changed, 28 insertions(+), 13 deletions(-)

diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py
index 143033eca8..4e79617dd8 100644
--- a/openmc/model/funcs.py
+++ b/openmc/model/funcs.py
@@ -395,9 +395,9 @@ def cylinder_from_points(p1, p2, r, **kwargs):
     dx = x2 - x1
     dy = y2 - y1
     dz = z2 - z1
-    cx = y1*z2 + y2*z1
-    cy = -(x1*z2 + x2*z1)
-    cz = x1*y2 + x2*y1
+    cx = y1*z2 - y2*z1
+    cy = x2*z1 - x1*z2
+    cz = x1*y2 - x2*y1
 
     # Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation for the
     # cylinder can be derived as r = |(p - p1) ⨯ (p - p2)| / |p2 - p1|.
diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py
index 6ed657a4c5..09bfa738e3 100644
--- a/tests/unit_tests/test_surface.py
+++ b/tests/unit_tests/test_surface.py
@@ -1,4 +1,5 @@
 import numpy as np
+from random import random
 import openmc
 import pytest
 
@@ -329,14 +330,28 @@ def test_quadric():
 
 
 def test_cylinder_from_points():
-    # Generate 45-degree rotated cylinder in x-y plane with radius 1
-    p1 = (0, 0, 0)
-    p2 = (1, 1, 0)
-    s = openmc.model.cylinder_from_points(p1, p2, 1)
+    for _ in range(10):
+        # Generate cylinder in random direction
+        p1 = np.array([random(), random(), random()])
+        p2 = np.array([random(), random(), random()])
+        r = random()
+        s = openmc.model.cylinder_from_points(p1, p2, r)
 
-    # Points p1 and p2 need to be inside cylinder
-    assert p1 in -s
-    assert p2 in -s
-    assert (-1, 1, 0) in +s
-    assert (1, -1, 0) in +s
-    assert (0, 0, 1.5) in +s
+        # Points p1 and p2 need to be inside cylinder
+        assert p1 in -s
+        assert p2 in -s
+
+        # Points further along the line should be inside cylinder as well
+        t = 100*random() - 200
+        p = p1 + t*(p2 - p1)
+        assert p in -s
+
+        # Check that a point outside cylinder is in positive half-space. We do
+        # this by constructing a plane that includes the cylinder's axis,
+        # finding the normal to the plane, and using it to find a point slightly
+        # more than one radius away from the axis.
+        plane = openmc.Plane.from_points(p1, p2, (0., 0., 0.))
+        n = np.array([plane.a, plane.b, plane.c])
+        n /= np.linalg.norm(n)
+        assert (p1 + 1.1*r*n) in +s
+        assert (p2 + 1.1*r*n) in +s

From 0db41d22810367ede353f128ed8b118ca226518b Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 10 Jul 2019 14:50:42 -0500
Subject: [PATCH 2/2] Update cylinder_from_points test based on @amandalund
 feedback

---
 tests/unit_tests/test_surface.py | 31 +++++++++++++++++++------------
 1 file changed, 19 insertions(+), 12 deletions(-)

diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py
index 09bfa738e3..1e3ae6f168 100644
--- a/tests/unit_tests/test_surface.py
+++ b/tests/unit_tests/test_surface.py
@@ -1,5 +1,7 @@
+from functools import partial
+from random import uniform, seed
+
 import numpy as np
-from random import random
 import openmc
 import pytest
 
@@ -330,11 +332,13 @@ def test_quadric():
 
 
 def test_cylinder_from_points():
-    for _ in range(10):
+    seed(1)  # Make random numbers reproducible
+    for _ in range(100):
         # Generate cylinder in random direction
-        p1 = np.array([random(), random(), random()])
-        p2 = np.array([random(), random(), random()])
-        r = random()
+        xi = partial(uniform, -10.0, 10.0)
+        p1 = np.array([xi(), xi(), xi()])
+        p2 = np.array([xi(), xi(), xi()])
+        r = uniform(1.0, 100.0)
         s = openmc.model.cylinder_from_points(p1, p2, r)
 
         # Points p1 and p2 need to be inside cylinder
@@ -342,16 +346,19 @@ def test_cylinder_from_points():
         assert p2 in -s
 
         # Points further along the line should be inside cylinder as well
-        t = 100*random() - 200
+        t = uniform(-100.0, 100.0)
         p = p1 + t*(p2 - p1)
         assert p in -s
 
-        # Check that a point outside cylinder is in positive half-space. We do
-        # this by constructing a plane that includes the cylinder's axis,
-        # finding the normal to the plane, and using it to find a point slightly
-        # more than one radius away from the axis.
+        # Check that points outside cylinder are in positive half-space and
+        # inside are in negative half-space. We do this by constructing a plane
+        # that includes the cylinder's axis, finding the normal to the plane,
+        # and using it to find a point slightly more/less than one radius away
+        # from the axis.
         plane = openmc.Plane.from_points(p1, p2, (0., 0., 0.))
         n = np.array([plane.a, plane.b, plane.c])
         n /= np.linalg.norm(n)
-        assert (p1 + 1.1*r*n) in +s
-        assert (p2 + 1.1*r*n) in +s
+        assert p1 + 1.1*r*n in +s
+        assert p2 + 1.1*r*n in +s
+        assert p1 + 0.9*r*n in -s
+        assert p2 + 0.9*r*n in -s