diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl
index 9bb0bd3ae..65a9307fe 100644
--- a/src/Equilibrium/Equilibrium.jl
+++ b/src/Equilibrium/Equilibrium.jl
@@ -625,8 +625,8 @@ function equilibrium_gse!(equil::PlasmaEquilibrium)
             gse_data = Array{Float32,3}(undef, mpsi + 1, mtheta + 1, 7)
             for ipsi in 0:mpsi
                 for itheta in 0:mtheta
-                    gse_data[ipsi+1, itheta+1, 1] = Float32(flux.ys[itheta+1])
-                    gse_data[ipsi+1, itheta+1, 2] = Float32(flux.xs[ipsi+1])
+                    gse_data[ipsi+1, itheta+1, 1] = Float32(equil.rzphi_ys[itheta+1])
+                    gse_data[ipsi+1, itheta+1, 2] = Float32(equil.rzphi_xs[ipsi+1])
                     gse_data[ipsi+1, itheta+1, 3] = Float32(flux_fs[ipsi+1, itheta+1, 1])
                     gse_data[ipsi+1, itheta+1, 4] = Float32(flux_fs[ipsi+1, itheta+1, 2])
                     gse_data[ipsi+1, itheta+1, 5] = Float32(source[ipsi+1, itheta+1])
@@ -639,7 +639,7 @@ function equilibrium_gse!(equil::PlasmaEquilibrium)
 
         # Write integrated error criterion
         h5open(joinpath(dirname(equil.config.eq_filename), "gsei.h5"), "w") do file
-            file["xs"] = Float32.(flux.xs)
+            file["xs"] = Float32.(equil.rzphi_xs)
             file["term"] = Float32.(term)
             file["totali"] = Float32.(totali)
             file["errori"] = Float32.(errori)
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index 2d33380c6..6468eaa18 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -1,3 +1,4 @@
+using HDF5
 
 @testset "Equilibrium Unit Tests" begin
 
@@ -594,4 +595,18 @@
         @test occursin("sol.retcode != ReturnCode.Success", src)
         @test occursin("isapprox(sol.t[end], 2π", src)
     end
+
+    @testset "equilibrium_gse! writes gse.h5 with diagnose_src" begin
+        dir = mktempdir()
+        cp(joinpath(data_dir, "EQDSK_COCOS_02"), joinpath(dir, "EQDSK_COCOS_02"))
+        config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(;
+            eq_filename=joinpath(dir, "EQDSK_COCOS_02"), eq_type="efit", jac_type="boozer",
+            grid_type="ldp", psilow=0.01, psihigh=0.994)
+        eq = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(config)
+        eq.params.diagnose_src = true
+        GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_gse!(eq)
+        @test all(isfile(joinpath(dir, f)) for f in ("gsec.h5", "gse.h5", "gsei.h5"))
+        # Column 2 is ψ, as plot_gse_by_theta reads it
+        @test h5read(joinpath(dir, "gse.h5"), "gse_data")[:, 1, 2] == Float32.(eq.rzphi_xs)
+    end
 end