diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst
index 66440a3c37..ec10fea4fa 100644
--- a/docs/source/usersguide/install.rst
+++ b/docs/source/usersguide/install.rst
@@ -404,7 +404,7 @@ to install the Python package in :ref:`"editable" mode <devguide_editable>`.
 Prerequisites
 -------------
 
-The Python API works with Python 3.4+. In addition to Python itself, the API
+The Python API works with Python 3.5+. In addition to Python itself, the API
 relies on a number of third-party packages. All prerequisites can be installed
 using Conda_ (recommended), pip_, or through the package manager in most Linux
 distributions. To run simulations in parallel using MPI, it is recommended to
diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py
index 3ee5052402..9999c383ff 100644
--- a/openmc/model/funcs.py
+++ b/openmc/model/funcs.py
@@ -373,52 +373,7 @@ def get_hexagonal_prism(*args, **kwargs):
     return hexagonal_prism(*args, **kwargs)
 
 
-def cylinder_from_points(p1, p2, r, **kwargs):
-    """Return cylinder defined by two points passing through its center.
-
-    Parameters
-    ----------
-    p1, p2 : 3-tuples
-        Coordinates of two points that pass through the center of the cylinder
-    r : float
-        Radius of the cylinder
-    kwargs : dict
-        Keyword arguments passed to the :class:`openmc.Quadric` constructor
-
-    Returns
-    -------
-    openmc.Quadric
-        Quadric surface representing the cylinder.
-
-    """
-    # Get x, y, z coordinates of two points
-    x1, y1, z1 = p1
-    x2, y2, z2 = p2
-
-    # Define intermediate terms
-    dx = x2 - x1
-    dy = y2 - y1
-    dz = z2 - z1
-    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|.
-    # Expanding out all terms and grouping according to what Quadric expects
-    # gives the following coefficients.
-    kwargs['a'] = dy*dy + dz*dz
-    kwargs['b'] = dx*dx + dz*dz
-    kwargs['c'] = dx*dx + dy*dy
-    kwargs['d'] = -2*dx*dy
-    kwargs['e'] = -2*dy*dz
-    kwargs['f'] = -2*dx*dz
-    kwargs['g'] = 2*(cy*dz - cz*dy)
-    kwargs['h'] = 2*(cz*dx - cx*dz)
-    kwargs['j'] = 2*(cx*dy - cy*dx)
-    kwargs['k'] = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r
-
-    return Quadric(**kwargs)
+cylinder_from_points = Cylinder.from_points
 
 
 def subdivide(surfaces):
diff --git a/openmc/surface.py b/openmc/surface.py
index 92c828d65f..567968b2c8 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -3,13 +3,14 @@ from collections import OrderedDict
 from copy import deepcopy
 from numbers import Real
 from xml.etree import ElementTree as ET
-from warnings import warn
+from warnings import warn, catch_warnings, simplefilter
+import math
 
 import numpy as np
 
 from openmc.checkvalue import check_type, check_value
 from openmc.region import Region, Intersection, Union
-from openmc.mixin import IDManagerMixin
+from openmc.mixin import IDManagerMixin, IDWarning
 
 
 _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic', 'white']
@@ -361,6 +362,59 @@ class PlaneMixin(metaclass=ABCMeta):
         self._periodic_surface = periodic_surface
         periodic_surface._periodic_surface = self
 
