diff --git a/hexrd/phase_transition/texture/__init__.py b/hexrd/phase_transition/texture/__init__.py
index 1b4c3d070..e3f046045 100644
--- a/hexrd/phase_transition/texture/__init__.py
+++ b/hexrd/phase_transition/texture/__init__.py
@@ -3,9 +3,23 @@
     SO3Kernel,
 )
 from hexrd.phase_transition.texture.uniform_odf import UniformODF
+from hexrd.phase_transition.texture.unimodal_odf import UnimodalODF
+from hexrd.phase_transition.texture.evaluation import (
+    validate_orientations,
+    eval_odf,
+    eval_odf_batch,
+    eval_at_identity,
+    eval_random_orientations,
+)
 
 __all__ = [
     'DeLaValleePoussinKernel',
     'SO3Kernel',
     'UniformODF',
+    'UnimodalODF',
+    'validate_orientations',
+    'eval_odf',
+    'eval_odf_batch',
+    'eval_at_identity',
+    'eval_random_orientations',
 ]
diff --git a/hexrd/phase_transition/texture/evaluation.py b/hexrd/phase_transition/texture/evaluation.py
new file mode 100644
index 000000000..4d19a9de4
--- /dev/null
+++ b/hexrd/phase_transition/texture/evaluation.py
@@ -0,0 +1,212 @@
+"""
+ODF Evaluation Functions
+
+Provides common evaluation functionality for orientation distribution functions.
+This module serves as a foundation for more complex ODF implementations.
+"""
+
+import numpy as np
+from scipy.spatial.transform import Rotation
+
+
+def validate_orientations(orientations):
+    """
+    Validate and normalize orientation input.
+
+    Parameters
+    ----------
+    orientations : array_like
+        Orientation matrices, expected shape (..., 3, 3)
+
+    Returns
+    -------
+    numpy.ndarray
+        Validated orientation matrices
+
+    Raises
+    ------
+    ValueError
+        If orientations don't have proper shape or aren't valid rotation matrices
+    """
+    orientations = np.asarray(orientations)
+
+    # Check basic shape requirements
+    if orientations.ndim < 2:
+        raise ValueError("Orientations must have at least 2 dimensions")
+
+    if orientations.shape[-2:] != (3, 3):
+        raise ValueError(
+            f"Orientations must have shape (..., 3, 3), got {orientations.shape}"
+        )
+
+    # Optional: Check if matrices are approximately orthogonal
+    # For performance, we'll skip this validation by default
+    # but could add a 'validate_rotation' parameter in the future
+
+    return orientations
+
+
+def eval_odf(odf, orientations, validate_input=True):
+    """
+    Generic ODF evaluation function.
+
+    Provides a unified interface for evaluating orientation distribution functions
+    at specified orientations. Supports both single orientations and batch evaluation.
+
+    Parameters
+    ----------
+    odf : ODF object
+        Any ODF object with an eval() method (e.g., UniformODF, UnimodalODF)
+    orientations : array_like
+        Orientation matrices of shape (..., 3, 3). Can be:
+        - Single 3x3 matrix
+        - Array of shape (N, 3, 3) for N orientations
+        - Higher dimensional arrays with shape (..., 3, 3)
+    validate_input : bool, optional
+        Whether to validate input orientations, default True
+
+    Returns
+    -------
+    numpy.ndarray
+        ODF values evaluated at the given orientations.
+        Shape matches the leading dimensions of input orientations.
+
+    Examples
+    --------
+    >>> from hexrd.phase_transition.texture import UniformODF, eval_odf
+    >>> odf = UniformODF('oh', 'triclinic')
+    >>>
+    >>> # Single orientation
+    >>> R = np.eye(3)
+    >>> value = eval_odf(odf, R)  # scalar result
+    >>>
+    >>> # Multiple orientations
+    >>> Rs = np.array([np.eye(3), np.eye(3)])  # shape (2, 3, 3)
+    >>> values = eval_odf(odf, Rs)  # shape (2,)
+    """
+    # Validate ODF object
+    if not hasattr(odf, 'eval'):
+        raise TypeError(
+            f"ODF object {type(odf).__name__} must have an eval() method"
+        )
+
+    # Validate and process orientations
+    if validate_input:
+        orientations = validate_orientations(orientations)
+
+    # Delegate to ODF's eval method
+    return odf.eval(orientations)
+
+
+def eval_odf_batch(odf, orientations, chunk_size=10000, validate_input=True):
+    """
+    Evaluate ODF on large batches of orientations with memory management.
+
+    For very large orientation datasets, this function processes orientations
+    in chunks to avoid memory issues while maintaining vectorization benefits.
+
+    Parameters
+    ----------
+    odf : ODF object
+        ODF object with eval() method
+    orientations : array_like
+        Large array of orientation matrices, shape (N, 3, 3)
+    chunk_size : int, optional
+        Number of orientations to process per chunk, default 10000
+    validate_input : bool, optional
+        Whether to validate input orientations, default True
+
+    Returns
+    -------
+    numpy.ndarray
+        ODF values, shape (N,)
+
+    Examples
+    --------
+    >>> odf = UniformODF('oh', 'triclinic')
+    >>> # Large batch of orientations
+    >>> Rs = np.random.normal(size=(50000, 3, 3))  # Not actual rotations!
+    >>> # In practice, use proper rotation matrices
+    >>> values = eval_odf_batch(odf, Rs, chunk_size=5000)
+    """
+    if validate_input:
+        orientations = validate_orientations(orientations)
+
+    # Handle single orientation or small batches normally
+    if orientations.ndim == 2 or orientations.shape[0] <= chunk_size:
+        return eval_odf(odf, orientations, validate_input=False)
+
+    # Process large batches in chunks
+    n_orientations = orientations.shape[0]
+    results = []
+
+    for start_idx in range(0, n_orientations, chunk_size):
+        end_idx = min(start_idx + chunk_size, n_orientations)
+        chunk = orientations[start_idx:end_idx]
+        chunk_results = eval_odf(odf, chunk, validate_input=False)
+        results.append(chunk_results)
+
+    return np.concatenate(results)
+
+
+def eval_at_identity(odf):
+    """
+    Evaluate ODF at the identity orientation.
+
+    Convenience function for getting ODF value at the identity rotation.
+    Useful for normalization checks and texture strength calculations.
+
+    Parameters
+    ----------
+    odf : ODF object
+        ODF to evaluate
+
+    Returns
+    -------
+    float
+        ODF value at identity orientation
+
+    Examples
+    --------
+    >>> odf = UniformODF('oh', 'triclinic')
+    >>> value = eval_at_identity(odf)  # 1.0 for a uniform ODF (MRD)
+    """
+    identity = np.eye(3)
+    return float(eval_odf(odf, identity, validate_input=False))
+
+
+def eval_random_orientations(odf, n_orientations=1000, seed=None):
+    """
+    Evaluate ODF at random orientations for statistical analysis.
+
+    Generates random rotation matrices and evaluates the ODF at these
+    orientations. Useful for normalization checks and Monte Carlo integration.
+
+    Parameters
+    ----------
+    odf : ODF object
+        ODF to evaluate
+    n_orientations : int, optional
+        Number of random orientations to generate, default 1000
+    seed : int, optional
+        Random seed for reproducibility
+
+    Returns
+    -------
+    tuple of (numpy.ndarray, numpy.ndarray)
+        orientations : Random rotation matrices, shape (n_orientations, 3, 3)
+        values : ODF values at these orientations, shape (n_orientations,)
+
+    Examples
+    --------
+    >>> odf = UniformODF('oh', 'triclinic')
+    >>> orientations, values = eval_random_orientations(odf, n_orientations=100)
+    >>> print(f"Mean ODF value: {np.mean(values)}")  # Should be ~1.0 (MRD)
+    """
+    # Haar-uniform random rotations on SO(3).
+    orientations = Rotation.random(n_orientations, random_state=seed).as_matrix()
+
+    # Evaluate ODF at these orientations
+    values = eval_odf(odf, orientations, validate_input=False)
+
+    return orientations, values
diff --git a/hexrd/phase_transition/texture/kernels.py b/hexrd/phase_transition/texture/kernels.py
index 826504b88..8b2391ec7 100644
--- a/hexrd/phase_transition/texture/kernels.py
+++ b/hexrd/phase_transition/texture/kernels.py
@@ -1,7 +1,7 @@
 """
 SO(3) Kernel Functions for Texture Analysis
 
-Implements radial basis functions on the rotation group SO(3) used for 
+Implements radial basis functions on the rotation group SO(3) used for
 constructing smooth orientation distribution functions.
 """
 
