d-burg · Oct 6, 2026
diff --git a/docs/src/conventions.md b/docs/src/conventions.md
@@ -106,11 +106,9 @@ helicity = bt_sign * Int(sign(crnt))   # src/ForcingTerms/CoilFourier.jl
 ## ``F = R B_\phi``
 
 The poloidal-current function ``F = R B_\phi`` from the Grad-Shafranov equation is **forced
-positive** via `abs` (`read_eq_efit` in `src/Equilibrium/ReadEquilibrium.jl`):
-
-```julia
-abs.(fpol_data)
-```
+positive**: the g-file and IMAS readers take its magnitude (`file_profiles` in
+`src/Equilibrium/ReadEquilibrium.jl`), whether it is integrated from the file's ``FF'`` or read
+from the tabulated values.
 
 The code always works with ``|F|``; the sign of ``B_t`` is carried separately (`fpol_sign`, and the
 `bt_sign`/`crnt` helicity inputs).

diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md
@@ -55,6 +55,31 @@ The Cartesian evaluation grid is clipped to the separatrix bounding box and its
 resolution is set adaptively from a bilinear interpolation error bound, so no manual
 tuning is needed.
 
+## F and P from a g-file or IMAS equilibrium
+
+A g-file tabulates both the profiles (`FPOL`, `PRES`) and their flux derivatives (`FFPRIM`,
+`PPRIME`); an IMAS equilibrium does the same (`f`, `pressure`, `f_df_dpsi`, `dpressure_dpsi`).
+The stability matrices use ``F'`` directly and respond to how ``F'`` varies across each rational
+surface, i.e. to the current gradient. ``F`` itself changes by only a few percent across the
+plasma, so differentiating a tabulated ``F`` amplifies its rounding error by roughly
+``1/(h^2\delta)`` in the slope of ``F'`` (``h`` the node spacing, ``\delta`` the fractional
+variation of ``F``). On a 257-node file written in single precision that is an error of order
+10 % in the current gradient, and several units in ``\Delta'``.
+
+`profile_source` in `[Equilibrium]` selects where the two profiles come from:
+
+| Value | Behaviour |
+|---|---|
+| `"derivatives"` (default) | ``F`` and ``P`` are the integrals of the file's ``FF'`` and ``p'``, anchored at the boundary values of the tabulated ``F`` and ``P`` |
+| `"values"` | The tabulated ``F`` and ``P`` are used as written |
+
+With `"derivatives"` both profiles are always rebuilt together, so the pair stays consistent in
+the force balance, and everything downstream still differentiates one spline. The reader logs the
+largest interior difference between the rebuilt and the tabulated profiles and warns when it
+exceeds `1e-4` of ``F`` or 5 % of the peak pressure. It falls back to `"values"` with a warning
+when the derivative arrays are absent, all zero, or opposite in sign for only one of the two
+profiles. The two end nodes are used as written.
+
 ## Radial grid packing
 
 With `grid_type = "auto"` and `mpsi = 0` (the defaults; `"log_asymptotic"` is a legacy alias), the radial grid is

diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl
@@ -41,6 +41,11 @@ specified in the input.
   - `etol::Float64` - Error tolerance for equilibrium solver
   - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations)
   - `use_galgrid::Bool` - Use the same grid as galerkin method
+  - `imas_cocos::Int` - COCOS convention of the input IMAS `dd.equilibrium` (11 = IMAS standard, 2 = internal)
+  - `profile_source::String` - Which 1D arrays of a g-file or IMAS equilibrium define F and P.
+    `"derivatives"` (default) integrates the file's FF′ and p′ inward from the boundary values of
+    F and P; `"values"` uses the tabulated F and P directly. The reader falls back to `"values"`
+    with a warning when the derivatives are absent or unusable.
 """
 @kwdef struct EquilibriumConfig
     eq_type::String = "efit"
@@ -72,8 +77,8 @@ specified in the input.
     force_termination::Bool = false
     use_galgrid::Bool = true
 
-    # IMAS-specific: expected COCOS convention of the input dd.equilibrium (11=IMAS standard, 2=GPEC internal)
     imas_cocos::Int = 11
+    profile_source::String = "derivatives"
 
     """
     Modified internal constructor that enforces self consistency within the inputs
