diff --git a/pyroomacoustics/directivities.py b/pyroomacoustics/directivities.py
index 8678c7d3..97149d96 100644
--- a/pyroomacoustics/directivities.py
+++ b/pyroomacoustics/directivities.py
@@ -331,54 +331,3 @@ def cardioid_func(x, direction, coef, gain=1.0, normalize=True, magnitude=False)
         return np.abs(resp)
     else:
         return resp
-
-
-def source_angle_shoebox(image_source_loc, wall_flips, mic_loc):
-    """
-    Determine outgoing angle for each image source for a ShoeBox configuration.
-
-    Implementation of the method described in the paper:
-    https://www2.ak.tu-berlin.de/~akgroup/ak_pub/2018/000458.pdf
-
-    Parameters
-    -----------
-    image_source_loc : array_like
-        Locations of image sources.
-    wall_flips: array_like
-        Number of x, y, z flips for each image source.
-    mic_loc: array_like
-        Microphone location.
-
-    Returns
-    -------
-    azimuth : :py:class:`~numpy.ndarray`
-        Azimith for each image source, in radians
-    colatitude : :py:class:`~numpy.ndarray`
-        Colatitude for each image source, in radians.
-
-    """
-
-    image_source_loc = np.array(image_source_loc)
-    wall_flips = np.array(wall_flips)
-    mic_loc = np.array(mic_loc)
-
-    dim, n_sources = image_source_loc.shape
-    assert wall_flips.shape[0] == dim
-    assert mic_loc.shape[0] == dim
-
-    p_vector_array = image_source_loc - np.array(mic_loc)[:, np.newaxis]
-    d_array = np.linalg.norm(p_vector_array, axis=0)
-
-    # Using (12) from the paper
-    power_array = np.ones_like(image_source_loc) * -1
-    power_array = np.power(power_array, (wall_flips + np.ones_like(image_source_loc)))
-    p_dash_array = p_vector_array * power_array
-
-    # Using (13) from the paper
-    azimuth = np.arctan2(p_dash_array[1], p_dash_array[0])
-    if dim == 2:
-        colatitude = np.ones(n_sources) * np.pi / 2
-    else:
-        colatitude = np.pi / 2 - np.arcsin(p_dash_array[2] / d_array)
-
-    return azimuth, colatitude
diff --git a/pyroomacoustics/room.py b/pyroomacoustics/room.py
index f01f1e8a..209f6f4f 100644
--- a/pyroomacoustics/room.py
+++ b/pyroomacoustics/room.py
@@ -698,7 +698,7 @@ def callback_mix(premix, snr=0, sir=0, ref_mic=0, n_src=None, n_tgt=None):
 from . import libroom
 from .acoustics import OctaveBandsFactory, rt60_eyring, rt60_sabine
 from .beamforming import MicrophoneArray
-from .directivities import CardioidFamily, source_angle_shoebox
+from .directivities import CardioidFamily
 from .experimental import measure_rt60
 from .libroom import Wall, Wall2D
 from .open_sofa_interpolate import MeasuredDirectivity
diff --git a/pyroomacoustics/simulation/ism.py b/pyroomacoustics/simulation/ism.py
index 1412dbca..60e7b830 100644
--- a/pyroomacoustics/simulation/ism.py
+++ b/pyroomacoustics/simulation/ism.py
@@ -4,7 +4,6 @@
 from scipy.signal import fftconvolve, hilbert
 
 from .. import libroom
-from ..directivities import source_angle_shoebox
 from ..parameters import constants
 from ..utilities import angle_function
 
@@ -82,6 +81,57 @@ def interpolate_octave_bands(
     return ir
 
 
+def source_angle_shoebox(image_source_loc, wall_flips, mic_loc):
+    """
+    Determine outgoing angle for each image source for a ShoeBox configuration.
+
+    Implementation of the method described in the paper:
+    https://www2.ak.tu-berlin.de/~akgroup/ak_pub/2018/000458.pdf
+
+    Parameters
+    -----------
+    image_source_loc : array_like
+        Locations of image sources.
+    wall_flips: array_like
+        Number of x, y, z flips for each image source.
+    mic_loc: array_like
+        Microphone location.
+
+    Returns
+    -------
+    azimuth : :py:class:`~numpy.ndarray`
+        Azimith for each image source, in radians
+    colatitude : :py:class:`~numpy.ndarray`
+        Colatitude for each image source, in radians.
+
+    """
+
+    image_source_loc = np.array(image_source_loc)
+    wall_flips = np.array(wall_flips)
+    mic_loc = np.array(mic_loc)
+
+    dim, n_sources = image_source_loc.shape
+    assert wall_flips.shape[0] == dim
+    assert mic_loc.shape[0] == dim
+
+    p_vector_array = image_source_loc - np.array(mic_loc)[:, np.newaxis]
+    d_array = np.linalg.norm(p_vector_array, axis=0)
+
+    # Using (12) from the paper
+    power_array = np.ones_like(image_source_loc) * -1
+    power_array = np.power(power_array, (wall_flips + np.ones_like(image_source_loc)))
+    p_dash_array = p_vector_array * power_array
+
+    # Using (13) from the paper
+    azimuth = np.arctan2(p_dash_array[1], p_dash_array[0])
+    if dim == 2:
+        colatitude = np.ones(n_sources) * np.pi / 2
+    else:
+        colatitude = np.pi / 2 - np.arcsin(p_dash_array[2] / d_array)
+
+    return azimuth, colatitude
+
+
 def compute_ism_rir(
     src,
     mic,