diff --git a/src/magic.jl b/src/magic.jl
index 880aea2..8247e6c 100644
--- a/src/magic.jl
+++ b/src/magic.jl
@@ -11,6 +11,8 @@
 using GGDUtils:
     get_prop_with_grid_subset_index, get_subset_centers, get_grid_subset, interp
 
+export magic_nesep
+
 """
     magic_nesep(
         dd::IMASDD.dd;
@@ -28,28 +30,32 @@ dd: a data dictionary instance from IMASDD
 grid_ggd_idx: index of the grid_ggd to use. Many cases will have only one grid_ggd.
 space_number: space index to use. You may only have one space.
 midplane_grid_subset: index of the grid subset for y-aligned faces at the outboard
-    midplane. This is normally 11 and should not deviate from the standard without a reason.
+midplane. This is normally 11 and should not deviate from the standard without a reason.
 separatrix_grid_subset: index of the subset for x-aligned faces at the separatrix. Leave
-    this alone at 16 unless there's a reason for changing it.
+this alone at 16 unless there's a reason for changing it.
 
 Output:
 nesep: an array of n_e,sep values in m^-3 vs time. If the dd has only one timeslice,
-    the array will only have one element but it will still be an array.
+the array will only have one element but it will still be an array.
 """
-
 function magic_nesep(
     dd::IMASDD.dd;
     grid_ggd_idx::Int64=1,
     space_number::Int64=1,
+    cell_grid_subset::Int64=5,
     midplane_grid_subset::Int64=11,
     separatrix_grid_subset::Int64=16,
 )::Array{Float64}
-    nslices = length(dd.edge_profiles.grid_ggd)
+    # Form empty output array
+    nslices = length(dd.edge_profiles.ggd)
     nesep = Array{Float64}(undef, nslices)
 
+    # Get out references to basic ggd objects
     grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx]
     space = grid_ggd.space[space_number]
+    cell_sub = get_grid_subset(grid_ggd, cell_grid_subset)
 
+    # Find the intersection of the outboard midplane slice with the separatrix
     midplane_sub = get_grid_subset(grid_ggd, midplane_grid_subset)
     separatrix_sub = get_grid_subset(grid_ggd, separatrix_grid_subset)
     intersect_element = subset_do(
@@ -60,15 +66,14 @@ function magic_nesep(
     rintersect = intersect_point[1]
     zintersect = intersect_point[2]
 
-    nesep = Array{Float64}(undef, nslices)
     for it ∈ 1:nslices
         ggd = dd.edge_profiles.ggd[it]
-        nemidplane =
+        ne =
             get_prop_with_grid_subset_index(
                 ggd.electrons.density,
-                midplane_grid_subset,
+                cell_grid_subset,
             ).values
-        nesep[it] = interp(nemidplane, space, midplane_sub)(rintersect, zintersect)
+        nesep[it] = interp(ne, space, cell_sub)(rintersect, zintersect)
     end
 
     return nesep
diff --git a/test/runtests.jl b/test/runtests.jl
index 724d0db..e83d1b7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,6 @@
 using SynthDiag: IMASDD, add_interferometer!, add_langmuir_probes!, add_gas_injection!,
     compute_gas_injection!, get_gas_injection_response, Noise, OverwriteAttemptError,
-    langmuir_probe_current
+    langmuir_probe_current, magic_nesep
 using IMASDD: json2imas
 using Test
 using Printf
@@ -24,6 +24,9 @@ function parse_commandline()
         ["--gas_injection"],
         Dict(:help => "Test only gas injection",
             :action => :store_true),
+        ["--magic"],
+        Dict(:help => "Test only magic diagnostics",
+            :action => :store_true),
     )
     args = ArgParse.parse_args(localARGS, s)
     if !any(values(args)) # If no flags are set, run all tests
@@ -337,3 +340,14 @@ if args["gas_injection"]
         @test true
     end
 end
+
+if args["magic"]
+    @testset "magic_diagnostics" begin
+        ids = json2imas(
+            "$(@__DIR__)/../samples/time_dep_edge_profiles_with_equilibrium.json",
+        )
+        nt = length(ids.edge_profiles.ggd)
+        nesep = magic_nesep(ids; cell_grid_subset=-5)
+        @test length(nesep) == nt
+    end
+end