@@ -84,7 +89,7 @@ specified in the input.
     function EquilibriumConfig(eq_type, eq_filename, r0exp, b0exp, jac_type, _, _, _, _,
         jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc,
         grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, etol,
-        force_termination, use_galgrid, imas_cocos)
+        force_termination, use_galgrid, imas_cocos, profile_source)
         if jac_type == "hamada"
             @info "Forcing hamada coordinate jacobian exponents: power_*"
             power_b = 0
@@ -147,10 +152,12 @@ specified in the input.
             @warn "psihigh = $psihigh exceeds 1.0 (separatrix); clamping to 1.0"
         end
         psihigh = min(psihigh, 1.0)
+        profile_source in ("derivatives", "values") ||
+            error("Cannot recognize profile_source = $(profile_source); use \"derivatives\" or \"values\"")
         return new(eq_type, eq_filename, r0exp, b0exp, jac_type, power_bp, power_b, power_r, power_rc,
             jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc,
             grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, etol,
-            force_termination, use_galgrid, imas_cocos)
+            force_termination, use_galgrid, imas_cocos, profile_source)
     end
 end
 

diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl
@@ -40,6 +40,90 @@ function _read_1d_gfile_format(lines_block::Vector{String}, num_values::Int)
     return parsed_values
 end
 