@@ -48,6 +48,22 @@ def _symmetry_quaternions(
         ) from error
 
 
+def _is_trivial_symmetry(quats: Optional[np.ndarray]) -> bool:
+    """
+    Return True if a symmetry group performs no reduction.
+
+    A group is trivial when it is absent (None) or contains only the
+    identity rotation (e.g. triclinic). The identity quaternion has a
+    scalar part w = ±1 with a zero vector part, while every proper,
+    non-identity rotation has |w| < 1. A group whose elements are all
+    ±identity therefore leaves misorientation angles unchanged, so it can
+    safely take the fast (no-symmetry) evaluation path.
+    """
+    if quats is None:
+        return True
+    return bool(np.allclose(np.abs(quats[0]), 1.0))
+
+
 class SO3Kernel(ABC):
     """
     Abstract base class for kernels on the SO(3) rotation group.
@@ -161,21 +177,53 @@ def __init__(
 
         # Normalization constant: C = B(3/2, 1/2) / B(3/2, κ + 1/2)
         self._C = betafn(1.5, 0.5) / betafn(1.5, self._kappa + 0.5)
-        self._crystal_symmetry = _symmetry_quaternions(
+        self._crystal_symmetry_quats = _symmetry_quaternions(
             crystal_symmetry,
             symtype='crystal',
         )
-        self._sample_symmetry = _symmetry_quaternions(
+        self._sample_symmetry_quats = _symmetry_quaternions(
             sample_symmetry,
             symtype='sample',
         )
 
+        # Retain the original symmetry labels (when given as strings) so the
+        # kernel can act as the single source of truth for symmetry. When
+        # symmetry is supplied as a quaternion array, no label is available.
+        self._crystal_symmetry_label = (
+            crystal_symmetry if isinstance(crystal_symmetry, str) else None
+        )
+        self._sample_symmetry_label = (
+            sample_symmetry if isinstance(sample_symmetry, str) else None
+        )
+
+    @property
+    def crystal_symmetry(self) -> Optional[str]:
+        """
+        str or None: Crystal symmetry label, if built from a label.
+
+        Returns None when the symmetry was supplied directly as a
+        quaternion array, even if symmetry reduction is active. Use
+        ``has_symmetry`` to test whether reduction is enabled.
+        """
+        return self._crystal_symmetry_label
+
+    @property
+    def sample_symmetry(self) -> Optional[str]:
+        """
+        str or None: Sample symmetry label, if built from a label.
+
+        Returns None when the symmetry was supplied directly as a
+        quaternion array, even if symmetry reduction is active. Use
+        ``has_symmetry`` to test whether reduction is enabled.
+        """
+        return self._sample_symmetry_label
+
     @property
     def has_symmetry(self) -> bool:
-        """bool: Whether symmetry reduction is enabled."""
-        return (
-            self._crystal_symmetry is not None
-            or self._sample_symmetry is not None
+        """bool: Whether non-trivial symmetry reduction is enabled."""
+        return not (
+            _is_trivial_symmetry(self._crystal_symmetry_quats)
+            and _is_trivial_symmetry(self._sample_symmetry_quats)
         )
 
     @property
@@ -251,6 +299,12 @@ def _symmetry_reduced_misorientation_angle(
 
         angles = np.empty(R1_flat.shape[0], dtype=float)
         symmetries = self._misorientation_symmetries()
+        # TODO(perf): this per-pair Python loop calls rotations.misorientation
+        # once per orientation and is ~1400x slower than the no-symmetry path
+        # for large grids (e.g. ~2 s for 3000 cubic-symmetry evaluations).
+        # Vectorize by applying the symmetry quaternion group with broadcasting
+        # and taking the minimum angle. Revisit on the texture-normalization
+        # branch before evaluating ODFs on dense (1e5-1e6) grids.
         for i, (r1, r2) in enumerate(zip(R1_flat, R2_flat)):
             q1 = rotations.quatOfRotMat(r1).reshape(4, 1)
             q2 = rotations.quatOfRotMat(r2).reshape(4, 1)
@@ -266,12 +320,12 @@ def _symmetry_reduced_misorientation_angle(
         return angles.reshape(output_shape)
 
     def _misorientation_symmetries(self) -> tuple[np.ndarray, ...]:
-        crystal_symmetry = self._crystal_symmetry
-        if crystal_symmetry is None and self._sample_symmetry is not None:
+        crystal_symmetry = self._crystal_symmetry_quats
+        if crystal_symmetry is None and self._sample_symmetry_quats is not None:
             crystal_symmetry = _IDENTITY_SYMMETRY
         symmetries = (crystal_symmetry,) if crystal_symmetry is not None else ()
-        if self._sample_symmetry is not None:
-            symmetries = symmetries + (self._sample_symmetry,)
+        if self._sample_symmetry_quats is not None:
+            symmetries = symmetries + (self._sample_symmetry_quats,)
         return symmetries
 
     def eval(
diff --git a/hexrd/phase_transition/texture/unimodal_odf.py b/hexrd/phase_transition/texture/unimodal_odf.py
new file mode 100644
index 000000000..abb7613ca
--- /dev/null
+++ b/hexrd/phase_transition/texture/unimodal_odf.py
@@ -0,0 +1,286 @@
+"""
+Unimodal Orientation Distribution Function (ODF)
+
+Implements a kernel-based ODF with a single or multiple preferred orientations.
+Uses a combination of modal orientations and kernel functions to create
+smooth, localized texture distributions.
+"""
+
+import numpy as np
+from .kernels import DeLaValleePoussinKernel
+
+
+class UnimodalODF:
+    """
+    Unimodal orientation distribution function.
+
+    Represents a texture with one or more preferred orientations around which
+    crystal orientations are concentrated. Uses kernel functions to create
+    smooth distributions around modal orientations.
+
+    Mathematical form: f(g) = Σᵢ wᵢ K(g, g₀ᵢ),  Σᵢ wᵢ = 1
+    Because the de la Vallée Poussin kernel is normalized so its mean over
+    SO(3) is 1 (MRD), this weighted sum is itself a valid ODF in MRD units
+    (uniform texture = 1).
+
+    Parameters
+    ----------
+    modal_orientations : array_like
+        Modal (preferred) orientation(s). Can be:
+        - Single 3x3 rotation matrix
+        - Array of shape (N, 3, 3) for N modal orientations
+    kernel : DeLaValleePoussinKernel
+        Kernel function defining the shape of the distribution. The kernel
+        is the single source of truth for symmetry: any crystal/sample
+        symmetry must be set on the kernel, and the ODF exposes it via the
+        ``crystal_symmetry``/``sample_symmetry`` properties.
+    weights : array_like, optional
+        Weights for multiple modal orientations. Must sum to 1.
+        If None, equal weights are used for multiple orientations.
+
+    Attributes
+    ----------
+    modal_orientations : numpy.ndarray
+        Array of modal orientations, shape (N, 3, 3)
+    kernel : DeLaValleePoussinKernel
+        Kernel function used for the distribution
+    weights : numpy.ndarray
+        Component weights, shape (N,)
+    crystal_symmetry : str or None
+        Crystal symmetry label, delegated from the kernel
+    sample_symmetry : str or None
+        Sample symmetry label, delegated from the kernel
+    n_components : int
+        Number of modal orientations
+
+    Examples
+    --------
+    >>> from hexrd.phase_transition.texture import UnimodalODF
+    >>> from hexrd.phase_transition.texture import DeLaValleePoussinKernel
+    >>>
+    >>> # Single modal orientation (cubic crystal symmetry on the kernel)
+    >>> kernel = DeLaValleePoussinKernel(
+    ...     halfwidth=np.radians(15), crystal_symmetry='oh'
+    ... )
+    >>> modal = np.eye(3)  # Identity orientation
+    >>> odf = UnimodalODF(modal, kernel)
+    >>>
+    >>> # Evaluate at modal orientation (should give maximum value)
+    >>> value = odf.eval(modal)
+    """
+
+    def __init__(self, modal_orientations, kernel, weights=None):
+        """
+        Initialize unimodal ODF.
+
+        Parameters
+        ----------
+        modal_orientations : array_like
+            Modal orientation(s) as rotation matrices
+        kernel : DeLaValleePoussinKernel
+            Kernel function for the distribution. Symmetry, if any, must be
+            configured on the kernel; the ODF delegates to it.
+        weights : array_like, optional
+            Component weights for multiple modal orientations
+        """
+        # Validate and store kernel. The kernel is the single source of
+        # truth for symmetry (validated when the kernel was constructed).
+        if not isinstance(kernel, DeLaValleePoussinKernel):
+            raise TypeError("kernel must be a DeLaValleePoussinKernel instance")
+
+        self._kernel = kernel
+
+        # Process modal orientations
+        modal_orientations = np.asarray(modal_orientations)
+
+        # Handle single vs multiple modal orientations
+        if modal_orientations.ndim == 2 and modal_orientations.shape == (3, 3):
+            # Single modal orientation
+            self._modal_orientations = modal_orientations.reshape(1, 3, 3)
+            self._n_components = 1
+        elif modal_orientations.ndim == 3 and modal_orientations.shape[-2:] == (3, 3):
+            # Multiple modal orientations
+            self._modal_orientations = modal_orientations
+            self._n_components = modal_orientations.shape[0]
+        else:
+            raise ValueError(
+                f"Modal orientations must have shape (3, 3) or (N, 3, 3), "
+                f"got {modal_orientations.shape}"
+            )
+
+        # Process weights
+        if weights is None:
+            # Equal weights for all components
+            self._weights = np.full(self._n_components, 1.0 / self._n_components)
+        else:
+            weights = np.asarray(weights)
+            if weights.shape != (self._n_components,):
+                raise ValueError(
+                    f"Weights must have shape ({self._n_components},), "
+                    f"got {weights.shape}"
+                )
+            if not np.allclose(np.sum(weights), 1.0):
+                raise ValueError("Weights must sum to 1.0")
+            if np.any(weights < 0):
+                raise ValueError("Weights must be non-negative")
+            self._weights = weights.copy()
+
+    @property
+    def modal_orientations(self):
+        """numpy.ndarray: Modal orientations, shape (N, 3, 3)."""
+        return self._modal_orientations.copy()
+
+    @property
+    def kernel(self):
+        """DeLaValleePoussinKernel: Kernel function."""
+        return self._kernel
+
+    @property
+    def weights(self):
+        """numpy.ndarray: Component weights, shape (N,)."""
+        return self._weights.copy()
+
+    @property
+    def crystal_symmetry(self):
+        """str or None: Crystal symmetry label, delegated from the kernel."""
+        return self._kernel.crystal_symmetry
+
+    @property
+    def sample_symmetry(self):
+        """str or None: Sample symmetry label, delegated from the kernel."""
+        return self._kernel.sample_symmetry
+
+    @property
+    def n_components(self):
+        """int: Number of modal orientations."""
+        return self._n_components
+
+    def eval(self, orientations):
+        """
+        Evaluate unimodal ODF at given orientations.
+
+        Computes f(g) = Σᵢ wᵢ K(g, g₀ᵢ) for all input orientations, in MRD.
+
+        Parameters
+        ----------
+        orientations : array_like
+            Orientation matrices of shape (..., 3, 3)
+
+        Returns
+        -------
+        numpy.ndarray
+            ODF values with shape matching leading dimensions of input
+
+        Examples
+        --------
+        >>> modal = np.eye(3)
+        >>> kernel = DeLaValleePoussinKernel(halfwidth=np.radians(10))
+        >>> odf = UnimodalODF(modal, kernel)
+        >>>
+        >>> # Single evaluation
+        >>> value = odf.eval(modal)  # Should give maximum
+        >>>
+        >>> # Batch evaluation
+        >>> Rs = np.array([np.eye(3), rotation_matrix_z(np.pi/4)])
+        >>> values = odf.eval(Rs)  # shape (2,)
+        """
+        orientations = np.asarray(orientations)
+
+        # Validate input shape
+        if orientations.shape[-2:] != (3, 3):
+            raise ValueError(
+                f"Orientation matrices must have shape (..., 3, 3), "
+                f"got {orientations.shape}"
+            )
+
+        # Determine output shape
+        output_shape = orientations.shape[:-2]
+        if output_shape == ():
+            # Single orientation
+            output_shape = ()
+            orientations = orientations.reshape(1, 3, 3)
+            squeeze_output = True
+        else:
+            # Multiple orientations - flatten for processing
+            n_orientations = int(np.prod(output_shape))
+            orientations = orientations.reshape(n_orientations, 3, 3)
+            squeeze_output = False
+
+        # Initialize results
+        n_query = orientations.shape[0]
+        results = np.zeros(n_query)
+
+        # Evaluate each component and sum with weights
+        for i in range(self.n_components):
+            modal_i = self._modal_orientations[i]  # Shape (3, 3)
+            weight_i = self._weights[i]
+
+            # Evaluate kernel between all query orientations and this modal orientation
+            # We need to broadcast: kernel.eval expects both args to have same leading dims
+            modal_broadcast = np.broadcast_to(
+                modal_i.reshape(1, 3, 3),
+                (n_query, 3, 3)
+            )
+
+            kernel_values = self.kernel.eval(orientations, modal_broadcast)
+
+            # Add weighted contribution to results
+            results += weight_i * kernel_values
+
+        # The de la Vallée Poussin kernel is normalized so its mean over
+        # SO(3) is 1 (MRD). With weights that sum to 1, this weighted sum is
+        # already a valid ODF in MRD units (uniform texture = 1), so no
+        # additional normalization is required.
+
+        # Reshape to match input
+        if squeeze_output:
+            return float(results[0])
+        else:
+            return results.reshape(output_shape)
+
+    def estimated_max_value(self):
+        """
+        Estimate the maximum ODF value, in MRD.
+
+        The ODF maxima occur at (or very near) the modal orientations, so
+        this evaluates the full ODF at each mode and returns the largest
+        value.
+
+        Returns
+        -------
+        float
+            Maximum ODF value in MRD (multiples of a random distribution)
+        """
+        # ODF maxima occur at (or very near) the modal orientations.
+        # Evaluating the full ODF at each mode accounts for overlap between
+        # nearby modes.
+        modal_values = np.atleast_1d(self.eval(self._modal_orientations))
+        return float(np.max(modal_values))
+
+    def __repr__(self):
+        """String representation of UnimodalODF."""
+        return (
+            f"UnimodalODF(n_components={self.n_components}, "
+            f"kernel_halfwidth={np.degrees(self.kernel.halfwidth):.1f}°, "
+            f"crystal_symmetry={self.crystal_symmetry!r}, "
+            f"sample_symmetry={self.sample_symmetry!r})"
+        )
+
+    def __str__(self):
+        """Human-readable description."""
+        desc = (
+            f"Unimodal ODF with {self.n_components} component(s)\n"
+            f"Crystal symmetry: {self.crystal_symmetry or 'none'}\n"
+            f"Sample symmetry: {self.sample_symmetry or 'none'}\n"
+            f"Kernel: {self.kernel}\n"
+        )
+
+        if self.n_components == 1:
+            desc += "Modal orientations: 1 component\n"
+        else:
+            desc += f"Modal orientations: {self.n_components} components\n"
+            desc += f"Weights: {self.weights}\n"
+
+        desc += f"Estimated max value: {self.estimated_max_value():.1f} MRD"
+
+        return desc
diff --git a/tests/test_delavalleepoussin_kernel.py b/tests/test_delavalleepoussin_kernel.py
index 894123f74..3b7c736a6 100644
--- a/tests/test_delavalleepoussin_kernel.py
+++ b/tests/test_delavalleepoussin_kernel.py
@@ -88,7 +88,7 @@ def test_kappa_from_halfwidth(self):
         self.assertAlmostEqual(
             kernel.kappa, expected_kappa, places=10
         )
-    
+
     def test_kernel_evaluation_at_identity(self):
         """Test that K(0) equals the normalization constant.
 
