Support Spheres

doc/source/api/modules.rst

@@ -76,6 +76,19 @@ Planes
    montepy.AxisPlane
 
 
+Spheres
+^^^^^^^
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myclass.rst
+
+   montepy.GeneralSphere
+   montepy.SphereAtOrigin
+   montepy.SphereOnAxis
+
+
 Data Inputs
 -----------
 

doc/source/changelog.rst

@@ -5,9 +5,17 @@ MontePy Changelog
 1.2 releases
 ============
 
-1.2.0
+#Next Version#
 --------------
 
+**Features Added**
+
+* Added support for ``Surface`` subtypes for spheres (:issue:`876`).
+
+
+1.2.0
+-----
+
 **Performance Improvement**
 
 * Optimized :math:`\mathcal{O}(N^2)` scaling in :func:`~montepy.numbered_object_collection.NumberedObjectCollection.request_number` by improving ``NumberedObjectCollection.check_number`` to :math:`\mathcal{O}(N)` (:issue:`786`). 

montepy/mcnp_object.py

@@ -225,13 +225,13 @@ def _update_values(self):
         """
         pass
 
-    def format_for_mcnp_input(self, mcnp_version: tuple[int]) -> list[str]:
+    def format_for_mcnp_input(self, mcnp_version: tuple[int, int, int]) -> list[str]:
         """Creates a list of strings representing this MCNP_Object that can be
         written to file.
 
         Parameters
         ----------
-        mcnp_version : tuple[int]
+        mcnp_version : tuple[int, int, int]
             The tuple for the MCNP version that must be exported to.
 
         Returns
@@ -242,7 +242,6 @@ def format_for_mcnp_input(self, mcnp_version: tuple[int]) -> list[str]:
         self.validate()
         self._update_values()
         self._tree.check_for_graveyard_comments()
-        message = None
         with warnings.catch_warnings(record=True) as ws:
             lines = self.wrap_string_for_mcnp(self._tree.format(), mcnp_version, True)
         self._flush_line_expansion_warning(lines, ws)

montepy/surfaces/__init__.py

@@ -1,4 +1,4 @@
-# Copyright 2024, Battelle Energy Alliance, LLC All Rights Reserved.
+# Copyright 2024 - 2026, Battelle Energy Alliance, LLC All Rights Reserved.
 from . import axis_plane
 from . import cylinder_par_axis
 from . import cylinder_on_axis
@@ -10,6 +10,9 @@
 from .axis_plane import AxisPlane
 from .cylinder_par_axis import CylinderParAxis
 from .cylinder_on_axis import CylinderOnAxis
+from .sphere_on_axis import SphereOnAxis
+from .sphere_at_origin import SphereAtOrigin
+from .general_sphere import GeneralSphere
 from .general_plane import GeneralPlane
 from .half_space import HalfSpace, UnitHalfSpace
 from .surface import Surface

montepy/surfaces/general_sphere.py

