diff --git a/CMakeLists.txt b/CMakeLists.txt
index 29eed4fb04..b761929523 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -42,8 +42,9 @@ if(dagmc)
   find_package(DAGMC REQUIRED)
   add_library(dagmc-imported INTERFACE IMPORTED)
   link_directories(${DAGMC_LIBRARY_DIRS})
-  target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES})
-  target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS})
+  set_target_properties(dagmc-imported PROPERTIES
+    INTERFACE_LINK_LIBRARIES "${DAGMC_LIBRARIES}"
+    INTERFACE_INCLUDE_DIRECTORIES "${DAGMC_INCLUDE_DIRS}")
 else()
   set(DAGMC_FOUND false)
 endif()
@@ -170,6 +171,8 @@ target_link_libraries(xtensor INTERFACE xtl)
 # GSL header-only library
 #===============================================================================
 
+set(GSL_LITE_OPT_INSTALL_COMPAT_HEADER ON CACHE BOOL
+  "Install MS-GSL compatibility header <gsl/gsl>")
 add_subdirectory(vendor/gsl-lite)
 
 # Make sure contract violations throw exceptions
diff --git a/openmc/region.py b/openmc/region.py
index 2b2fce14d0..22718ac7b2 100644
--- a/openmc/region.py
+++ b/openmc/region.py
@@ -225,7 +225,6 @@ class Region(metaclass=ABCMeta):
         # at the end
         return output[0]
 
-    @abstractmethod
     def clone(self, memo=None):
         """Create a copy of this region - each of the surfaces in the
         region's nodes will be cloned and will have new unique IDs.
@@ -242,9 +241,14 @@ class Region(metaclass=ABCMeta):
             The clone of this region
 
         """
-        pass
 
-    @abstractmethod
+        if memo is None:
+            memo = {}
+
+        clone = deepcopy(self)
+        clone[:] = [n.clone(memo) for n in self]
+        return clone
+
     def translate(self, vector, memo=None):
         """Translate region in given direction
 
@@ -262,7 +266,52 @@ class Region(metaclass=ABCMeta):
             Translated region
 
         """
-        pass
+
+        if memo is None:
+            memo = {}
+        return type(self)(n.translate(vector, memo) for n in self)
+
+    def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False,
+               memo=None):
+        r"""Rotate surface by angles provided or by applying matrix directly.
+
+        Parameters
+        ----------
+        rotation : 3-tuple of float, or 3x3 iterable
+            A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where
+            the first element is the rotation about the x-axis in the fixed
+            laboratory frame, the second element is the rotation about the
+            y-axis in the fixed laboratory frame, and the third element is the
+            rotation about the z-axis in the fixed laboratory frame. The
+            rotations are active rotations. Additionally a 3x3 rotation matrix
+            can be specified directly either as a nested iterable or array.
+        pivot : iterable of float, optional
+            (x, y, z) coordinates for the point to rotate about. Defaults to
+            (0., 0., 0.)
+        order : str, optional
+            A string of 'x', 'y', and 'z' in some order specifying which
+            rotation to perform first, second, and third. Defaults to 'xyz'
+            which means, the rotation by angle :math:`\phi` about x will be
+            applied first, followed by :math:`\theta` about y and then
+            :math:`\psi` about z. This corresponds to an x-y-z extrinsic
+            rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan
+            angles :math:`(\phi, \theta, \psi)`.
+        inplace : boolean
+            Whether or not to return a new instance of Surface or to modify the
+            coefficients of this Surface in place. Defaults to False.
+        memo : dict or None
+            Dictionary used for memoization
+
+        Returns
+        -------
+        openmc.Region
+            Translated region
+
+        """
+        if memo is None:
+            memo = {}
+        return type(self)(n.rotate(rotation, pivot=pivot, order=order,
+                                   inplace=inplace, memo=memo) for n in self)
 
 
 class Intersection(Region, MutableSequence):
@@ -354,51 +403,6 @@ class Intersection(Region, MutableSequence):
             upper_right[:] = np.minimum(upper_right, upper_right_n)
         return lower_left, upper_right
 
