Nonintuitive indices on positions of head_related_impulse_responses
## General

* pyfar version: 0.5.1
* Python version: 3.10
* Operating System: Windows 11
* Did you install pyfar via pip: yes

## Description

When using `pf.signals.files.head_related_impulse_responses` to get HRIRs for multiple positions, I expected the indices to align with the order of items in the passed list. 

E.g. `[270,0],[0,0],[90,0]` => same order indices `0,1,2`.

Instead:
- index 0 => position [0,0]
- index 1 => position [90,0]
- index 2 => position [270,0]

## What I Did

```
import pyfar as pf
import matplotlib.pyplot as plt
hrir, _ = pf.signals.files.head_related_impulse_responses([[270,0],[0,0],[90,0]],diffuse_field_compensation=True)
pf.plot.time_freq(hrir[0])
pf.plot.time_freq(hrir[1])
pf.plot.time_freq(hrir[2])
plt.show()
# -> looking at the plots reveals the "non-matching" order

# alternatively:
import sounddevice
out=pf.dsp.convolve(hrir[0],pf.signals.files.speech())
sounddevice.play(out.time.T,out.sampling_rate,blocking=True)
# -> listening reveals the same
```

## Labels

Label your issue to make it easier for us to assign and track:
**hot**

Use a second label to let us know to which module the issue is related:
**signals**

## Project

**Do not assign your issue to a project.** This is done automatically.

diff --git a/pyfar/classes/coordinates.py b/pyfar/classes/coordinates.py
index a1c2e25a..75033e2f 100644
--- a/pyfar/classes/coordinates.py
+++ b/pyfar/classes/coordinates.py
@@ -960,7 +960,7 @@ class Coordinates():
 
         Returns
         -------
-        index : numpy array of ints
+        index : tuple of numpy arrays
             The indices of the selected points as a tuple of arrays. The length
             of the tuple matches :py:func:`~cdim`. The length of each array
             matches the number of selected points.
@@ -1028,7 +1028,7 @@ class Coordinates():
         if show:
             self.show(mask)
 
-        index = np.asarray(mask).nonzero()
+        index = np.where(mask)
 
         return index, mask
 