+# Largest interior disagreement between the tabulated and the integrated profiles the reader accepts
+# without a warning: relative for F, as a fraction of the peak pressure for P.
+const PROFILE_F_MISMATCH_WARN = 1e-4
+const PROFILE_P_MISMATCH_WARN = 5e-2
+
+"""
+    integrate_profile_derivatives(xs, f, p, ffprime, pprime) -> Union{Nothing,NamedTuple}
+
+Rebuild F and P on the nodes `xs` of a normalized flux coordinate (0 at the axis, 1 at the
+boundary) by integrating the tabulated `ffprime` = d(F²/2)/dx and `pprime` = dP/dx inward from the
+boundary values of `f` and `p`:
+
+    F²(x) = F(1)² − 2∫ₓ¹ FF′ dx,    P(x) = P(1) − ∫ₓ¹ P′ dx.
+
+F varies by only a few percent across the plasma, so rounding in a tabulated F is amplified by
+roughly 1/(h²δ) in the slope of F′ (h the node spacing, δ the fractional variation of F), which is
+the current gradient the stability matrices respond to. The tabulated derivatives carry that
+information directly, and they are the source terms the Grad-Shafranov solution ψ(R,Z) was computed
+from. F and P are rebuilt together so that the pair stays consistent in the force balance.
+
+Returns `(; f, p, f_mismatch, p_mismatch)` with `f` positive, where the mismatches are the largest
+interior differences from the tabulated values (relative for F, as a fraction of the peak pressure
+for P; the two end nodes are excluded). Returns `nothing` when the derivatives cannot be used: a
+derivative array that is all zero or non-finite while its profile varies, a non-increasing `xs`, a
+non-positive F², or derivative signs that match the tabulated profiles for only one of the two.
+"""
+function integrate_profile_derivatives(xs::AbstractVector{<:Real}, f::AbstractVector{<:Real}, p::AbstractVector{<:Real},
+    ffprime::AbstractVector{<:Real}, pprime::AbstractVector{<:Real})
+    n = length(xs)
+    (n >= 4 && length(f) == n && length(p) == n && length(ffprime) == n && length(pprime) == n) || return nothing
+    (all(isfinite, ffprime) && all(isfinite, pprime) && all(>(0), diff(xs))) || return nothing
+    f_abs = abs.(f)
+    f_has, p_has = any(!iszero, ffprime), any(!iszero, pprime)
+    # A zero derivative array is only believable for a profile that is itself flat
+    (f_has || maximum(f_abs) - minimum(f_abs) <= eps(maximum(f_abs))) || return nothing
+    (p_has || maximum(p) == minimum(p)) || return nothing
+    (f_has || p_has) || return nothing
+
+    nodes = collect(Float64, xs)
+    tail(y) = (c = FastInterpolations.cumulative_integrate(cubic_interp(nodes, collect(Float64, y))); c[end] .- c)
+    f_tail, p_tail = tail(ffprime), tail(pprime)
+    interior = 2:(n-1)
+    p_scale = maximum(abs, p)
+    f_error(sgn) = maximum(i -> abs(sqrt(max(f_abs[end]^2 - 2 * sgn * f_tail[i], 0.0)) - f_abs[i]) / f_abs[i], interior)
+    p_error(sgn) = p_scale == 0 ? 0.0 : maximum(i -> abs(p[end] - sgn * p_tail[i] - p[i]), interior) / p_scale
+
+    # Writers differ in the sign convention of the flux derivative, so take the sign that reproduces
+    # the tabulated profiles and require F and P to agree on it.
+    f_sign = f_error(1) <= f_error(-1) ? 1 : -1
+    p_sign = p_error(1) <= p_error(-1) ? 1 : -1
+    f_has && p_has && f_sign != p_sign && return nothing
+    sgn = f_has ? f_sign : p_sign
+
+    f2 = f_abs[end]^2 .- 2 .* sgn .* f_tail
+    all(>(0), f2) || return nothing
+    return (; f=sqrt.(f2), p=p[end] .- sgn .* p_tail, f_mismatch=f_error(sgn), p_mismatch=p_error(sgn))
+end
+
+"""
+    file_profiles(config, xs, f, p, ffprime, pprime) -> (f_nodes, p_nodes)
+
+F (as a magnitude) and P on the 1D nodes of a file-based direct equilibrium, from the arrays that
+`config.profile_source` selects. With `"derivatives"` the profiles come from
+`integrate_profile_derivatives`, and the tabulated values are used with a warning when
+that is not possible.
+"""
+function file_profiles(config::EquilibriumConfig, xs::AbstractVector{<:Real}, f::AbstractVector{<:Real}, p::AbstractVector{<:Real},
+    ffprime::AbstractVector{<:Real}, pprime::AbstractVector{<:Real})
+    config.profile_source == "values" && return abs.(f), collect(Float64, p)
+    built = integrate_profile_derivatives(xs, f, p, ffprime, pprime)
+    if built === nothing
+        @warn "profile_source = \"derivatives\", but the file's FF′ and p′ are absent or unusable; using its tabulated F and P"
+        return abs.(f), collect(Float64, p)
+    end
+    msg = "F and P integrated from the file's FF′ and p′; largest interior difference from the tabulated values: " *
+          "$(@sprintf("%.1e", built.f_mismatch)) of F, $(@sprintf("%.1e", built.p_mismatch)) of the peak pressure"
+    if built.f_mismatch > PROFILE_F_MISMATCH_WARN || built.p_mismatch > PROFILE_P_MISMATCH_WARN
+        @warn msg * ". The file's profiles and derivatives disagree; profile_source = \"values\" selects the tabulated profiles."
+    else
+        @info msg
+    end
+    return built.f, built.p
+end
+
 """
     _read_efit(equil_in)
 
@@ -91,18 +175,20 @@ function read_efit(config::EquilibriumConfig)
     psi_rz = reshape(psi_flat_vec, nw, nh)
 
     # --- Create 1D Profile Spline (sq_in) ---
+    psio_signed = sibry - simag
     psi_norm_grid = range(0.0, 1.0; length=nw)
+    # FFPRIM and PPRIME are derivatives in the file's ψ; psio_signed converts them to ψ_norm
+    f_nodes, p_nodes = file_profiles(config, psi_norm_grid, fpol_data, pres_data, ffprime_data .* psio_signed, pprime_data .* psio_signed)
     sq_fs_nodes = hcat(
-        abs.(fpol_data),
-        max.(pres_data .* mu0, 0.0),
+        f_nodes,
+        max.(p_nodes .* mu0, 0.0),
         qprof_data,
         sqrt.(psi_norm_grid)
     )
     sq_xs = collect(psi_norm_grid)
     sq_in = cubic_interp(sq_xs, Series(sq_fs_nodes); extrap=ExtendExtrap())
 
     # --- Process and Normalize 2D Psi Data ---
-    psio_signed = sibry - simag
     psi_proc = (sibry .- psi_rz)
     psio = abs(psio_signed)
     # Ensure psi at the magnetic axis is positive relative to the boundary
@@ -495,11 +581,17 @@ function read_imas(config::EquilibriumConfig, dd)
     psi_norm_grid = range(0.0, 1.0; length=nw)
 
     # Build equilibrium source terms for spline interpolation
-    # abs(f_1d): F(ψ) can be negative depending on toroidal field direction convention;
-    #            take absolute value since GPEC uses magnitude F = R·|Bt|
+    # F(ψ) can be negative depending on toroidal field direction convention; GPEC uses the
+    # magnitude F = R·|Bt|. f_df_dpsi and dpressure_dpsi are derivatives in the stored ψ, so the
+    # stored ψ span converts them to its normalized coordinate (absent arrays read as zero).
+    psi_stored = eqt.profiles_1d.psi
+    psi_span = psi_stored[end] - psi_stored[1]
+    stored_or_zero(name) = (v = getproperty(eqt.profiles_1d, name, Float64[]); length(v) == nw ? v .* psi_span : zeros(nw))
+    f_nodes, p_nodes = file_profiles(config, (psi_stored .- psi_stored[1]) ./ psi_span, f_1d, p_1d,
+        stored_or_zero(:f_df_dpsi), stored_or_zero(:dpressure_dpsi))
     sq_fs_nodes = hcat(
-        abs.(f_1d),
-        max.(p_1d .* mu0, 0.0),
+        f_nodes,
+        max.(p_nodes .* mu0, 0.0),
         q_1d,
         sqrt.(psi_norm_grid)
     )

diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
@@ -147,6 +147,56 @@
             (Symbol(k) => v for (k, v) in ffs_table)...) isa GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl
     end
 
+    @testset "F and P from file derivatives" begin
+        Eq = GeneralizedPerturbedEquilibrium.Equilibrium
+        # Smooth profiles with exact derivatives; F varies by ~2 % like a tokamak F
+        xs = range(0.0, 1.0; length=257)
+        f_exact = @. sqrt(9.0 + 0.4 * (1 - xs)^2 - 0.1 * (1 - xs)^3)
+        ffprime = @. -0.4 * (1 - xs) + 0.15 * (1 - xs)^2
+        f2_exact = @. (ffprime^2 / f_exact^2 * -1 + 0.4 - 0.3 * (1 - xs)) / f_exact   # d²F/dx²
+        p_exact = @. 5e4 * (1 - xs)^2 + 300.0
+        pprime = @. -1e5 * (1 - xs)
+        f_single = Float64.(Float32.(f_exact))   # a table written in single precision
+
+        built = Eq.integrate_profile_derivatives(xs, f_single, p_exact, ffprime, pprime)
+        @test maximum(abs.(built.f .- f_exact) ./ f_exact) < 1e-7   # boundary anchor carries the rounding
+        @test maximum(abs.(built.p .- p_exact)) < 1e-6 * 5e4
+        @test built.f_mismatch < 1e-6 && built.p_mismatch < 1e-9
+
+        # The slope of F′ is what single-precision rounding ruins and the derivative route keeps
+        slope_error(f) = maximum(abs.(deriv2(cubic_interp(collect(xs), f)).(xs[20:end-20]) .- f2_exact[20:end-20]))
+        @test slope_error(built.f) < 1e-2 * slope_error(f_single)
+
+        # Either sign convention of the flux derivative gives the same profiles
+        flipped = Eq.integrate_profile_derivatives(xs, f_single, p_exact, -ffprime, -pprime)
+        @test flipped.f ≈ built.f && flipped.p ≈ built.p
+        @test Eq.integrate_profile_derivatives(xs, -f_single, p_exact, ffprime, pprime).f ≈ built.f
+
+        # Unusable derivatives: absent, non-finite, or signs that suit only one of the two profiles
+        @test Eq.integrate_profile_derivatives(xs, f_single, p_exact, zero(ffprime), zero(pprime)) === nothing
+        @test Eq.integrate_profile_derivatives(xs, f_single, p_exact, ffprime, zero(pprime)) === nothing
+        @test Eq.integrate_profile_derivatives(xs, f_single, p_exact, fill(NaN, 257), pprime) === nothing
+        @test Eq.integrate_profile_derivatives(xs, f_single, p_exact, ffprime, -pprime) === nothing
+
+        derivs = Eq.EquilibriumConfig(; profile_source="derivatives")
+        values = Eq.EquilibriumConfig(; profile_source="values")
+        @test_throws ErrorException Eq.EquilibriumConfig(; profile_source="spline")
+        @test Eq.file_profiles(values, xs, -f_single, p_exact, ffprime, pprime) == (f_single, p_exact)
+        @test (@test_logs (:warn, r"absent or unusable") Eq.file_profiles(derivs, xs, f_single, p_exact, zero(ffprime), zero(pprime))) == (f_single, p_exact)
+        @test (@test_logs (:info, r"integrated from") Eq.file_profiles(derivs, xs, f_single, p_exact, ffprime, pprime))[1] ≈ built.f
+        # A pressure table 20 % off its own derivative is reported, not silently accepted
+        @test_logs (:warn, r"disagree") Eq.file_profiles(derivs, xs, f_single, p_exact, ffprime, 1.2 .* pprime)
+
+        # The g-file reader keeps the boundary values and moves the interior only at the rounding level
+        gfile = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example", "TkMkr_D3Dlike_Hmode.geqdsk")
+        table(source) = Eq.read_efit(Eq.EquilibriumConfig(; eq_type="efit", eq_filename=gfile, profile_source=source)).ingest.sq_fs
+        tabulated, integrated = table("values"), table("derivatives")
+        @test integrated[end, 1:2] == tabulated[end, 1:2]
+        @test integrated[:, 1] != tabulated[:, 1]
+        @test maximum(abs.(integrated[:, 1] .- tabulated[:, 1]) ./ tabulated[:, 1]) < 1e-5
+        @test integrated[:, 3:4] == tabulated[:, 3:4]
+    end
+
     @testset "EFIT Method Consistency" begin
         # All three methods solve the same equilibrium — q-profiles should broadly agree.
         # Tolerance is 10% to allow for method-specific discretisation differences.

diff --git a/test/runtests_imas.jl b/test/runtests_imas.jl
@@ -72,6 +72,35 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
         @test isapprox(result.psio, psi_bnd_int; rtol=1e-6)
     end
 
+    # read_imas — F and P come from the stored derivatives when the dd carries them
+    @testset "read_imas: profiles from stored derivatives" begin
+        for cocos in (11, 2)
+            dd, _ = make_mock_dd(; cocos=cocos)
+            prof = dd.equilibrium.time_slice[1].profiles_1d
+            nw = length(prof.psi)
+            # Derivatives in the stored ψ that match the mock's flat F and linear pressure
+            prof.f_df_dpsi = zeros(nw)
+            prof.dpressure_dpsi = fill(-1e4 / (prof.psi[end] - prof.psi[1]), nw)
+            # A pressure table carrying an error the derivatives do not have
+            exact_pressure = copy(prof.pressure)
+            prof.pressure = exact_pressure .* (1 .+ 1e-3 .* sinpi.(range(0, 8; length=nw)))
+
+            mu0 = 4π * 1e-7
+            read_mu0p(source) = Equilibrium.read_imas(
+                EquilibriumConfig(; eq_type="imas", eq_filename="mock", imas_cocos=cocos, profile_source=source), dd).ingest.sq_fs[:, 2]
+            @test read_mu0p("derivatives") ≈ mu0 .* exact_pressure rtol = 1e-9 atol = 1e-12
+            @test read_mu0p("values") == mu0 .* prof.pressure
+        end
+    end
+
+    # read_imas — a dd without the derivative arrays falls back to the tabulated profiles
+    @testset "read_imas: fallback without stored derivatives" begin
+        dd, _ = make_mock_dd()
+        config = EquilibriumConfig(; eq_type="imas", eq_filename="mock")
+        result = @test_logs (:warn, r"absent or unusable") match_mode = :any Equilibrium.read_imas(config, dd)
+        @test result.ingest.sq_fs[:, 1] == fill(5.0, 64)
+    end
+
     # Test 3: read_imas — splines are callable after construction
     @testset "read_imas: splines callable" begin
         dd, _ = make_mock_dd()