-    def clone(self, memo=None):
-        """Create a copy of this region - each of the surfaces in the
-        intersection's nodes will be cloned and will have new unique IDs.
-
-        Parameters
-        ----------
-        memo : dict or None
-            A nested dictionary of previously cloned objects. This parameter
-            is used internally and should not be specified by the user.
-
-        Returns
-        -------
-        clone : openmc.Intersection
-            The clone of this intersection
-
-        """
-
-        if memo is None:
-            memo = {}
-
-        clone = deepcopy(self)
-        clone[:] = [n.clone(memo) for n in self]
-        return clone
-
-    def translate(self, vector, memo=None):
-        """Translate region in given direction
-
-        Parameters
-        ----------
-        vector : iterable of float
-            Direction in which region should be translated
-        memo : dict or None
-            Dictionary used for memoization. This parameter is used internally
-            and should not be specified by the user.
-
-        Returns
-        -------
-        openmc.Intersection
-            Translated region
-
-        """
-        if memo is None:
-            memo = {}
-        return type(self)(n.translate(vector, memo) for n in self)
-
 
 class Union(Region, MutableSequence):
     r"""Union of two or more regions.
@@ -487,51 +491,6 @@ class Union(Region, MutableSequence):
             upper_right[:] = np.maximum(upper_right, upper_right_n)
         return lower_left, upper_right
 
-    def clone(self, memo=None):
-        """Create a copy of this region - each of the surfaces in the
-        union's nodes will be cloned and will have new unique IDs.
-
-        Parameters
-        ----------
-        memo : dict or None
-            A nested dictionary of previously cloned objects. This parameter
-            is used internally and should not be specified by the user.
-
-        Returns
-        -------
-        clone : openmc.Union
-            The clone of this union
-
-        """
-
-        if memo is None:
-            memo = {}
-
-        clone = deepcopy(self)
-        clone[:] = [n.clone(memo) for n in self]
-        return clone
-
-    def translate(self, vector, memo=None):
-        """Translate region in given direction
-
-        Parameters
-        ----------
-        vector : iterable of float
-            Direction in which region should be translated
-        memo : dict or None
-            Dictionary used for memoization. This parameter is used internally
-            and should not be specified by the user.
-
-        Returns
-        -------
-        openmc.Union
-            Translated region
-
-        """
-        if memo is None:
-            memo = {}
-        return type(self)(n.translate(vector, memo) for n in self)
-
 
 class Complement(Region):
     """Complement of a region.
@@ -642,22 +601,6 @@ class Complement(Region):
             region.remove_redundant_surfaces(redundant_surfaces)
 
     def clone(self, memo=None):
-        """Create a copy of this region - each of the surfaces in the
-        complement's node will be cloned and will have new unique IDs.
-
-        Parameters
-        ----------
-        memo : dict or None
-            A nested dictionary of previously cloned objects. This parameter
-            is used internally and should not be specified by the user.
-
-        Returns
-        -------
-        clone : openmc.Complement
-            The clone of this complement
-
-        """
-
         if memo is None:
             memo = {}
 
@@ -666,22 +609,13 @@ class Complement(Region):
         return clone
 
     def translate(self, vector, memo=None):
-        """Translate region in given direction
-
-        Parameters
-        ----------
-        vector : iterable of float
-            Direction in which region should be translated
-        memo : dict or None
-            Dictionary used for memoization. This parameter is used internally
-            and should not be specified by the user.
-
-        Returns
-        -------
-        openmc.Complement
-            Translated region
-
-        """
         if memo is None:
             memo = {}
         return type(self)(self.node.translate(vector, memo))
+
+    def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False,
+               memo=None):
+        if memo is None:
+            memo = {}
+        return type(self)(self.node.rotate(rotation, pivot=pivot, order=order,
+                                           inplace=inplace, memo=memo))
diff --git a/openmc/surface.py b/openmc/surface.py
index 567968b2c8..021f134b20 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -1,5 +1,6 @@
 from abc import ABCMeta, abstractmethod
 from collections import OrderedDict
+from collections.abc import Iterable
 from copy import deepcopy
 from numbers import Real
 from xml.etree import ElementTree as ET
@@ -8,7 +9,7 @@ import math
 
 import numpy as np
 
-from openmc.checkvalue import check_type, check_value
+from openmc.checkvalue import check_type, check_value, check_length
 from openmc.region import Region, Intersection, Union
 from openmc.mixin import IDManagerMixin, IDWarning
 
