diff --git a/imagesets/multibandmsimageset.cpp b/imagesets/multibandmsimageset.cpp
index 3fb64ecdec7fd600f1d77b69e189246f14797237..c2ea070db84794f17bcebb8851bb9d1f23f920a9 100644
--- a/imagesets/multibandmsimageset.cpp
+++ b/imagesets/multibandmsimageset.cpp
@@ -197,51 +197,62 @@ static std::vector<const MSMetaData*> GetInitializedMetaData(
 }
 
 namespace {
-// These two helper functions are copied from
+// These two helper functions below are based on
 // https://stackoverflow.com/questions/17074324/how-can-i-sort-two-vectors-in-the-same-way-with-criteria-that-uses-only-one-of
+
+// They are used in MultiBandMsImageSet::ProcessMetaData() to sort multiple
+// vectors simulatenously by frequency. There is a nicer ranges based solution,
+// supported from C++23 onwards, see the comment where the functions are used.
+// The functions here can be removed when moving to to ranges based solution
+
+/**
+ * Returns an indexing vector that sorts vector \p input according to comparator
+ * \p compare
+ */
 template <typename T, typename Compare>
-std::vector<std::size_t> sort_permutation(const std::vector<T>& vec,
-                                          const Compare& compare) {
-  std::vector<std::size_t> p(vec.size());
-  std::iota(p.begin(), p.end(), 0);
-  std::sort(p.begin(), p.end(), [&](std::size_t i, std::size_t j) {
-    return compare(vec[i], vec[j]);
-  });
-  return p;
+std::vector<std::size_t> MakeSortingPermutation(const std::vector<T>& input,
+                                                const Compare& compare) {
+  std::vector<std::size_t> permutation(input.size());
+  std::iota(permutation.begin(), permutation.end(), 0);
+  std::sort(permutation.begin(), permutation.end(),
+            [&](std::size_t i, std::size_t j) {
+              return compare(input[i], input[j]);
+            });
+  return permutation;
 }
 
+/**
+ * Permute the vector \p input using indexing vector \p permutation.
+ */
 template <typename T>
-void apply_permutation_in_place(std::vector<T>& vec,
-                                const std::vector<std::size_t>& p) {
-  std::vector<bool> done(vec.size());
-  for (std::size_t i = 0; i < vec.size(); ++i) {
+void ApplyPermutation(std::vector<T>& input,
+                      const std::vector<std::size_t>& permutation) {
+  std::vector<bool> done(input.size());
+  for (std::size_t i = 0; i < input.size(); ++i) {
     if (done[i]) {
       continue;
     }
     done[i] = true;
     std::size_t prev_j = i;
-    std::size_t j = p[i];
+    std::size_t j = permutation[i];
     while (i != j) {
-      std::swap(vec[prev_j], vec[j]);
+      std::swap(input[prev_j], input[j]);
       done[j] = true;
       prev_j = j;
-      j = p[j];
+      j = permutation[j];
     }
   }
 }
 
-}  // anonymous namespace
-
-void MultiBandMsImageSet::ProcessMetaData() {
-  std::vector<const MSMetaData*> meta_data =
-      GetInitializedMetaData(readers_.begin(), readers_.end());
-
-  // The ms's should be combined in frequency order
-  // Either ascending or descending depending on the order of the channels in
-  // the ms
-
-  // Check whether frequency is ascending/descending
-  // Initially no order has been set
+/**
+ * Find the order in which the channels in the spws are sorted.
+ *
+ * Returns std::optional<bool>(true) when all spws are in ascending frequency
+ * order, std::optiona<bool>(false) when in descending frequency order.
+ * std::optional<bool>() when in mixed order
+ */
+std::optional<bool> IsInAscendingFrequencyOrder(
+    const std::vector<const MSMetaData*>& meta_data) {
   std::optional<bool> frequency_is_ascending;
   for (const MSMetaData* ms_meta_data : meta_data) {
     // Only check the first band, assuming the rest is in the same order
@@ -257,9 +268,7 @@ void MultiBandMsImageSet::ProcessMetaData() {
       // If an ordering has already been established, the current ms should be
       // conformal
       if (*frequency_is_ascending != ms_frequency_is_ascending) {
-        throw std::runtime_error(
-            "Trying to concatenate MeasurementSets that have not the same "
-            "ordering in frequency.");
+        return std::optional<bool>();
       }
     } else {
       // There was no ordering established, the current ms determines the
@@ -271,10 +280,42 @@ void MultiBandMsImageSet::ProcessMetaData() {
   // Set the ordering arbitrarily to ascending
   if (!frequency_is_ascending.has_value()) frequency_is_ascending = true;
 
+  return frequency_is_ascending;
+}
+
+}  // anonymous namespace
+
+void MultiBandMsImageSet::ProcessMetaData() {
+  std::vector<const MSMetaData*> meta_data =
+      GetInitializedMetaData(readers_.begin(), readers_.end());
+
+  // The ms's should be combined in frequency order
+  // Either ascending or descending depending on the order of the channels in
+  // the ms
+  std::optional<bool> frequency_is_ascending =
+      IsInAscendingFrequencyOrder(meta_data);
+  if (!frequency_is_ascending.has_value()) {
+    throw std::runtime_error(
+        "Trying to concatenate MeasurementSets that have not the same "
+        "ordering in frequency.");
+  }
+
+  // In C++23 multiple vectors can be sorted simultaneously using ranges:
+  //
+  //   std::ranges::sort(std::views::zip(ms_names_, readers_, meta_data),
+  //   [frequency_is_ascending](auto&& a, auto&& b) {
+  //     return (std::get<2>(a)->GetBandInfo(0).channels[0].frequencyHz <
+  //             std::get<2>(b)->GetBandInfo(0).channels[0].frequencyHz) !=
+  //             !(*frequency_is_ascending);
+  //   });
+  //
+  // Without ranges support, the helper functions MakeSortingPermutation and
+  // ApplyPermutation are needed
+
   // Get the permutation to sort the ms's in ascending/descending frequency
   // order The != operator is equivalent to a xor operation, toggling the
   // comparison
-  auto permutation = sort_permutation(
+  std::vector<size_t> permutation = MakeSortingPermutation(
       meta_data,
       [frequency_is_ascending](const MSMetaData* a, const MSMetaData* b) {
         return (a->GetBandInfo(0).channels[0].frequencyHz <
@@ -282,10 +323,10 @@ void MultiBandMsImageSet::ProcessMetaData() {
                !(*frequency_is_ascending);
       });
 
-  // Apply the permutation to ms_names_ and it's derivatives
-  apply_permutation_in_place(ms_names_, permutation);
-  apply_permutation_in_place(readers_, permutation);
-  apply_permutation_in_place(meta_data, permutation);
+  // Apply the permutation to ms_names_ and its derivatives
+  ApplyPermutation(ms_names_, permutation);
+  ApplyPermutation(readers_, permutation);
+  ApplyPermutation(meta_data, permutation);
 
   // These fields are only validated.
   ExtractField(meta_data, GetBaselines, "baselines");
diff --git a/test/integration/tConcatenateFrequency.py b/test/integration/tConcatenateFrequency.py
index 5e5edfd20787f11dcdf7f4e052e02c7fcc7b0c08..7bcaf180bd49b54baa7e49a4cb9eaf1fabbad49c 100644
--- a/test/integration/tConcatenateFrequency.py
+++ b/test/integration/tConcatenateFrequency.py
@@ -40,6 +40,8 @@ def test(read_mode):
         assert_taql(f"select from {MS}/HISTORY")
 
     # Run Aoflagger and validate output.
+    # The MSs are passed in reversed order to check that aoflagger sorts
+    # them along frequency before concatenating
     if read_mode != "":
         check_call(
             [config.aoflagger, "-concatenate-frequency", read_mode] + MSs[::-1]
@@ -80,7 +82,7 @@ def test_two_chunks(chunk_size):
 
     check_call(
         [config.aoflagger, "-concatenate-frequency", "-chunk-size", chunk_size]
-        + MSs[::-1]
+        + MSs
     )
 
     for MS in MSs:
@@ -117,7 +119,7 @@ def test_three_chunks_limited_interval(chunk_size):
             "-chunk-size",
             chunk_size,
         ]
-        + MSs[::-1]
+        + MSs
     )
 
     for MS in MSs: