              , d_pops
              , sizeof(float) * n_rows * n_radii
              , cudaMemcpyDeviceToHost);
-
-    //TODO sort tmp_pops into pops
-
+    // sort tmp_pops into pops
+    Pops pops;
+    for (unsigned int j=0; j < n_radii; ++j) {
+      pops[radii[j]].resize(n_rows);
+      for (unsigned int i=0; i < n_rows; ++i) {
+        pops[radii[j]][i] = tmp_pops[j*n_rows+i];
+    }
     cudaFree(d_sorted_coords);
     cudaFree(d_radii2);
     cudaFree(d_pops);

                                         , d_radii2
                                         , n_radii
                                         , d_in_radius);
      // compute pops / radius
      for (unsigned int r=0; r < n_radii; ++r) {
        //TODO define range / offset
        reduce_sum<BSIZE> <<< rng
                            , BSIZE
                            , 0
                            , streams[i_stream] >>> (offset

      cudaFree(d_in_radius[s]);
      cudaFree(d_coords_ref[s]);
    }
    return pops;
  }

  Pops