@@ -26,6 +27,53 @@ will not accept positional parameters for superclass arguments.\
 """
 
 
+def _future_kwargs_warning_helper(cls, *args, **kwargs):
+    # Warn if Surface parameters are passed by position, not by keyword
+    argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args))
+    for k in argsdict:
+        warn(_WARNING_KWARGS.format(cls.__name__, k), FutureWarning)
+    kwargs.update(argsdict)
+    return kwargs
+
+
+def get_rotation_matrix(rotation, order='xyz'):
+    r"""Generate a 3x3 rotation matrix from input angles
+
+    Parameters
+    ----------
+    rotation : 3-tuple of float
+        A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where the
+        first element is the rotation about the x-axis in the fixed laboratory
+        frame, the second element is the rotation about the y-axis in the fixed
+        laboratory frame, and the third element is the rotation about the
+        z-axis in the fixed laboratory frame. The rotations are active
+        rotations.
+    order : str, optional
+        A string of 'x', 'y', and 'z' in some order specifying which rotation
+        to perform first, second, and third. Defaults to 'xyz' which means, the
+        rotation by angle :math:`\phi` about x will be applied first, followed
+        by :math:`\theta` about y and then :math:`\psi` about z. This
+        corresponds to an x-y-z extrinsic rotation as well as a z-y'-x''
+        intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`.
+
+    """
+    check_type('surface rotation', rotation, Iterable, Real)
+    check_length('surface rotation', rotation, 3)
+
+    phi, theta, psi = np.array(rotation)*(math.pi/180.)
+    cx, sx = math.cos(phi), math.sin(phi)
+    cy, sy = math.cos(theta), math.sin(theta)
+    cz, sz = math.cos(psi), math.sin(psi)
+    R = {
+        'x': np.array([[1., 0., 0.], [0., cx, -sx], [0., sx, cx]]),
+        'y': np.array([[cy, 0., sy], [0., 1., 0.], [-sy, 0., cy]]),
+        'z': np.array([[cz, -sz, 0.], [sz, cz, 0.], [0., 0., 1.]]),
+    }
+
+    R1, R2, R3 = (R[xi] for xi in order)
+    return R3 @ R2 @ R1
+
+
 class Surface(IDManagerMixin, metaclass=ABCMeta):
     """An implicit surface with an associated boundary condition.
 
@@ -266,6 +314,42 @@ class Surface(IDManagerMixin, metaclass=ABCMeta):
 
         """
 
+    @abstractmethod
+    def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False):
+        r"""Rotate surface by angles provided or by applying matrix directly.
+
+        Parameters
+        ----------
+        rotation : 3-tuple of float, or 3x3 iterable
+            A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where
+            the first element is the rotation about the x-axis in the fixed
+            laboratory frame, the second element is the rotation about the
+            y-axis in the fixed laboratory frame, and the third element is the
+            rotation about the z-axis in the fixed laboratory frame. The
+            rotations are active rotations. Additionally a 3x3 rotation matrix
+            can be specified directly either as a nested iterable or array.
+        pivot : iterable of float, optional
+            (x, y, z) coordinates for the point to rotate about. Defaults to
+            (0., 0., 0.)
+        order : str, optional
+            A string of 'x', 'y', and 'z' in some order specifying which
+            rotation to perform first, second, and third. Defaults to 'xyz'
+            which means, the rotation by angle :math:`\phi` about x will be
+            applied first, followed by :math:`\theta` about y and then
+            :math:`\psi` about z. This corresponds to an x-y-z extrinsic
+            rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan
+            angles :math:`(\phi, \theta, \psi)`.
+        inplace : boolean
+            Whether or not to return a new instance of Surface or to modify the
+            coefficients of this Surface in place. Defaults to False.
+
+        Returns
+        -------
+        openmc.Surface
+            Rotated surface
+
+        """
+
     def to_xml_element(self):
         """Return XML representation of the surface
 
@@ -452,9 +536,11 @@ class PlaneMixin(metaclass=ABCMeta):
             Translated surface
 
         """
-        vx, vy, vz = vector
+        if np.allclose(vector, 0., rtol=0., atol=self._atol):
+            return self
+
         a, b, c, d = self._get_base_coeffs()
-        d = d + a*vx + b*vy + c*vz
+        d = d + np.dot([a, b, c], vector)
 
         surf = self if inplace else self.clone()
 
@@ -462,6 +548,32 @@ class PlaneMixin(metaclass=ABCMeta):
 
         return surf
 
+    def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False):
+        pivot = np.asarray(pivot)
+        rotation = np.asarray(rotation, dtype=float)
+
+        # Allow rotation matrix to be passed in directly, otherwise build it
+        if rotation.ndim == 2:
+            check_length('surface rotation', rotation.ravel(), 9)
+            Rmat = rotation
+        else:
+            Rmat = get_rotation_matrix(rotation, order=order)
+
+        # Translate surface to pivot
+        surf = self.translate(-pivot, inplace=inplace)
+
+        a, b, c, d = surf._get_base_coeffs()
+        # Compute new rotated coefficients a, b, c
+        a, b, c = Rmat @ [a, b, c]
+
+        kwargs = {'boundary_type': surf.boundary_type, 'name': surf.name}
+        if inplace:
+            kwargs['surface_id'] = surf.id
+
+        surf = Plane(a=a, b=b, c=c, d=d, **kwargs)
+
+        return surf.translate(pivot, inplace=inplace)
+
     def to_xml_element(self):
         """Return XML representation of the surface
 
@@ -538,30 +650,23 @@ class Plane(PlaneMixin, Surface):
         # *args should ultimately be limited to a, b, c, d as specified in
         # __init__, but to preserve the API it is allowed to accept Surface
         # parameters for now, but will raise warnings if this is done.
-        argtup = ('a', 'b', 'c', 'd', 'boundary_type', 'name', 'surface_id')
-        kwargs.update(dict(zip(argtup, args)))
-
-        # Warn if Surface parameters are passed by position, not by keyword
-        superkwargs = {}
-        for k in ('boundary_type', 'name', 'surface_id'):
-            val = kwargs.get(k, None)
-            if val is not None:
-                superkwargs[k] = val
-                warn(_WARNING_KWARGS.format(type(self), k),
-                     FutureWarning)
-
-        super().__init__(**superkwargs)
-
-        for key, val in zip(self._coeff_keys, (a, b, c, d)):
-            setattr(self, key, val)
-
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         # Warn if capital letter arguments are passed
+        capdict = {}
         for k in 'ABCD':
             val = kwargs.pop(k, None)
             if val is not None:
                 warn(_WARNING_UPPER.format(type(self), k.lower(), k),
                      FutureWarning)
-                setattr(self, k.lower(), val)
+                capdict[k.lower()] = val
+
+        super().__init__(**kwargs)
+
+        for key, val in zip(self._coeff_keys, (a, b, c, d)):
+            setattr(self, key, val)
+
+        for key, val in capdict.items():
+            setattr(self, key, val)
 
     @classmethod
     def __subclasshook__(cls, c):
@@ -683,11 +788,7 @@ class XPlane(PlaneMixin, Surface):
     def __init__(self, x0=0., *args, **kwargs):
         # work around for accepting Surface kwargs as positional parameters
         # until they are deprecated
-        argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args))
-        for k in argsdict:
-            warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning)
-        kwargs.update(argsdict)
-
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
         self.x0 = x0
 
@@ -765,11 +866,7 @@ class YPlane(PlaneMixin, Surface):
     def __init__(self, y0=0., *args, **kwargs):
         # work around for accepting Surface kwargs as positional parameters
         # until they are deprecated
-        argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args))
-        for k in argsdict:
-            warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning)
-        kwargs.update(argsdict)
-
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
         self.y0 = y0
 
@@ -847,11 +944,7 @@ class ZPlane(PlaneMixin, Surface):
     def __init__(self, z0=0., *args, **kwargs):
         # work around for accepting Surface kwargs as positional parameters
         # until they are deprecated
-        argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args))
-        for k in argsdict:
-            warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning)
-        kwargs.update(argsdict)
-
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
         self.z0 = z0
 
@@ -976,6 +1069,8 @@ class QuadricMixin(metaclass=ABCMeta):
 
         """
         vector = np.asarray(vector)
+        if np.allclose(vector, 0., rtol=0., atol=self._atol):
+            return self
 
         surf = self if inplace else self.clone()
 
@@ -999,6 +1094,55 @@ class QuadricMixin(metaclass=ABCMeta):
 
         return surf
 
+    def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False):
+        # Get pivot and rotation matrix 
+        pivot = np.asarray(pivot)
+        rotation = np.asarray(rotation, dtype=float)
+
+        # Allow rotaiton matrix to be passed in directly, otherwise build it
+        if rotation.ndim == 2:
+            check_length('surface rotation', rotation.ravel(), 9)
+            Rmat = rotation
+        else:
+            Rmat = get_rotation_matrix(rotation, order=order)
+
+        # Translate surface to the pivot point
+        tsurf = self.translate(-pivot, inplace=inplace)
+
+        # If the surface is already generalized just clone it
+        if type(tsurf) is tsurf._virtual_base:
+            surf = tsurf if inplace else tsurf.clone()
+        else:
+            base_cls = type(tsurf)._virtual_base
+            # Copy necessary surface attributes to new kwargs dictionary
+            kwargs = {'boundary_type': tsurf.boundary_type, 'name': tsurf.name}
+            if inplace:
+                kwargs['surface_id'] = tsurf.id
+            kwargs.update({k: getattr(tsurf, k) for k in base_cls._coeff_keys})
+            # Create new instance of the virtual base class
+            surf = base_cls(**kwargs)
+
+        # Perform rotations on axis, origin, or quadric coefficients
+        if hasattr(surf, 'dx'):
+            for key, val in zip(('dx', 'dy', 'dz'), Rmat @ tsurf._axis):
+                setattr(surf, key, val)
+        if hasattr(surf, 'x0'):
+            for key, val in zip(('x0', 'y0', 'z0'), Rmat @ tsurf._origin):
+                setattr(surf, key, val)
+        else:
+            A, bvec, k = surf.get_Abc()
+            Arot = Rmat @ A @ Rmat.T
+
+            a, b, c = np.diagonal(Arot)
+            d, e, f = 2*Arot[0, 1], 2*Arot[1, 2], 2*Arot[0, 2]
+            g, h, j = Rmat @ bvec
+
+            for key, val in zip(surf._coeff_keys, (a, b, c, d, e, f, g, h, j, k)):
+                setattr(surf, key, val)
+
+        # translate back to the original frame and return the surface
+        return surf.translate(pivot, inplace=inplace)
+
 
 class Cylinder(QuadricMixin, Surface):
     """A cylinder with radius r, centered on the point (x0, y0, z0) with an
@@ -1067,7 +1211,9 @@ class Cylinder(QuadricMixin, Surface):
     _type = 'cylinder'
     _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):
+    def __init__(self, x0=0., y0=0., z0=0., r=1., dx=0., dy=0., dz=1., *args,
+                 **kwargs):
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)):
@@ -1142,6 +1288,19 @@ class Cylinder(QuadricMixin, Surface):
         check_type('dz coefficient', dz, Real)
         self._coefficients['dz'] = dz
 
+    def bounding_box(self, side):
+        if side == '-':
+            r = self.r
+            ll = [xi - r if np.isclose(dxi, 0., rtol=0., atol=self._atol)
+                  else -np.inf for xi, dxi in zip(self._origin, self._axis)]
+            ur = [xi + r if np.isclose(dxi, 0., rtol=0., atol=self._atol)
+                  else np.inf for xi, dxi in zip(self._origin, self._axis)]
+            return (np.array(ll), np.array(ur))
+
+        elif side == '+':
+            return (np.array([-np.inf, -np.inf, -np.inf]),
+                    np.array([np.inf, np.inf, np.inf]))
+
     def _get_base_coeffs(self):
         # Get x, y, z coordinates of two points
         x1, y1, z1 = self._origin
@@ -1269,12 +1428,13 @@ class XCylinder(QuadricMixin, Surface):
     _type = 'x-cylinder'
     _coeff_keys = ('y0', 'z0', 'r')
 
-    def __init__(self, y0=0., z0=0., r=1., **kwargs):
+    def __init__(self, y0=0., z0=0., r=1., *args, **kwargs):
         R = kwargs.pop('R', None)
         if R is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'),
                  FutureWarning)
             r = R
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (y0, z0, r)):
@@ -1394,12 +1554,13 @@ class YCylinder(QuadricMixin, Surface):
     _type = 'y-cylinder'
     _coeff_keys = ('x0', 'z0', 'r')
 
-    def __init__(self, x0=0., z0=0., r=1., **kwargs):
+    def __init__(self, x0=0., z0=0., r=1., *args, **kwargs):
         R = kwargs.pop('R', None)
         if R is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'),
                  FutureWarning)
             r = R
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, z0, r)):
@@ -1519,12 +1680,13 @@ class ZCylinder(QuadricMixin, Surface):
     _type = 'z-cylinder'
     _coeff_keys = ('x0', 'y0', 'r')
 
-    def __init__(self, x0=0., y0=0., r=1., **kwargs):
+    def __init__(self, x0=0., y0=0., r=1., *args, **kwargs):
         R = kwargs.pop('R', None)
         if R is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'),
                  FutureWarning)
             r = R
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, y0, r)):
@@ -1646,12 +1808,13 @@ class Sphere(QuadricMixin, Surface):
     _type = 'sphere'
     _coeff_keys = ('x0', 'y0', 'z0', 'r')
 
-    def __init__(self, x0=0., y0=0., z0=0., r=1., **kwargs):
+    def __init__(self, x0=0., y0=0., z0=0., r=1., *args, **kwargs):
         R = kwargs.pop('R', None)
         if R is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'),
                  FutureWarning)
             r = R
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, y0, z0, r)):
@@ -1784,12 +1947,14 @@ class Cone(QuadricMixin, Surface):
     _type = 'cone'
     _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):
+    def __init__(self, x0=0., y0=0., z0=0., r2=1., dx=0., dy=0., dz=1., *args,
+                 **kwargs):
         R2 = kwargs.pop('R2', None)
         if R2 is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'),
                  FutureWarning)
             r2 = R2
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, y0, z0, r2, dx, dy, dz)):
@@ -1965,12 +2130,13 @@ class XCone(QuadricMixin, Surface):
     _type = 'x-cone'
     _coeff_keys = ('x0', 'y0', 'z0', 'r2')
 
-    def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs):
+    def __init__(self, x0=0., y0=0., z0=0., r2=1., *args, **kwargs):
         R2 = kwargs.pop('R2', None)
         if R2 is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'),
                  FutureWarning)
             r2 = R2
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)):
@@ -2093,12 +2259,13 @@ class YCone(QuadricMixin, Surface):
     _type = 'y-cone'
     _coeff_keys = ('x0', 'y0', 'z0', 'r2')
 
-    def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs):
+    def __init__(self, x0=0., y0=0., z0=0., r2=1., *args, **kwargs):
         R2 = kwargs.pop('R2', None)
         if R2 is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'),
                  FutureWarning)
             r2 = R2
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)):
@@ -2221,12 +2388,13 @@ class ZCone(QuadricMixin, Surface):
     _type = 'z-cone'
     _coeff_keys = ('x0', 'y0', 'z0', 'r2')
 
-    def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs):
+    def __init__(self, x0=0., y0=0., z0=0., r2=1., *args, **kwargs):
         R2 = kwargs.pop('R2', None)
         if R2 is not None:
             warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'),
                  FutureWarning)
             r2 = R2
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)):
@@ -2338,8 +2506,8 @@ class Quadric(QuadricMixin, Surface):
     _coeff_keys = ('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k')
 
     def __init__(self, a=0., b=0., c=0., d=0., e=0., f=0., g=0., h=0., j=0.,
-                 k=0., **kwargs):
-
+                 k=0., *args, **kwargs):
+        kwargs = _future_kwargs_warning_helper(type(self), *args, **kwargs)
         super().__init__(**kwargs)
 
         for key, val in zip(self._coeff_keys, (a, b, c, d, e, f, g, h, j, k)):
@@ -2623,7 +2791,74 @@ class Halfspace(Region):
         if key not in memo:
             memo[key] = self.surface.translate(vector)
 
-        # Return translated surface
+        # Return translated half-space
         return type(self)(memo[key], self.side)
 
+    def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False,
+               memo=None):
+        r"""Rotate surface by angles provided or by applying matrix directly.
+
+        Parameters
+        ----------
+        rotation : 3-tuple of float, or 3x3 iterable
+            A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where
+            the first element is the rotation about the x-axis in the fixed
+            laboratory frame, the second element is the rotation about the
+            y-axis in the fixed laboratory frame, and the third element is the
+            rotation about the z-axis in the fixed laboratory frame. The
+            rotations are active rotations. Additionally a 3x3 rotation matrix
+            can be specified directly either as a nested iterable or array.
+        pivot : iterable of float, optional
+            (x, y, z) coordinates for the point to rotate about. Defaults to
+            (0., 0., 0.)
+        order : str, optional
+            A string of 'x', 'y', and 'z' in some order specifying which
+            rotation to perform first, second, and third. Defaults to 'xyz'
+            which means, the rotation by angle :math:`\phi` about x will be
+            applied first, followed by :math:`\theta` about y and then
+            :math:`\psi` about z. This corresponds to an x-y-z extrinsic
+            rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan
+            angles :math:`(\phi, \theta, \psi)`.
+        inplace : boolean
+            Whether or not to return a new instance of Surface or to modify the
+            coefficients of this Surface in place. Defaults to False.
+        memo : dict or None
+            Dictionary used for memoization
+
+        Returns
+        -------
+        openmc.Halfspace
+            Translated half-space
+
+        """
+        if memo is None:
+            memo = {}
+
+        # If rotated surface not in memo, add it
+        key = (self.surface, tuple(rotation), tuple(pivot), order, inplace)
+        if key not in memo:
+            memo[key] = self.surface.rotate(rotation, pivot=pivot, order=order,
+                                            inplace=inplace)
+
+        # Return rotated half-space
+        return type(self)(memo[key], self.side)
+
+
 _SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()}
+
+
+# Set virtual base classes for "casting" up the heirarchy
+Plane._virtual_base = Plane
+XPlane._virtual_base = Plane
+YPlane._virtual_base = Plane
+ZPlane._virtual_base = Plane
+Cylinder._virtual_base = Cylinder
+XCylinder._virtual_base = Cylinder
+YCylinder._virtual_base = Cylinder
+ZCylinder._virtual_base = Cylinder
+Cone._virtual_base = Cone
+XCone._virtual_base = Cone
+YCone._virtual_base = Cone
+ZCone._virtual_base = Cone
+Sphere._virtual_base = Sphere
+Quadric._virtual_base = Quadric
diff --git a/tests/unit_tests/test_region.py b/tests/unit_tests/test_region.py
index f9656ebb4b..8b2e06ad1b 100644
--- a/tests/unit_tests/test_region.py
+++ b/tests/unit_tests/test_region.py
@@ -40,6 +40,13 @@ def test_union(reset):
     assert (6, 0, 0) not in regt
     assert (8, 0, 0) in regt
 
+    # rotate method
+    regr = region.rotate((0., 90., 0.))
+    assert (-4, 0, 0) not in regr
+    assert (0, 0, 6) in regr
+    assert (0, 0, -6) in regr
+    assert (0, 0, 3) not in regr
+
 
 def test_intersection(reset):
     s1 = openmc.XPlane(x0=5, surface_id=1)
@@ -73,6 +80,13 @@ def test_intersection(reset):
     assert (6, 0, 0) in regt
     assert (8, 0, 0) not in regt
 
+    # rotate method
+    regr = region.rotate((0., 90., 0.))
+    assert (-4, 0, 0) in regr
+    assert (0, 0, 6) not in regr
+    assert (0, 0, -6) not in regr
+    assert (0, 0, 3) in regr
+
 
 def test_complement(reset):
     zcyl = openmc.ZCylinder(r=1., surface_id=1)
@@ -106,6 +120,14 @@ def test_complement(reset):
     assert ll == pytest.approx((0., 0., -4.))
     assert ur == pytest.approx((2., 2., 6.))
 
+    # rotate method
+    inside_r = inside.rotate((90., 0., 0.))
+    ll, ur = inside_r.bounding_box
+    assert (.5, 2, 0) in inside_r
+    assert (0, 0, 6) not in inside_r
+    assert ll == pytest.approx((-1., -5., -1.))
+    assert ur == pytest.approx((1., 5., 1.))
+
 
 def test_get_surfaces():
     s1 = openmc.XPlane()
diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py
index 749ae96212..b37dc8f310 100644
--- a/tests/unit_tests/test_surface.py
+++ b/tests/unit_tests/test_surface.py
@@ -37,6 +37,17 @@ def test_plane():
     st = s.translate((1.0, 0.0, 0.0))
     assert (st.a, st.b, st.c, st.d) == (s.a, s.b, s.c, 4)
 
+    # rotate method
+    yp = openmc.YPlane(abs(s.d)/math.sqrt(s.a**2 + s.b**2 + s.c**2))
+    psi = math.degrees(math.atan2(1, 2))
+    phi = math.degrees(math.atan2(1, math.sqrt(5)))
+    sr = s.rotate((phi, 0., psi), order='zyx')
+    assert yp.normalize() == pytest.approx(sr.normalize())
+    # test rotation ordering
+    phi = math.degrees(math.atan2(1, math.sqrt(2)))
+    sr = s.rotate((0., -45., phi), order='xyz')
+    assert yp.normalize() == pytest.approx(sr.normalize())
+
     # Make sure repr works
     repr(s)
 
@@ -81,6 +92,14 @@ def test_xplane():
     st = s.translate((1.0, 0.0, 0.0))
     assert st.x0 == s.x0 + 1
 
+    # rotate method
+    # make sure rotating around x axis does nothing to coefficients
+    sr = s.rotate((37.4, 0., 0.))
+    assert s._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+    # rotating around z by 90 deg then x by -90 deg should give negative z-plane
+    sr = s.rotate((-90., 0., 90), order='zyx')
+    assert (0., 0., -1., 3.) == pytest.approx(sr._get_base_coeffs())
+
     # Make sure repr works
     repr(s)
 
@@ -110,6 +129,17 @@ def test_yplane():
     st = s.translate((0.0, 1.0, 0.0))
     assert st.y0 == s.y0 + 1
 
+    # rotate method
+    # make sure rotating around y axis does nothing to coefficients
+    sr = s.rotate((0., -12.4, 0.), order='yxz')
+    assert s._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+    # rotate around x by -90 deg and y by 90 deg should give negative x-plane
+    sr = s.rotate((-90, 90, 0.))
+    assert (-1, 0., 0., 3.) == pytest.approx(sr._get_base_coeffs())
+
+    # Make sure repr works
+    repr(s)
+
 
 def test_zplane():
     s = openmc.ZPlane(z0=3.)
@@ -136,6 +166,14 @@ def test_zplane():
     st = s.translate((0.0, 0.0, 1.0))
     assert st.z0 == s.z0 + 1
 
+    # rotate method
+    # make sure rotating around z axis does nothing to coefficients
+    sr = s.rotate((0., 0., 123), order='zxy')
+    assert s._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+    # rotate around x by -90 deg and y by 90 deg should give negative x-plane
+    sr = s.rotate((-90, 0., 90.))
+    assert (-1., 0., 0., 3.) == pytest.approx(sr._get_base_coeffs())
+
     # Make sure repr works
     repr(s)
 
@@ -181,6 +219,15 @@ def test_cylinder():
     assert st.dz == s.dz
     assert st.r == s.r
 
+    # rotate method
+    sr = s.rotate((90, 90, 90))
+    R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
+    assert sr._origin == pytest.approx(R @ s._origin)
+    assert sr._axis == pytest.approx(R @ s._axis)
+    # test passing in rotation matrix
+    sr2 = s.rotate(R)
+    assert sr2.is_equal(sr)
+
     # Make sure repr works
     repr(s)
 
@@ -209,6 +256,15 @@ def test_xcylinder():
     assert st.z0 == s.z0 + 1
     assert st.r == s.r
 
+    # rotate method
+    sr = s.rotate((90, 90, 90))
+    R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
+    assert sr._origin == pytest.approx(R @ s._origin)
+    assert sr._axis == pytest.approx(R @ s._axis)
+    # test passing in rotation matrix
+    sr2 = s.rotate(R)
+    assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+
     # Make sure repr works
     repr(s)
 
@@ -246,6 +302,18 @@ def test_ycylinder():
     assert st.z0 == s.z0 + 1
     assert st.r == s.r
 
+    # rotate method
+    sr = s.rotate((90, 90, 90))
+    R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
+    assert sr._origin == pytest.approx(R @ s._origin)
+    assert sr._axis == pytest.approx(R @ s._axis)
+    # test passing in rotation matrix
+    sr2 = s.rotate(R)
+    assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+
+    # Make sure repr works
+    repr(s)
+
 
 def test_zcylinder():
     x, y, r = 3, 5, 2
@@ -271,6 +339,15 @@ def test_zcylinder():
     assert st.y0 == s.y0 + 1
     assert st.r == s.r
 
+    # rotate method
+    sr = s.rotate((90, 90, 90))
+    R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
+    assert sr._origin == pytest.approx(R @ s._origin)
+    assert sr._axis == pytest.approx(R @ s._axis)
+    # test passing in rotation matrix
+    sr2 = s.rotate(R)
+    assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+
     # Make sure repr works
     repr(s)
 
@@ -301,6 +378,15 @@ def test_sphere():
     assert st.z0 == s.z0 + 1
     assert st.r == s.r
 
+    # rotate method
+    pivot = np.array([1, -2, 3])
+    sr = s.rotate((90, 90, 90), pivot=pivot)
+    R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
+    assert sr._origin == pytest.approx((R @ (s._origin - pivot)) + pivot)
+    # test passing in rotation matrix
+    sr2 = s.rotate(R, pivot=pivot)
+    assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+
     # Make sure repr works
     repr(s)
 
@@ -326,6 +412,15 @@ def cone_common(apex, r2, cls):
     assert st.z0 == s.z0 + 1
     assert st.r2 == s.r2
 
+    # rotate method
+    sr = s.rotate((90, 90, 90))
+    R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
+    assert sr._origin == pytest.approx(R @ s._origin)
+    assert sr._axis == pytest.approx(R @ s._axis)
+    # test passing in rotation matrix
+    sr2 = s.rotate(R)
+    assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+
     # Make sure repr works
     repr(s)
 
@@ -378,6 +473,15 @@ def test_cone():
     assert st.dz == s.dz
     assert st.r2 == s.r2
 
+    # rotate method
+    sr = s.rotate((90, 90, 90))
+    R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
+    assert sr._origin == pytest.approx(R @ s._origin)
+    assert sr._axis == pytest.approx(R @ s._axis)
+    # test passing in rotation matrix
+    sr2 = s.rotate(R)
+    assert sr2._get_base_coeffs() == pytest.approx(sr._get_base_coeffs())
+
     # Make sure repr works
     repr(s)
 
@@ -409,6 +513,7 @@ def test_quadric():
     assert s.b == coeffs['b']
     assert s.c == coeffs['c']
     assert s.k == coeffs['k']
+    assert openmc.Sphere(r=10).is_equal(s)
 
     # All other coeffs should be zero
     for coeff in ('d', 'e', 'f', 'g', 'h', 'j'):
@@ -428,6 +533,15 @@ def test_quadric():
     assert (st.g, st.h, st.j) == (-2, -2, -2)
     assert st.k == s.k + 3
 
+    # rotate method
+    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)
+    q = openmc.Quadric(*s._get_base_coeffs())
+    qr = q.rotate((45, 60, 30))
+    sr = s.rotate((45, 60, 30))
+    assert qr.is_equal(sr)
+
 
 def test_cylinder_from_points():
     seed(1)  # Make random numbers reproducible