diff --git a/pyfar/signals/files.py b/pyfar/signals/files.py
index a57abb50..4a62d0c0 100644
--- a/pyfar/signals/files.py
+++ b/pyfar/signals/files.py
@@ -391,25 +391,28 @@ def head_related_impulse_responses(
 
     # get indices of source positions
     if position == "horizontal":
-        _, mask = sources.find_slice('elevation', 'deg', 0)
+        idx, _ = sources.find_slice('elevation', 'deg', 0)
     elif position == "median":
-        _, mask = sources.find_slice('lateral', 'deg', 0)
+        idx, _ = sources.find_slice('lateral', 'deg', 0)
+        # sort positions according to polar angle
+        polar = sources.get_sph("side", "deg")[idx, 1].flatten()
+        idx = (idx[0][np.argsort(polar)], )
     else:
-        mask = np.full((358, ), False)
+        idx = []
         for pos in position:
-            _, mask_current = sources.find_nearest_sph(
+            idx_current, _ = sources.find_nearest_sph(
                 pos[0], pos[1], 1.7, distance=0,
                 domain="sph", convention="top_elev", unit="deg")
-            if np.any(mask_current):
-                mask = mask | mask_current
+            if idx_current:
+                idx.append(idx_current[0])
             else:
                 raise ValueError((
                     f"HRIR for azimuth={pos[0]} and elevation={pos[1]} degrees"
                     " is not available. See help for more information."))
 
     # select data for desired source positions
-    hrirs.time = hrirs.time[mask]
-    sources = sources[mask]
+    hrirs.time = hrirs.time[idx]
+    sources = sources[idx]
 
     # diffuse field compensation
     if diffuse_field_compensation:
diff --git a/tests/test_coordinates.py b/tests/test_coordinates.py
index 8fca2ce5..7b59a75a 100644
--- a/tests/test_coordinates.py
+++ b/tests/test_coordinates.py
@@ -565,33 +565,39 @@ def test_find_slice():
 
     c = Coordinates(d, 0, 0)
     index, mask = c.find_slice('x', 'met', 0, 1)
-    npt.assert_allclose(index[0], np.array([1, 2, 3]))
+    npt.assert_allclose(index, (np.array([1, 2, 3]), ))
     npt.assert_allclose(mask, np.array([0, 1, 1, 1, 0]))
 
     c = Coordinates(0, d, 0)
     index, mask = c.find_slice('y', 'met', 0, 1)
-    npt.assert_allclose(index[0], np.array([1, 2, 3]))
+    npt.assert_allclose(index, (np.array([1, 2, 3]), ))
     npt.assert_allclose(mask, np.array([0, 1, 1, 1, 0]))
 
     c = Coordinates(0, 0, d)
     index, mask = c.find_slice('z', 'met', 0, 1)
-    npt.assert_allclose(index[0], np.array([1, 2, 3]))
+    npt.assert_allclose(index, (np.array([1, 2, 3]), ))
     npt.assert_allclose(mask, np.array([0, 1, 1, 1, 0]))
 
+    # cartesian grid, multi-dimensional coordinates
+    c = Coordinates([[0, 1], [1, 0]], 2, 3)
+    index, mask = c.find_slice('x', 'met', 0)
+    npt.assert_allclose(index, ([0, 1], [0, 1]))
+    npt.assert_allclose(mask, np.array([[1, 0], [0, 1]]))
+
     # spherical grid
     d = [358, 359, 0, 1, 2]
     c = Coordinates(d, 0, 1, 'sph', 'top_elev', 'deg')
     # cyclic query for lower bound
     index, mask = c.find_slice('azimuth', 'deg', 0, 1)
-    npt.assert_allclose(index[0], np.array([1, 2, 3]))
+    npt.assert_allclose(index, (np.array([1, 2, 3]), ))
     npt.assert_allclose(mask, np.array([0, 1, 1, 1, 0]))
     # cyclic query for upper bound
     index, mask = c.find_slice('azimuth', 'deg', 359, 2)
-    npt.assert_allclose(index[0], np.array([0, 1, 2, 3]))
+    npt.assert_allclose(index, (np.array([0, 1, 2, 3]), ))
     npt.assert_allclose(mask, np.array([1, 1, 1, 1, 0]))
     # non-cyclic query
     index, mask = c.find_slice('azimuth', 'deg', 1, 1)
-    npt.assert_allclose(index[0], np.array([2, 3, 4]))
+    npt.assert_allclose(index, (np.array([2, 3, 4]), ))
     npt.assert_allclose(mask, np.array([0, 0, 1, 1, 1]))
     # out of range query
     with raises(AssertionError):
@@ -608,6 +614,17 @@ def test_find_slice():
     plt.close("all")
 
 
+@pytest.mark.parametrize("coordinates,desired", [
+    (Coordinates([0, 1], 2, 3), [0, 2, 3]),
+    (Coordinates([[0, 1], [1, 0]], 2, 3), [[0, 2, 3], [0, 2, 3]])])
+def test_find_slice_slicing(coordinates, desired):
+    """Test if return values can be used for slicing"""
+
+    index, mask = coordinates.find_slice('x', 'met', 0)
+    assert coordinates[index] == coordinates[mask]
+    npt.assert_equal(coordinates[index].get_cart(), np.atleast_2d(desired))
+
+
 @pytest.mark.parametrize("rot_type,rot", [
     ('quat', [0, 0, 1 / np.sqrt(2), 1 / np.sqrt(2)]),
     ('matrix',  [[0, -1, 0], [1, 0, 0], [0, 0, 1]]),
diff --git a/tests/test_signals_files.py b/tests/test_signals_files.py
index 3700fa10..84c41a21 100644
--- a/tests/test_signals_files.py
+++ b/tests/test_signals_files.py
@@ -53,7 +53,8 @@ def test_speech(sampling_rate):
 
 @pytest.mark.parametrize('position,convention,first,second', [
     ([[0, 20]], 'top_elev', [0], [20]),
-    ([[0, 20], [10, 0]], 'top_elev', [0, 10], [0, 20]),
+    ([[0, 20], [10, 0]], 'top_elev', [0, 10], [20, 0]),
+    ([[90, 0], [0, 0]], 'top_elev', [90, 0], [0, 0]),
     ('horizontal', 'top_elev', np.arange(0, 180)*2, np.zeros(180)),
     ('median', 'side', np.zeros(180), np.arange(-45, 135)*2)])
 def test_hrirs_position(position, convention, first, second):
@@ -67,8 +68,8 @@ def test_hrirs_position(position, convention, first, second):
 
     # test source positions
     sg = sources.get_sph(convention, 'deg')
-    npt.assert_allclose(first, np.sort(sg[..., 0].flatten()), atol=1e-12)
-    npt.assert_allclose(second, np.sort(sg[..., 1].flatten()), atol=1e-12)
+    npt.assert_allclose(first, sg[..., 0].flatten(), atol=1e-12)
+    npt.assert_allclose(second, sg[..., 1].flatten(), atol=1e-12)
 
 
 def test_hrirs_diffuse_field_compensation():
['pyfar/classes/coordinates.py', 'tests/test_coordinates.py', 'tests/test_signals_files.py', 'pyfar/signals/files.py']