diff --git a/pyfstat/optimal_setup_functions.py b/pyfstat/optimal_setup_functions.py
index dda36c233bebb4cb9d95daee0819c34ca83d19d2..dc0c5fb65cf76b8d756a31885bdbab28aa7db34f 100644
--- a/pyfstat/optimal_setup_functions.py
+++ b/pyfstat/optimal_setup_functions.py
@@ -226,13 +226,21 @@ def get_Nstar_estimate(
     lalpulsar.ConvertPhysicalToSuperskyPoints(
         out_rssky, in_phys, SSkyMetric.semi_rssky_transf)
 
-    parallelepiped = (out_rssky.data[i:, 1:].T - out_rssky.data[i:, 0]).T
+    d = out_rssky.data
 
-    sqrtdetG = np.sqrt(np.linalg.det(
-        SSkyMetric.semi_rssky_metric.data[i:, i:]))
+    g = SSkyMetric.semi_rssky_metric.data
 
-    dV = np.abs(np.linalg.det(parallelepiped))
+    d[2:] = d[2:][::-1]  # Convert to Alpha, Delta, F0, F1.. ordering
+    g[2:] = g[2:][::-1]  # Convert to Alpha, Delta, F0, F1.. ordering
 
-    Nstar = sqrtdetG * dV
+    g = g[i:, i:]  # Remove sky if required
+    parallelepiped = (out_rssky.data[i:, 1:].T - out_rssky.data[i:, 0]).T
 
-    return Nstar
+    Nstars = []
+    for j in range(1, len(g)+1):
+        dV = np.abs(np.linalg.det(parallelepiped[:j, :j]))
+        sqrtdetG = np.sqrt(np.abs(np.linalg.det(g[:j, :j])))
+        Nstars.append(sqrtdetG * dV)
+    logging.debug('Nstar for each dimension = {}'.format(
+        ', '.join(["{:1.1e}".format(n) for n in Nstars])))
+    return np.max(Nstars)