+    def _get_base_coeffs(self):
+        return (self.a, self.b, self.c, self.d)
+
+    def _get_normal(self):
+        a, b, c = self._get_base_coeffs()[:3]
+        return np.array((a, b, c)) / math.sqrt(a*a + b*b + c*c)
+
+    def bounding_box(self, side):
+        """Determine an axis-aligned bounding box.
+
+        An axis-aligned bounding box for Plane half-spaces is represented by
+        its lower-left and upper-right coordinates. If the half-space is
+        unbounded in a particular direction, numpy.inf is used to represent
+        infinity.
+
+        Parameters
+        ----------
+        side : {'+', '-'}
+            Indicates the negative or positive half-space
+
+        Returns
+        -------
+        numpy.ndarray
+            Lower-left coordinates of the axis-aligned bounding box for the
+            desired half-space
+        numpy.ndarray
+            Upper-right coordinates of the axis-aligned bounding box for the
+            desired half-space
+
+        """
+        # Compute the bounding box based on the normal vector to the plane
+        nhat = self._get_normal()
+        ll = np.array([-np.inf, -np.inf, -np.inf])
+        ur = np.array([np.inf, np.inf, np.inf])
+        # If the plane is axis aligned, find the proper bounding box
+        if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)):
+            sign = nhat.sum()
+            a, b, c, d = self._get_base_coeffs()
+            vals = [d/val if not np.isclose(val, 0., rtol=0., atol=self._atol) 
+                    else np.nan for val in (a, b, c)]
+            if side == '-':
+                if sign > 0:
+                    ur = np.array([v if not np.isnan(v) else np.inf for v in vals])
+                else:
+                    ll = np.array([v if not np.isnan(v) else -np.inf for v in vals])
+            elif side == '+':
+                if sign > 0:
+                    ll = np.array([v if not np.isnan(v) else -np.inf for v in vals])
+                else:
+                    ur = np.array([v if not np.isnan(v) else np.inf for v in vals])
+
+        return (ll, ur)
+
     def evaluate(self, point):
         """Evaluate the surface equation at a given point.
 
@@ -509,6 +563,12 @@ class Plane(PlaneMixin, Surface):
                      FutureWarning)
                 setattr(self, k.lower(), val)
 
+    @classmethod
+    def __subclasshook__(cls, c):
+        if cls is Plane and c in (XPlane, YPlane, ZPlane):
+            return True
+        return NotImplemented
+
     @property
     def a(self):
         return self.coefficients['a']
@@ -545,9 +605,6 @@ class Plane(PlaneMixin, Surface):
         check_type('D coefficient', d, Real)
         self._coefficients['d'] = d
 
-    def _get_base_coeffs(self):
-        return (self.a, self.b, self.c, self.d)
-
     @classmethod
     def from_points(cls, p1, p2, p3, **kwargs):
         """Return a plane given three points that pass through it.
@@ -638,29 +695,31 @@ class XPlane(PlaneMixin, Surface):
     def x0(self):
         return self.coefficients['x0']
 
+    @property
+    def a(self):
+        return 1.
+
+    @property
+    def b(self):
+        return 0.
+
+    @property
+    def c(self):
+        return 0.
+
+    @property
+    def d(self):
+        return self.x0
+
     @x0.setter
     def x0(self, x0):
         check_type('x0 coefficient', x0, Real)
         self._coefficients['x0'] = x0
 