@@ -123,7 +123,7 @@ def test_string_representations(self):
         self.assertIn('de la Vallée Poussin kernel', str_repr)
         self.assertIn('half-width', str_repr)
         self.assertIn('K(ω)', str_repr)
-        
+
     def test_kernel_decreases_with_misorientation(self):
         """Test that kernel value decreases as misorientation grows."""
         halfwidth_rad = np.radians(4.0)
@@ -361,3 +361,50 @@ def test_batch_evaluation(self):
         for val in values[1:]:
             self.assertAlmostEqual(float(val), 0.0, places=5)
 
+    # --- symmetry label properties ---
+
+    def test_symmetry_labels_retained(self):
+        """String symmetry labels are exposed via read-only properties."""
+        kernel = DeLaValleePoussinKernel(
+            halfwidth=np.radians(10),
+            crystal_symmetry='d4h',
+            sample_symmetry='orthorhombic',
+        )
+        self.assertEqual(kernel.crystal_symmetry, 'd4h')
+        self.assertEqual(kernel.sample_symmetry, 'orthorhombic')
+        self.assertTrue(kernel.has_symmetry)
+
+    def test_symmetry_labels_none_without_symmetry(self):
+        """A kernel built without symmetry reports None labels."""
+        kernel = DeLaValleePoussinKernel(halfwidth=np.radians(10))
+        self.assertIsNone(kernel.crystal_symmetry)
+        self.assertIsNone(kernel.sample_symmetry)
+        self.assertFalse(kernel.has_symmetry)
+
+    def test_symmetry_label_none_for_array_input(self):
+        """Symmetry given as a quaternion array has no recoverable label."""
+        kernel = DeLaValleePoussinKernel(
+            halfwidth=np.radians(10),
+            crystal_symmetry=rotations.quatOfLaueGroup('d4h'),
+        )
+        self.assertIsNone(kernel.crystal_symmetry)
+        self.assertTrue(kernel.has_symmetry)
+
+    def test_trivial_symmetry_takes_fast_path(self):
+        """Triclinic symmetry is the identity, so it enables no reduction."""
+        kernel = DeLaValleePoussinKernel(
+            halfwidth=np.radians(10),
+            crystal_symmetry='ci',
+            sample_symmetry='triclinic',
+        )
+        # Labels are still retained for the single-source-of-truth contract.
+        self.assertEqual(kernel.crystal_symmetry, 'ci')
+        self.assertEqual(kernel.sample_symmetry, 'triclinic')
+        # But the group is trivial, so the fast (no-reduction) path is used.
+        self.assertFalse(kernel.has_symmetry)
+
+        # A 90° misorientation is unchanged by the trivial symmetry group.
+        r90z = _rotation_about_z(np.radians(90))
+        angle = kernel.misorientation_angle(np.eye(3), r90z)
+        self.assertAlmostEqual(float(angle), np.radians(90), places=8)
+
diff --git a/tests/test_evaluation.py b/tests/test_evaluation.py
new file mode 100644
index 000000000..ecf8ddbc8
--- /dev/null
+++ b/tests/test_evaluation.py
@@ -0,0 +1,145 @@
+"""Tests for the generic ODF evaluation helpers."""
+
+import numpy as np
+import unittest
+from scipy.spatial.transform import Rotation
+
+from hexrd.phase_transition.texture import (
+    UniformODF,
+    UnimodalODF,
+    DeLaValleePoussinKernel,
+    validate_orientations,
+    eval_odf,
+    eval_odf_batch,
+    eval_at_identity,
+    eval_random_orientations,
+)
+
+
+class TestValidateOrientations(unittest.TestCase):
+    """Test validate_orientations."""
+
+    def test_valid_single(self):
+        result = validate_orientations(np.eye(3))
+        self.assertEqual(result.shape, (3, 3))
+
+    def test_valid_batch(self):
+        batch = np.tile(np.eye(3), (5, 1, 1))
+        result = validate_orientations(batch)
+        self.assertEqual(result.shape, (5, 3, 3))
+
+    def test_too_few_dimensions(self):
+        with self.assertRaises(ValueError):
+            validate_orientations(np.array([1.0, 2.0, 3.0]))
+
+    def test_wrong_trailing_shape(self):
+        with self.assertRaises(ValueError):
+            validate_orientations(np.zeros((4, 3, 2)))
+
+
+class TestEvalOdf(unittest.TestCase):
+    """Test eval_odf against the real ODF classes."""
+
+    def setUp(self):
+        self.uniform = UniformODF('oh', 'triclinic')
+        kernel = DeLaValleePoussinKernel(halfwidth=np.radians(15))
+        self.unimodal = UnimodalODF(np.eye(3), kernel)
+
+    def test_single_uniform_is_one_mrd(self):
+        self.assertEqual(eval_odf(self.uniform, np.eye(3)), 1.0)
+
+    def test_batch_uniform_is_ones(self):
+        batch = np.tile(np.eye(3), (4, 1, 1))
+        np.testing.assert_array_equal(eval_odf(self.uniform, batch), np.ones(4))
+
+    def test_matches_odf_eval(self):
+        orientations = Rotation.random(50, random_state=1).as_matrix()
+        np.testing.assert_allclose(
+            eval_odf(self.unimodal, orientations),
+            self.unimodal.eval(orientations),
+        )
+
+    def test_requires_eval_method(self):
+        with self.assertRaises(TypeError):
+            eval_odf(object(), np.eye(3))
+
+
+class TestEvalOdfBatch(unittest.TestCase):
+    """Test eval_odf_batch chunking."""
+
+    def setUp(self):
+        self.uniform = UniformODF('oh', 'triclinic')
+        kernel = DeLaValleePoussinKernel(halfwidth=np.radians(15))
+        self.unimodal = UnimodalODF(np.eye(3), kernel)
+
+    def test_small_batch_matches_eval_odf(self):
+        batch = np.tile(np.eye(3), (3, 1, 1))
+        np.testing.assert_array_equal(
+            eval_odf_batch(self.uniform, batch),
+            eval_odf(self.uniform, batch),
+        )
+
+    def test_chunking_matches_direct_eval(self):
+        orientations = Rotation.random(2500, random_state=2).as_matrix()
+        chunked = eval_odf_batch(self.unimodal, orientations, chunk_size=1000)
+        self.assertEqual(chunked.shape, (2500,))
+        np.testing.assert_allclose(chunked, self.unimodal.eval(orientations))
+
+    def test_chunking_preserves_uniform_value(self):
+        batch = np.tile(np.eye(3), (2500, 1, 1))
+        chunked = eval_odf_batch(self.uniform, batch, chunk_size=500)
+        np.testing.assert_array_equal(chunked, np.ones(2500))
+
+
+class TestEvalAtIdentity(unittest.TestCase):
+    """Test eval_at_identity (and MRD convention)."""
+
+    def test_uniform_identity_is_one_mrd(self):
+        odf = UniformODF('oh', 'triclinic')
+        self.assertEqual(eval_at_identity(odf), 1.0)
+
+    def test_unimodal_identity_is_kernel_peak(self):
+        kernel = DeLaValleePoussinKernel(halfwidth=np.radians(15))
+        odf = UnimodalODF(np.eye(3), kernel)
+        self.assertAlmostEqual(
+            eval_at_identity(odf), kernel.norm_constant, places=6
+        )
+
+
+class TestEvalRandomOrientations(unittest.TestCase):
+    """Test eval_random_orientations."""
+
+    def setUp(self):
+        self.uniform = UniformODF('oh', 'triclinic')
+
+    def test_output_shapes(self):
+        orientations, values = eval_random_orientations(
+            self.uniform, n_orientations=200, seed=0
+        )
+        self.assertEqual(orientations.shape, (200, 3, 3))
+        self.assertEqual(values.shape, (200,))
+
+    def test_returns_valid_rotations(self):
+        orientations, _ = eval_random_orientations(
+            self.uniform, n_orientations=50, seed=0
+        )
+        identity = np.tile(np.eye(3), (50, 1, 1))
+        products = np.matmul(orientations, np.swapaxes(orientations, -2, -1))
+        np.testing.assert_allclose(products, identity, atol=1e-10)
+        dets = np.linalg.det(orientations)
+        np.testing.assert_allclose(dets, np.ones(50), atol=1e-10)
+
+    def test_reproducible_with_seed(self):
+        o1, _ = eval_random_orientations(self.uniform, n_orientations=20, seed=42)
+        o2, _ = eval_random_orientations(self.uniform, n_orientations=20, seed=42)
+        np.testing.assert_array_equal(o1, o2)
+
+    def test_uniform_mean_is_one_mrd(self):
+        _, values = eval_random_orientations(
+            self.uniform, n_orientations=1000, seed=3
+        )
+        self.assertAlmostEqual(np.mean(values), 1.0)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_unimodal_odf.py b/tests/test_unimodal_odf.py
new file mode 100644
index 000000000..65fb67cf1
--- /dev/null
+++ b/tests/test_unimodal_odf.py
@@ -0,0 +1,247 @@
+"""Tests for UnimodalODF."""
+
+import numpy as np
+import unittest
+from scipy.spatial.transform import Rotation
+
+from hexrd.phase_transition.texture.kernels import DeLaValleePoussinKernel
+from hexrd.phase_transition.texture.unimodal_odf import UnimodalODF
+
+
+def _rotation_about_z(angle: float) -> np.ndarray:
+    cos_angle = np.cos(angle)
+    sin_angle = np.sin(angle)
+    return np.array([
+        [cos_angle, -sin_angle, 0.0],
+        [sin_angle, cos_angle, 0.0],
+        [0.0, 0.0, 1.0],
+    ])
+
+
+def _cubic_kernel(halfwidth_deg: float = 15.0) -> DeLaValleePoussinKernel:
+    """A kernel carrying cubic crystal symmetry (single source of truth)."""
+    return DeLaValleePoussinKernel(
+        halfwidth=np.radians(halfwidth_deg),
+        crystal_symmetry='oh',
+        sample_symmetry='triclinic',
+    )
+
+
+class TestUnimodalODF(unittest.TestCase):
+    """Test the real UnimodalODF against the real DeLaValleePoussinKernel."""
+
+    def setUp(self):
+        # Default kernel carries NO symmetry, so modes stay distinct.
+        self.kernel = DeLaValleePoussinKernel(halfwidth=np.radians(15))
+        self.modal = np.eye(3)
+
+    # --- construction ---
+
+    def test_package_export(self):
+        """UnimodalODF is importable from the texture package."""
+        from hexrd.phase_transition.texture import UnimodalODF as Cls
+        self.assertIs(Cls, UnimodalODF)
+
+    def test_initialization_single_modal(self):
+        """A single (3, 3) modal orientation is stored as one component."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        self.assertEqual(odf.n_components, 1)
+        self.assertEqual(odf.modal_orientations.shape, (1, 3, 3))
+        np.testing.assert_array_equal(odf.modal_orientations[0], self.modal)
+        np.testing.assert_array_almost_equal(odf.weights, [1.0])
+        self.assertIs(odf.kernel, self.kernel)
+
+    def test_initialization_multiple_modals(self):
+        """Multiple modal orientations and explicit weights are stored."""
+        modals = np.array([np.eye(3), _rotation_about_z(np.radians(90))])
+        odf = UnimodalODF(modals, self.kernel, weights=[0.6, 0.4])
+        self.assertEqual(odf.n_components, 2)
+        np.testing.assert_array_equal(odf.modal_orientations, modals)
+        np.testing.assert_array_almost_equal(odf.weights, [0.6, 0.4])
+
+    def test_equal_weights_default(self):
+        """Equal weights are assigned when none are provided."""
+        modals = np.array([
+            np.eye(3),
+            _rotation_about_z(np.radians(45)),
+            _rotation_about_z(np.radians(90)),
+        ])
+        odf = UnimodalODF(modals, self.kernel)
+        np.testing.assert_array_almost_equal(odf.weights, np.full(3, 1.0 / 3))
+
+    # --- symmetry delegation (Option A: kernel is the single source) ---
+
+    def test_symmetry_delegates_to_kernel(self):
+        """crystal/sample symmetry are read straight from the kernel."""
+        kernel = _cubic_kernel()
+        odf = UnimodalODF(self.modal, kernel)
+        self.assertEqual(odf.crystal_symmetry, 'oh')
+        self.assertEqual(odf.sample_symmetry, 'triclinic')
+        self.assertEqual(odf.crystal_symmetry, kernel.crystal_symmetry)
+        self.assertEqual(odf.sample_symmetry, kernel.sample_symmetry)
+
+    def test_no_symmetry_kernel_reports_none(self):
+        """A kernel built without symmetry reports None for both labels."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        self.assertIsNone(odf.crystal_symmetry)
+        self.assertIsNone(odf.sample_symmetry)
+
+    def test_symmetry_makes_equivalent_modes_equal(self):
+        """Symmetry now flows through the kernel into eval.
+
+        Under cubic symmetry, a 90 deg rotation about z is equivalent to the
+        identity, so the ODF evaluates equally at both. Without symmetry the
+        two orientations are distinct.
+        """
+        cubic = UnimodalODF(self.modal, _cubic_kernel())
+        nosym = UnimodalODF(self.modal, self.kernel)
+        r90 = _rotation_about_z(np.radians(90))
+        self.assertAlmostEqual(cubic.eval(r90), cubic.eval(np.eye(3)), places=6)
+        self.assertLess(nosym.eval(r90), nosym.eval(np.eye(3)))
+
+    # --- validation ---
+
+    def test_invalid_kernel_type(self):
+        """A kernel that is not a DeLaValleePoussinKernel raises TypeError."""
+        with self.assertRaises(TypeError):
+            UnimodalODF(self.modal, object())
+
+    def test_invalid_modal_orientation_shape(self):
+        """Modal orientations must have shape (3, 3) or (N, 3, 3)."""
+        with self.assertRaises(ValueError):
+            UnimodalODF(np.array([1, 2, 3]), self.kernel)
+        with self.assertRaises(ValueError):
+            UnimodalODF(np.random.random((3, 2)), self.kernel)
+
+    def test_invalid_weights(self):
+        """Weights must match component count, be non-negative, and sum to 1."""
+        modals = np.array([np.eye(3), _rotation_about_z(np.radians(90))])
+        with self.assertRaises(ValueError):
+            UnimodalODF(modals, self.kernel, weights=[1.0])
+        with self.assertRaises(ValueError):
+            UnimodalODF(modals, self.kernel, weights=[0.3, 0.3])
+        with self.assertRaises(ValueError):
+            UnimodalODF(modals, self.kernel, weights=[1.2, -0.2])
+
+    # --- evaluation ---
+
+    def test_single_orientation_evaluation(self):
+        """A single orientation evaluates to a positive scalar, peaked at the mode."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        value_modal = odf.eval(self.modal)
+        value_far = odf.eval(_rotation_about_z(np.radians(90)))
+        self.assertIsInstance(value_modal, float)
+        self.assertIsInstance(value_far, float)
+        self.assertGreater(value_modal, value_far)
+        self.assertGreater(value_far, 0.0)
+
+    def test_eval_at_mode_equals_kernel_peak(self):
+        """For a single mode (weight 1), eval at the mode equals the kernel peak C."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        self.assertAlmostEqual(
+            odf.eval(self.modal), self.kernel.norm_constant, places=6
+        )
+
+    def test_multiple_orientations_evaluation(self):
+        """A batch of orientations returns a 1-D array peaked at the mode."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        orientations = np.array([
+            np.eye(3),
+            _rotation_about_z(np.radians(45)),
+            _rotation_about_z(np.radians(90)),
+        ])
+        values = odf.eval(orientations)
+        self.assertEqual(values.shape, (3,))
+        self.assertEqual(np.argmax(values), 0)
+        self.assertTrue(np.all(values > 0))
+
+    def test_batch_orientation_evaluation(self):
+        """Multi-dimensional batches preserve leading dimensions."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        batch = np.tile(self.modal, (2, 3, 1, 1))  # (2, 3, 3, 3)
+        values = odf.eval(batch)
+        self.assertEqual(values.shape, (2, 3))
+        np.testing.assert_allclose(values, self.kernel.norm_constant)
+
+    def test_invalid_evaluation_shapes(self):
+        """Evaluation rejects inputs that are not (..., 3, 3)."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        with self.assertRaises(ValueError):
+            odf.eval(np.array([1, 2, 3]))
+        with self.assertRaises(ValueError):
+            odf.eval(np.random.random((3, 2)))
+
+    def test_weighted_components(self):
+        """At a well-separated mode the value is approximately w_i * C.
+
+        The default kernel carries no crystal symmetry, so the two modes
+        (identity and 90 deg about z) are treated as distinct orientations.
+        """
+        modals = np.array([np.eye(3), _rotation_about_z(np.radians(90))])
+        odf = UnimodalODF(modals, self.kernel, weights=[0.7, 0.3])
+        peak = self.kernel.norm_constant
+        value_0 = odf.eval(np.eye(3))
+        value_1 = odf.eval(_rotation_about_z(np.radians(90)))
+        self.assertGreater(value_0, value_1)
+        np.testing.assert_allclose(value_0, 0.7 * peak, rtol=1e-2)
+        np.testing.assert_allclose(value_1, 0.3 * peak, rtol=1e-2)
+
+    def test_mrd_normalization(self):
+        """The ODF integrates to mean 1 (MRD) over SO(3).
+
+        A weighted de la Vallee Poussin kernel sum is already normalized,
+        so the Monte Carlo mean over uniformly random orientations is ~1.
+        """
+        kernel = DeLaValleePoussinKernel(halfwidth=np.radians(40))
+        odf = UnimodalODF(np.eye(3), kernel)
+        samples = Rotation.random(60000, random_state=7).as_matrix()
+        mean_mrd = np.mean(odf.eval(samples))
+        self.assertAlmostEqual(mean_mrd, 1.0, delta=0.03)
+
+    # --- derived quantities ---
+
+    def test_estimated_max_value(self):
+        """estimated_max_value equals the peak (kernel C) for a single mode."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        max_val = odf.estimated_max_value()
+        self.assertIsInstance(max_val, float)
+        self.assertAlmostEqual(max_val, self.kernel.norm_constant, places=6)
+        self.assertAlmostEqual(max_val, odf.eval(self.modal), places=6)
+
+    # --- properties / representations ---
+
+    def test_property_getters_return_copies(self):
+        """modal_orientations and weights getters return independent copies."""
+        modals = np.array([np.eye(3), _rotation_about_z(np.radians(90))])
+        odf = UnimodalODF(modals, self.kernel, weights=[0.6, 0.4])
+
+        modals_view = odf.modal_orientations
+        modals_view[0] *= 2
+        np.testing.assert_array_equal(odf.modal_orientations[0], np.eye(3))
+
+        weights_view = odf.weights
+        weights_view[0] = 0.9
+        self.assertAlmostEqual(odf.weights[0], 0.6)
+
+    def test_string_representations(self):
+        """__repr__ and __str__ contain the key descriptors."""
+        odf = UnimodalODF(self.modal, _cubic_kernel())
+
+        repr_str = repr(odf)
+        self.assertIn('UnimodalODF', repr_str)
+        self.assertIn('oh', repr_str)
+        self.assertIn('triclinic', repr_str)
+
+        str_repr = str(odf)
+        self.assertIn('Unimodal ODF', str_repr)
+        self.assertIn('MRD', str_repr)
+
+    def test_string_representations_without_symmetry(self):
+        """repr/str render gracefully when the kernel has no symmetry."""
+        odf = UnimodalODF(self.modal, self.kernel)
+        self.assertIn('crystal_symmetry=None', repr(odf))
+        self.assertIn('none', str(odf))
+
+
+if __name__ == '__main__':
+    unittest.main()