@@ -0,0 +1,95 @@
+# Copyright 2026, Battelle Energy Alliance, LLC All Rights Reserved.
+from .surface_type import SurfaceType
+from .surface import Surface, InitInput
+from montepy.exceptions import *
+from montepy.utilities import *
+
+from numbers import Real
+from typing import Union
+
+
+def _enforce_positive_radius(self, value):
+    if value < 0.0:
+        raise ValueError(f"Radius must be positive. {value} given")
+
+
+class GeneralSphere(Surface):
+    """Represents surface S: a general sphere
+
+    .. versionadded:: 1.3.0
+
+    Parameters
+    ----------
+    input : Union[Input, str]
+        The Input object representing the input
+    number : int
+        The number to set for this object.
+    """
+
+    def __init__(
+        self,
+        input: InitInput = None,
+        number: int = None,
+    ):
+        self._coordinates = [
+            self._generate_default_node(float, None),
+            self._generate_default_node(float, None),
+            self._generate_default_node(float, None),
+        ]
+        self._radius = self._generate_default_node(float, None)
+        super().__init__(input, number)
+        if input and self.surface_type != SurfaceType.S:
+            raise ValueError(f"A {self.__class__.__name__} must be a surface of type S")
+        if len(self.surface_constants) != 4:
+            raise ValueError(
+                f"A {self.__class__.__name__} must have exactly 4 surface_constants"
+            )
+        self._coordinates = self._surface_constants[:3]
+        self._radius = self._surface_constants[3]
+
+    @staticmethod
+    def _number_of_params():
+        return 4
+
+    @make_prop_val_node(
+        "_radius", (float, int), float, validator=_enforce_positive_radius
+    )
+    def radius(self):
+        """The radius of the sphere
+
+        Returns
+        -------
+        float
+        """
+        pass
+
+    @property
+    def coordinates(self):
+        """The three coordinates for the sphere center
+
+        :rytpe: tuple
+        """
+        return tuple(c.value for c in self._coordinates)
+
+    @coordinates.setter
+    def coordinates(self, coordinates):
+        if not isinstance(coordinates, (list, tuple)):
+            raise TypeError("coordinates must be a list")
+        if len(coordinates) != 3:
+            raise ValueError("coordinates must have exactly three elements")
+        for val in coordinates:
+            if not isinstance(val, Real):
+                raise TypeError(f"Coordinate must be a number. {val} given.")
+        for i, val in enumerate(coordinates):
+            self._coordinates[i].value = val
+
+    def validate(self):
+        super().validate()
+        if self.radius is None:  # pragma: no cover
+            raise IllegalState(f"Surface: {self.number} does not have a radius set.")
+        if any({c is None for c in self.coordinates}):  # pragma: no cover
+            raise IllegalState(f"Surface: {self.number} does not have coordinates set.")
+
+    def find_duplicate_surfaces(self, surfaces, tolerance):
+        """Duplicate sphere finding is not yet implemented"""
+        return []

montepy/surfaces/sphere_at_origin.py

@@ -0,0 +1,61 @@
+# Copyright 2026, Battelle Energy Alliance, LLC All Rights Reserved.
+from .surface_type import SurfaceType
+from .surface import Surface, InitInput
+from montepy.exceptions import *
+from montepy.utilities import *
+
+from typing import Union
+
+
+def _enforce_positive_radius(self, value):
+    if value < 0.0:
+        raise ValueError(f"Radius must be positive. {value} given")
+
+
+class SphereAtOrigin(Surface):
+    """Represents surface SO: a sphere at the origin
+
+    .. versionadded:: 1.3.0
+
+    Parameters
+    ----------
+    input : Union[Input, str]
+        The Input object representing the input
+    number : int
+        The number to set for this object.
+    """
+
+    def __init__(
+        self,
+        input: InitInput = None,
+        number: int = None,
+    ):
+        self._location = self._generate_default_node(float, None)
+        self._radius = self._generate_default_node(float, None)
+        super().__init__(input, number)
+        if len(self.surface_constants) != 1:
+            raise ValueError(
+                f"{self.__class__.__name__} must have exactly 1 surface_constant"
+            )
+        self._radius = self._surface_constants[0]
+
+    @make_prop_val_node(
+        "_radius", (float, int), float, validator=_enforce_positive_radius
+    )
+    def radius(self):
+        """The radius of the sphere
+
+        Returns
+        -------
+        float
+        """
+        pass
+
+    def validate(self):
+        super().validate()
+        if self.radius is None:  # pragma: no cover
+            raise IllegalState(f"Surface: {self.number} does not have a radius set.")
+
+    def find_duplicate_surfaces(self, surfaces, tolerance):
+        """Duplicate sphere finding is not yet implemented"""
+        return []

montepy/surfaces/sphere_on_axis.py