-    def _get_base_coeffs(self):
-        return (1., 0., 0., self.x0)
-
-    def bounding_box(self, side):
-        if side == '-':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([self.x0, np.inf, np.inf]))
-        elif side == '+':
-            return (np.array([self.x0, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
-
     def evaluate(self, point):
         return point[0] - self.x0
 
 
-Plane.register(XPlane)
-
-
 class YPlane(PlaneMixin, Surface):
     """A plane perpendicular to the y axis of the form :math:`y - y_0 = 0`
 
@@ -718,29 +777,31 @@ class YPlane(PlaneMixin, Surface):
     def y0(self):
         return self.coefficients['y0']
 
+    @property
+    def a(self):
+        return 0.
+
+    @property
+    def b(self):
+        return 1.
+
+    @property
+    def c(self):
+        return 0.
+
+    @property
+    def d(self):
+        return self.y0
+
     @y0.setter
     def y0(self, y0):
         check_type('y0 coefficient', y0, Real)
         self._coefficients['y0'] = y0
 
-    def _get_base_coeffs(self):
-        return (0., 1., 0., self.y0)
-
-    def bounding_box(self, side):
-        if side == '-':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, self.y0, np.inf]))
-        elif side == '+':
-            return (np.array([-np.inf, self.y0, -np.inf]),
-                    np.array([np.inf, np.inf, np.inf]))
-
     def evaluate(self, point):
         return point[1] - self.y0
 
 
-Plane.register(YPlane)
-
-
 class ZPlane(PlaneMixin, Surface):
     """A plane perpendicular to the z axis of the form :math:`z - z_0 = 0`
 
@@ -798,32 +859,43 @@ class ZPlane(PlaneMixin, Surface):
     def z0(self):
         return self.coefficients['z0']
 
+    @property
+    def a(self):
+        return 0.
+
+    @property
+    def b(self):
+        return 0.
+
+    @property
+    def c(self):
+        return 1.
+
+    @property
+    def d(self):
+        return self.z0
+
     @z0.setter
     def z0(self, z0):
         check_type('z0 coefficient', z0, Real)
         self._coefficients['z0'] = z0
 
-    def _get_base_coeffs(self):
-        return (0., 0., 1., self.z0)
-
-    def bounding_box(self, side):
-        if side == '-':
-            return (np.array([-np.inf, -np.inf, -np.inf]),
-                    np.array([np.inf, np.inf, self.z0]))
-        elif side == '+':
-            return (np.array([-np.inf, -np.inf, self.z0]),
-                    np.array([np.inf, np.inf, np.inf]))
-
     def evaluate(self, point):
         return point[2] - self.z0
 
 
-Plane.register(ZPlane)
-
-
 class QuadricMixin(metaclass=ABCMeta):
     """A Mixin class implementing common functionality for quadric surfaces"""
 
+    @property
+    def _origin(self):
+        return np.array((self.x0, self.y0, self.z0))
+
+    @property
+    def _axis(self):
+        axis = np.array((self.dx, self.dy, self.dz))
+        return axis / np.linalg.norm(axis)
+
     def get_Abc(self, coeffs=None):
         """Compute matrix, vector, and scalar coefficients for this surface or
         for a specified set of coefficients.
@@ -884,7 +956,7 @@ class QuadricMixin(metaclass=ABCMeta):
         """
         x = np.asarray(point)
         A, b, c = self.get_Abc()
-        return np.matmul(x.T, np.matmul(A, x)) + np.matmul(b.T, x) + c
+        return x.T @ A @ x + b.T @ x + c
 
     def translate(self, vector, inplace=False):
         """Translate surface in given direction
@@ -907,17 +979,20 @@ class QuadricMixin(metaclass=ABCMeta):
 
         surf = self if inplace else self.clone()
 
-        if set(('x0', 'y0', 'z0')).intersection(set(surf._coeff_keys)):
+        if hasattr(self, 'x0'):
             for vi, xi in zip(vector, ('x0', 'y0', 'z0')):
-                val = getattr(surf, xi, None)
-                if val is not None:
+                val = getattr(surf, xi)
+                try:
                     setattr(surf, xi, val + vi)
+                except AttributeError:
+                    # That attribute is read only i.e x0 for XCylinder
+                    pass
+
         else:
             A, bvec, cnst = self.get_Abc()
 
-            g, h, j = bvec - 2*np.matmul(vector.T, A)
-            k = cnst + np.matmul(vector.T, np.matmul(A, vector)) \
-                - np.matmul(bvec.T, vector)
+            g, h, j = bvec - 2*vector.T @ A
+            k = cnst + vector.T @ A @ vector - bvec.T @ vector
 
             for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)):
                 setattr(surf, key, val)
@@ -993,16 +1068,17 @@ class Cylinder(QuadricMixin, Surface):
     _coeff_keys = ('x0', 'y0', 'z0', 'r', 'dx', 'dy', 'dz')
 
     def __init__(self, x0=0., y0=0., z0=0., r=1., dx=0., dy=0., dz=1., **kwargs):
-        raise NotImplementedError('There is no C++ implementation for general '
-                                  'Cylinders yet, please use '
-                                  'openmc.model.funcs.cylinder_from_points to '
-                                  'return a Quadric instance instead for now')
-
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)):
             setattr(self, key, val)
 
+    @classmethod
+    def __subclasshook__(cls, c):
+        if cls is Cylinder and c in (XCylinder, YCylinder, ZCylinder):
+            return True
+        return NotImplemented
+
     @property
     def x0(self):
         return self.coefficients['x0']
@@ -1068,8 +1144,8 @@ class Cylinder(QuadricMixin, Surface):
 
     def _get_base_coeffs(self):
         # Get x, y, z coordinates of two points
-        x1, y1, z1 = self.x0, self.y0, self.z0
-        x2, y2, z2 = x1 + self.dx, y1 + self.dy, z1 + self.dz
+        x1, y1, z1 = self._origin
+        x2, y2, z2 = self._origin + self._axis
         r = self.r 
 
         # Define intermediate terms
@@ -1118,10 +1194,6 @@ class Cylinder(QuadricMixin, Surface):
             radius r.
 
         """
-        raise NotImplementedError('There is no C++ implementation for general '
-                                  'Cylinders yet, please use '
-                                  'openmc.model.funcs.cylinder_from_points to '
-                                  'return a Quadric instance instead for now')
         # Convert to numpy arrays
         p1 = np.asarray(p1)
         p2 = np.asarray(p2)
@@ -1130,6 +1202,24 @@ class Cylinder(QuadricMixin, Surface):
 
         return cls(x0=x0, y0=y0, z0=z0, r=r, dx=dx, dy=dy, dz=dz, **kwargs)
 
+    def to_xml_element(self):
+        """Return XML representation of the surface
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing source data
+
+        """
+        # This method overrides Surface.to_xml_element to generate a Quadric
+        # since the C++ layer doesn't support Cylinders right now
+        with catch_warnings():
+            simplefilter('ignore', IDWarning)
+            kwargs = {'boundary_type': self.boundary_type, 'name': self.name,
+                      'surface_id': self.id} 
+            quad_rep = Quadric(*self._get_base_coeffs(), **kwargs)
+        return quad_rep.to_xml_element()
+
 
 class XCylinder(QuadricMixin, Surface):
     """An infinite cylinder whose length is parallel to the x-axis of the form
@@ -1202,6 +1292,22 @@ class XCylinder(QuadricMixin, Surface):
     def r(self):
         return self.coefficients['r']
 
+    @property
+    def x0(self):
+        return 0.
+
+    @property
+    def dx(self):
+        return 1.
+
+    @property
+    def dy(self):
+        return 0.
+
+    @property
+    def dz(self):
+        return 0.
+
     @y0.setter
     def y0(self, y0):
         check_type('y0 coefficient', y0, Real)
@@ -1240,9 +1346,6 @@ class XCylinder(QuadricMixin, Surface):
         return y*y + z*z - self.r**2
 
 
-Cylinder.register(XCylinder)
-
-
 class YCylinder(QuadricMixin, Surface):
     """An infinite cylinder whose length is parallel to the y-axis of the form
     :math:`(x - x_0)^2 + (z - z_0)^2 = r^2`.
@@ -1314,6 +1417,22 @@ class YCylinder(QuadricMixin, Surface):
     def r(self):
         return self.coefficients['r']
 
+    @property
+    def y0(self):
+        return 0.
+
+    @property
+    def dx(self):
+        return 0.
+
+    @property
+    def dy(self):
+        return 1.
+
+    @property
+    def dz(self):
+        return 0.
+
     @x0.setter
     def x0(self, x0):
         check_type('x0 coefficient', x0, Real)
@@ -1352,9 +1471,6 @@ class YCylinder(QuadricMixin, Surface):
         return x*x + z*z - self.r**2
 
 
-Cylinder.register(YCylinder)
-
-
 class ZCylinder(QuadricMixin, Surface):
     """An infinite cylinder whose length is parallel to the z-axis of the form
     :math:`(x - x_0)^2 + (y - y_0)^2 = r^2`.
@@ -1426,6 +1542,22 @@ class ZCylinder(QuadricMixin, Surface):
     def r(self):
         return self.coefficients['r']
 
+    @property
+    def z0(self):
+        return 0.
+
+    @property
+    def dx(self):
+        return 0.
+
+    @property
+    def dy(self):
+        return 0.
+
+    @property
+    def dz(self):
+        return 1.
+
     @x0.setter
     def x0(self, x0):
         check_type('x0 coefficient', x0, Real)
@@ -1464,9 +1596,6 @@ class ZCylinder(QuadricMixin, Surface):
         return x*x + y*y - self.r**2
 
 
-Cylinder.register(ZCylinder)
-
-
 class Sphere(QuadricMixin, Surface):
     """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2`.
 
@@ -1656,9 +1785,6 @@ class Cone(QuadricMixin, Surface):
     _coeff_keys = ('x0', 'y0', 'z0', 'r2', 'dx', 'dy', 'dz')
 
     def __init__(self, x0=0., y0=0., z0=0., r2=1., dx=0., dy=0., dz=1., **kwargs):
-        raise NotImplementedError('There is no C++ implementation for general '
-                                  'Cones yet, this functionality should be '
-                                  'added soon.')
         R2 = kwargs.pop('R2', None)
         if R2 is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'),
@@ -1669,6 +1795,12 @@ class Cone(QuadricMixin, Surface):
         for key, val in zip(self._coeff_keys, (x0, y0, z0, r2, dx, dy, dz)):
             setattr(self, key, val)
 
+    @classmethod
+    def __subclasshook__(cls, c):
+        if cls is Cone and c in (XCone, YCone, ZCone):
+            return True
+        return NotImplemented
+
     @property
     def x0(self):
         return self.coefficients['x0']
@@ -1745,27 +1877,42 @@ class Cone(QuadricMixin, Surface):
         # The argument r2 for cones is actually tan^2(theta) so that
         # cos^2(theta) = 1 / (1 + r2)
 
-        x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2
-        dx, dy, dz = self.dx, self.dy, self.dz
-        dnorm = dx*dx + dy*dy + dz*dz
-        dx /= dnorm
-        dy /= dnorm
-        dz /= dnorm
-        cos2 = 1 / (1 + r2)
+        x0, y0, z0 = self._origin
+        dx, dy, dz = self._axis
+        cos2 = 1 / (1 + self.r2)
 
-        a = dx*dx - cos2
-        b = dy*dy - cos2
-        c = dz*dz - cos2
-        d = 2*dx*dy
-        e = 2*dy*dz
-        f = 2*dx*dz
-        g = -2*(dx*dx*x0 + dx*dy*y0 + dx*dz*z0 - cos2)
-        h = -2*(dy*dy*y0 + dx*dy*x0 + dy*dz*z0 - cos2)
-        j = -2*(dz*dz*y0 + dx*dz*x0 + dy*dz*y0 - cos2)
-        k = (dx*x0 + dy*y0 + dz*z0)**2 - cos2*(x0*x0 + y0*y0 + z0*z0)
+        a = cos2 - dx*dx
+        b = cos2 - dy*dy
+        c = cos2 - dz*dz
+        d = -2*dx*dy
+        e = -2*dy*dz
+        f = -2*dx*dz
+        g = 2*(dx*(dy*y0 + dz*z0) - a*x0)
+        h = 2*(dy*(dx*x0 + dz*z0) - b*y0)
+        j = 2*(dz*(dx*x0 + dy*y0) - c*z0)
+        k = a*x0*x0 + b*y0*y0 + c*z0*z0 - 2*(dx*dy*x0*y0 + dy*dz*y0*z0 +
+                                             dx*dz*x0*z0)
 
         return (a, b, c, d, e, f, g, h, j, k)
 
+    def to_xml_element(self):
+        """Return XML representation of the surface
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing source data
+
+        """
+        # This method overrides Surface.to_xml_element to generate a Quadric
+        # since the C++ layer doesn't support Cones right now
+        with catch_warnings():
+            simplefilter('ignore', IDWarning)
+            kwargs = {'boundary_type': self.boundary_type, 'name': self.name,
+                      'surface_id': self.id} 
+            quad_rep = Quadric(*self._get_base_coeffs(), **kwargs)
+        return quad_rep.to_xml_element()
+
 
 class XCone(QuadricMixin, Surface):
     """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 =
@@ -1845,6 +1992,18 @@ class XCone(QuadricMixin, Surface):
     def r2(self):
         return self.coefficients['r2']
 
+    @property
+    def dx(self):
+        return 1.
+
+    @property
+    def dy(self):
+        return 0.
+
+    @property
+    def dz(self):
+        return 0.
+
     @x0.setter
     def x0(self, x0):
         check_type('x0 coefficient', x0, Real)
@@ -1883,9 +2042,6 @@ class XCone(QuadricMixin, Surface):
         return y*y + z*z - self.r2*x*x
 
 
-Cone.register(XCone)
-
-
 class YCone(QuadricMixin, Surface):
     """A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 =
     r^2 (y - y_0)^2`.
@@ -1964,6 +2120,18 @@ class YCone(QuadricMixin, Surface):
     def r2(self):
         return self.coefficients['r2']
 
+    @property
+    def dx(self):
+        return 0.
+
+    @property
+    def dy(self):
+        return 1.
+
+    @property
+    def dz(self):
+        return 0.
+
     @x0.setter
     def x0(self, x0):
         check_type('x0 coefficient', x0, Real)
@@ -2002,9 +2170,6 @@ class YCone(QuadricMixin, Surface):
         return x*x + z*z - self.r2*y*y
 
 
-Cone.register(YCone)
-
-
 class ZCone(QuadricMixin, Surface):
     """A cone parallel to the x-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 =
     r^2 (z - z_0)^2`.
@@ -2083,6 +2248,18 @@ class ZCone(QuadricMixin, Surface):
     def r2(self):
         return self.coefficients['r2']
 
+    @property
+    def dx(self):
+        return 0.
+
+    @property
+    def dy(self):
+        return 0.
+
+    @property
+    def dz(self):
+        return 1.
+
     @x0.setter
     def x0(self, x0):
         check_type('x0 coefficient', x0, Real)
@@ -2121,9 +2298,6 @@ class ZCone(QuadricMixin, Surface):
         return x*x + y*y - self.r2*z*z
 
 
-Cone.register(ZCone)
-
-
 class Quadric(QuadricMixin, Surface):
     """A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy +
     Jz + K = 0`.
@@ -2452,5 +2626,4 @@ class Halfspace(Region):
         # Return translated surface
         return type(self)(memo[key], self.side)
 
-
 _SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()}
diff --git a/setup.py b/setup.py
index 99db831484..c193645c35 100755
--- a/setup.py
+++ b/setup.py
@@ -56,14 +56,13 @@ kwargs = {
         'Topic :: Scientific/Engineering'
         'Programming Language :: C++',
         'Programming Language :: Python :: 3',
-        'Programming Language :: Python :: 3.4',
         'Programming Language :: Python :: 3.5',
         'Programming Language :: Python :: 3.6',
         'Programming Language :: Python :: 3.7',
     ],
 
     # Dependencies
-    'python_requires': '>=3.4',
+    'python_requires': '>=3.5',
     'install_requires': [
         'numpy>=1.9', 'h5py', 'scipy', 'ipython', 'matplotlib',
         'pandas', 'lxml', 'uncertainties'
diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py
index 703010a059..749ae96212 100644
--- a/tests/unit_tests/test_surface.py
+++ b/tests/unit_tests/test_surface.py
@@ -2,6 +2,7 @@ from functools import partial
 from random import uniform, seed
 
 import numpy as np
+import math
 import openmc
 import pytest
 
@@ -139,6 +140,51 @@ def test_zplane():
     repr(s)
 
 
+def test_cylinder():
+    x0, y0, z0, r = 2, 3, 4, 2
+    dx, dy, dz = 1, -1, 1
+    s = openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=r)
+    assert s.x0 == 2
+    assert s.y0 == 3
+    assert s.z0 == 4
+    assert s.dx == 1
+    assert s.dy == -1
+    assert s.dz == 1
+    assert s.r == 2
+
+    # Check bounding box
+    assert_infinite_bb(s)
+
+    # evaluate method
+    # |(p - p1) ⨯ (p - p2)|^2 / |p2 - p1|^2 - r^2
+    p1 = s._origin
+    p2 = p1 + s._axis
+    perp = np.array((1, -2, 1))*(1 / s._axis)
+    divisor = np.linalg.norm(p2 - p1)
+    pin = p1 + 5*s._axis # point inside cylinder
+    pout = np.array((4., 0., 2.5)) # point outside the cylinder
+    pon = p1 + s.r*perp / np.linalg.norm(perp) # point on cylinder
+    for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)):
+        c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / divisor
+        val = c1*c1 - s.r*s.r
+        p_eval = s.evaluate(p)
+        assert fn(p_eval, 0.)
+        assert p_eval == pytest.approx(val)
+
+    # translate method
+    st = s.translate((1.0, 1.0, 1.0))
+    assert st.x0 == s.x0 + 1
+    assert st.y0 == s.y0 + 1
+    assert st.z0 == s.z0 + 1
+    assert st.dx == s.dx
+    assert st.dy == s.dy
+    assert st.dz == s.dz
+    assert st.r == s.r
+
+    # Make sure repr works
+    repr(s)
+
+
 def test_xcylinder():
     y, z, r = 3, 5, 2
     s = openmc.XCylinder(y0=y, z0=z, r=r)
@@ -284,6 +330,58 @@ def cone_common(apex, r2, cls):
     repr(s)
 
 
+def test_cone():
+    x0, y0, z0, r2 = 2, 3, 4, 4
+    dx, dy, dz = 1, -1, 1
+    s = openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=r2)
+    assert s.x0 == 2
+    assert s.y0 == 3
+    assert s.z0 == 4
+    assert s.dx == 1
+    assert s.dy == -1
+    assert s.dz == 1
+    assert s.r2 == 4
+
+    # Check bounding box
+    assert_infinite_bb(s)
+
+    # evaluate method
+    # cos^2(theta) * ((p - p1))**2 - (d @ (p - p1))^2 
+    # The argument r2 for cones is actually tan^2(theta) so that
+    # cos^2(theta) = 1 / (1 + r2)
+    #
+    # This makes the evaluation equation shown below where p is the evaluation
+    # point (x, y, z) p1 is the apex (origin) of the cone and r2 is related to
+    # the aperature of the cone as described above
+    # (p - p1) @ (p - p1) / (1 + r2) - (d @ (p - p1))^2
+    # point inside
+    p1 = s._origin
+    d = s._axis
+    perp = np.array((1, -2, 1))*(1 / d)
+    perp /= np.linalg.norm(perp)
+    pin = p1 + 5*d # point inside cone
+    pout = p1 + 3.2*perp # point outside cone
+    pon = p1 + 3.2*d + 3.2*math.sqrt(s.r2)*perp # point on cone
+    for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)):
+        val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2)
+        p_eval = s.evaluate(p)
+        assert fn(p_eval, 0.)
+        assert p_eval == pytest.approx(val)
+
+    # translate method
+    st = s.translate((1.0, 1.0, 1.0))
+    assert st.x0 == s.x0 + 1
+    assert st.y0 == s.y0 + 1
+    assert st.z0 == s.z0 + 1
+    assert st.dx == s.dx
+    assert st.dy == s.dy
+    assert st.dz == s.dz
+    assert st.r2 == s.r2
+
+    # Make sure repr works
+    repr(s)
+
+
 def test_xcone():
     apex = (10, 0, 0)
     r2 = 4
@@ -339,7 +437,7 @@ def test_cylinder_from_points():
         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)
+        s = openmc.Cylinder.from_points(p1, p2, r)
 
         # Points p1 and p2 need to be inside cylinder
         assert p1 in -s
@@ -369,24 +467,27 @@ def test_cylinder_from_points_axis():
 
     # (x - 3)^2 + (y - 4)^2 = 2^2
     # x^2 + y^2 - 6x - 8y + 21 = 0
-    s = openmc.model.cylinder_from_points((3., 4., 0.), (3., 4., 1.), 2.)
-    assert (s.a, s.b, s.c) == pytest.approx((1., 1., 0.))
-    assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
-    assert (s.g, s.h, s.j) == pytest.approx((-6., -8., 0.))
-    assert s.k == pytest.approx(21.)
+    s = openmc.Cylinder.from_points((3., 4., 0.), (3., 4., 1.), 2.)
+    a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
+    assert (a, b, c) == pytest.approx((1., 1., 0.))
+    assert (d, e, f) == pytest.approx((0., 0., 0.))
+    assert (g, h, j) == pytest.approx((-6., -8., 0.))
+    assert k == pytest.approx(21.)
 
     # (y + 7)^2 + (z - 1)^2 = 3^2
     # y^2 + z^2 + 14y - 2z + 41 = 0
-    s = openmc.model.cylinder_from_points((0., -7, 1.), (1., -7., 1.), 3.)
-    assert (s.a, s.b, s.c) == pytest.approx((0., 1., 1.))
-    assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
-    assert (s.g, s.h, s.j) == pytest.approx((0., 14., -2.))
-    assert s.k == 41.
+    s = openmc.Cylinder.from_points((0., -7, 1.), (1., -7., 1.), 3.)
+    a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
+    assert (a, b, c) == pytest.approx((0., 1., 1.))
+    assert (d, e, f) == pytest.approx((0., 0., 0.))
+    assert (g, h, j) == pytest.approx((0., 14., -2.))
+    assert k == 41.
 
     # (x - 2)^2 + (z - 5)^2 = 4^2
     # x^2 + z^2 - 4x - 10z + 13 = 0
-    s = openmc.model.cylinder_from_points((2., 0., 5.), (2., 1., 5.), 4.)
-    assert (s.a, s.b, s.c) == pytest.approx((1., 0., 1.))
-    assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
-    assert (s.g, s.h, s.j) == pytest.approx((-4., 0., -10.))
-    assert s.k == pytest.approx(13.)
+    s = openmc.Cylinder.from_points((2., 0., 5.), (2., 1., 5.), 4.)
+    a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
+    assert (a, b, c) == pytest.approx((1., 0., 1.))
+    assert (d, e, f) == pytest.approx((0., 0., 0.))
+    assert (g, h, j) == pytest.approx((-4., 0., -10.))
+    assert k == pytest.approx(13.)