@@ -0,0 +1,86 @@
+# Copyright 2026, Battelle Energy Alliance, LLC All Rights Reserved.
+from .surface_type import SurfaceType
+from .surface import Surface, InitInput
+from montepy.exceptions import *
+from montepy.utilities import *
+
+from typing import Union
+
+
+def _enforce_positive_radius(self, value):
+    if value < 0.0:
+        raise ValueError(f"Radius must be positive. {value} given")
+
+
+class SphereOnAxis(Surface):
+    """Represents surfaces SX, SY, and SZ: spheres on axes
+
+    .. versionadded:: 1.3.0
+
+    Parameters
+    ----------
+    input : Union[Input, str]
+        The Input object representing the input
+    number : int
+        The number to set for this object.
+    surface_type: Union[SurfaceType, str]
+        The surface_type to set for this object
+    """
+
+    COORDINATE = {SurfaceType.SX: "x", SurfaceType.SY: "y", SurfaceType.SZ: "z"}
+
+    def __init__(
+        self,
+        input: InitInput = None,
+        number: int = None,
+        surface_type: Union[SurfaceType, str] = None,
+    ):
+        self._location = self._generate_default_node(float, None)
+        self._radius = self._generate_default_node(float, None)
+        super().__init__(input, number, surface_type)
+        if len(self.surface_constants) != 2:
+            raise ValueError(
+                f"{self.__class__.__name__} must have exactly 2 surface_constants"
+            )
+        self._location, self._radius = self._surface_constants
+
+    @staticmethod
+    def _number_of_params():
+        return 2
+
+    @make_prop_val_node(
+        "_radius", (float, int), float, validator=_enforce_positive_radius
+    )
+    def radius(self):
+        """The radius of the sphere
+
+        Returns
+        -------
+        float
+        """
+        pass
+
+    @make_prop_val_node("_location", (float, int), float)
+    def location(self):
+        """The location of the center of the sphere in space
+
+        Returns
+        -------
+        float
+        """
+        pass
+
+    @staticmethod
+    def _allowed_surface_types():
+        return {SurfaceType.SX, SurfaceType.SY, SurfaceType.SZ}
+
+    def validate(self):
+        super().validate()
+        if self.radius is None:
+            raise IllegalState(f"Surface: {self.number} does not have a radius set.")
+        if self.location is None:  # pragma: no cover
+            raise IllegalState(f"Surface: {self.number} does not have a location set.")
+
+    def find_duplicate_surfaces(self, surfaces, tolerance):
+        """Duplicate sphere finding is not yet implemented"""
+        return []

montepy/surfaces/surface_builder.py

@@ -5,6 +5,9 @@
 from montepy.surfaces.cylinder_on_axis import CylinderOnAxis
 from montepy.surfaces.cylinder_par_axis import CylinderParAxis
 from montepy.surfaces.general_plane import GeneralPlane
+from montepy.surfaces.general_sphere import GeneralSphere
+from montepy.surfaces.sphere_at_origin import SphereAtOrigin
+from montepy.surfaces.sphere_on_axis import SphereOnAxis
 
 
 def parse_surface(input: InitInput):
@@ -32,6 +35,12 @@ def parse_surface(input: InitInput):
         return AxisPlane(input)
     elif type_of_surface == ST.P:
         return GeneralPlane(input)
+    elif type_of_surface == ST.S:
+        return GeneralSphere(input)
+    elif type_of_surface == ST.SO:
+        return SphereAtOrigin(input)
+    elif type_of_surface in [ST.SX, ST.SY, ST.SZ]:
+        return SphereOnAxis(input)
     else:
         return buffer_surface
 

tests/inputs/test_surfaces.imcnp

@@ -1,8 +1,8 @@
 1 0  -1
 
-1 -2 SO -5
+1 SO -5
 2 PZ 0.0
-3 PZ 5.0
+3 -2 PZ 5.0
 4 1 PZ 3.0
 
 TR1 0 0 1.0

tests/test_integration.py

@@ -322,10 +322,10 @@ def test_problem_duplicate_surface_remover():
 
 def test_surface_periodic():
     problem = montepy.read_input("tests/inputs/test_surfaces.imcnp")
-    surf = problem.surfaces[1]
     periodic = problem.surfaces[2]
+    surf = problem.surfaces[3]
     assert surf.periodic_surface == periodic
-    assert "1 -2 SO" in surf.format_for_mcnp_input((6, 2, 0))[0]
+    assert "3 -2 PZ" in surf.format_for_mcnp_input((6, 2, 0))[0]
     surf.periodic_surface = problem.surfaces[3]
     assert surf.periodic_surface == problem.surfaces[3]
     del surf.periodic_surface
@@ -379,11 +379,11 @@ def test_surface_card_pass_through():
     problem = montepy.read_input("tests/inputs/test_surfaces.imcnp")
     surf = problem.surfaces[1]
     # Test input pass through
-    answer = ["1 -2 SO -5"]
+    answer = ["1 SO -5"]
     assert surf.format_for_mcnp_input((6, 2, 0)) == answer
     # Test changing periodic surface
     new_prob = copy.deepcopy(problem)
-    surf = new_prob.surfaces[1]
+    surf = new_prob.surfaces[3]
     new_prob.surfaces[2].number = 5
     assert int(surf.format_for_mcnp_input((6, 2, 0))[0].split()[1]) == -5
     # Test changing transform