diff --git a/src/Fields.jl b/src/Fields.jl
index 33663795..d7437b9b 100644
--- a/src/Fields.jl
+++ b/src/Fields.jl
@@ -493,6 +493,7 @@ julia> fields[1].fields[1].energy/energy ≈ 0.98071312
 true
 julia> fields[2].fields[1].energy/energy ≈ 0.0061826217
 true
+```
 """
 function gauss_beam_init(modes, k, ω0, fieldfunc; energy, kwargs...)
     gauss = normalised_gauss_beam(k, ω0)
@@ -797,7 +798,7 @@ function optcomp_taylor(Eω::AbstractVecOrMat, grid, λ0; order=2, boundfac=8)
     τ0FTL = Maths.fwhm(grid.t, ItFTL)/(2*sqrt(log(2)))
     τ0 = Maths.fwhm(grid.t, _It(iFT(Eω, grid), grid))/(2*sqrt(log(2)))
 
-    ϕ2_0 = τ0FTL*sqrt(τ0^2 - τ0FTL^2) # GDD to stretch Gaussian from FTL to actual duration
+    ϕ2_0 = τ0FTL*sqrt(abs(τ0^2 - τ0FTL^2)) # GDD to stretch Gaussian from FTL to actual duration
 
     # for Gaussian with pure GDD, sqrt(ϕ2_0) is the FTL duration, so use that as guide
     bounds = boundfac*(sqrt(ϕ2_0) .^(2:order))
diff --git a/src/Output.jl b/src/Output.jl
index af2bf0fd..95fb6643 100644
--- a/src/Output.jl
+++ b/src/Output.jl
@@ -252,32 +252,40 @@ end
 
 # more indices -> read slice of data
 function getindex(o::HDF5Output, ds::AbstractString,
-                  I::Union{AbstractRange, Int, Colon, Ellipsis}...)
+                  I::Union{AbstractRange, Int, Colon, Ellipsis, AbstractString, Array}...)
     HDF5.h5open(o.fpath, "r") do file
-        file[ds][to_indices(file[ds], I)...]
+        _getindex(file[ds], I...)
     end
 end
 
+function _getindex(parent::HDF5.Group, dss::AbstractString, I::Union{AbstractRange, Int, Colon, Ellipsis, Array}...)
+    _getindex(parent[dss], I...)
+end
+
+function _getindex(parent::HDF5.Group, dss::AbstractString)
+    read(parent[dss])
+end
+
+function _getindex(parent::HDF5.Dataset, I::Union{AbstractRange, Int, Colon, Ellipsis}...)
+    parent[to_indices(parent, I)...]
+end
+
 # indexing with an array, e.g. o["Eω", :, [1, 2, 3]] has to be handled separately
-function getindex(o::HDF5Output, ds::AbstractString,
-                  I::Union{AbstractRange, Int, Colon, Array, Ellipsis}...)
+function _getindex(parent::HDF5.Dataset, I::Union{AbstractRange, Int, Colon, Array, Ellipsis}...)
     if count(isa.(I, Array)) > 1
         error("Only one dimension can be index with an array.")
     end
-    HDF5.h5open(o.fpath, "r") do file
-        dset = file[ds]
-        idcs = to_indices(dset, I)
-        adim = findfirst(isa.(idcs, Array)) # which of the indices is the array
-        arr = idcs[adim] # the array itself
-        T = eltype(dset)
-        ret = Array{T}(undef, map(length, idcs))
-        Ilo = idcs[1:adim-1]
-        Ihi = idcs[adim+1:end]
-        for ii in eachindex(arr)
-            ret[Ilo..., ii, Ihi...] .= dset[Ilo..., arr[ii], Ihi...]
-        end
-        ret
-    end
+    idcs = to_indices(parent, I)
+    adim = findfirst(isa.(idcs, Array)) # which of the indices is the array
+    arr = idcs[adim] # the array itself
+    T = eltype(parent)
+    ret = Array{T}(undef, map(length, idcs))
+    Ilo = idcs[1:adim-1]
+    Ihi = idcs[adim+1:end]
+    for ii in eachindex(arr)
+        ret[Ilo..., ii, Ihi...] .= parent[Ilo..., arr[ii], Ihi...]
+    end
+    ret
 end
 
 function show(io::IO, o::HDF5Output)
diff --git a/src/Processing.jl b/src/Processing.jl
index 020b24ad..7ec50404 100644
--- a/src/Processing.jl
+++ b/src/Processing.jl
@@ -9,6 +9,7 @@ import Luna.Output: AbstractOutput, HDF5Output
 import Cubature: hcubature
 import ProgressLogging: @progress
 import Logging: @warn
+import DSP: unwrap
 
 """
     Common(val)
@@ -623,6 +624,41 @@ end
 fftnorm(grid::RealGrid) = Maths.rfftnorm(grid.t[2] - grid.t[1])
 fftnorm(grid::EnvGrid) = Maths.fftnorm(grid.t[2] - grid.t[1])
 
+
+"""
+    getφ(grid, Eω)
+    getφ(ω, Eω, τ)
+
+Extract the unwrapped spectral phase from the field `Eω`, subtracting the linear phase ramp corresponding
+to a pulse in the middle of the time window defined by the `grid`. 
+"""
+function getφ(grid::AbstractGrid, Eω)
+    ω = grid.ω
+    t = grid.t
+    τ = length(t) * (t[2] - t[1])/2 # middle of time window
+    getφ(ω, Eω, τ)
+end
+
+function getφ(ω::AbstractVector, Eω, τ)
+    φ = unwrap(angle.(Eω); dims=1)
+    φ .- ω*τ
+end
+
+"""
+    getφ(output, args...)
+
+Extract the frequency-domain `Eω` from the `output` (additional `args...` are passed to `getEω`) and
+extract the spectral phase, subtracting the linear phase ramp corresponding
+to a pulse in the middle of the time window defined by the frequency grid.
+"""
+function getφ(output, args...)
+    ω, Eω = getEω(output, args...)
+    grid = makegrid(output)
+    t = grid.t
+    τ = length(t) * (t[2] - t[1])/2 # middle of time window
+    getφ(ω, Eω, τ)
+end
+
 """
     getEt(output[, zslice]; kwargs...)
 
diff --git a/src/RK45.jl b/src/RK45.jl
index b4711ecc..2f60f9d8 100644
--- a/src/RK45.jl
+++ b/src/RK45.jl
@@ -47,13 +47,13 @@ function solve(s, tmax; stepfun=donothing!, output=false, outputN=201,
             speed = s.tn/(Dates.value(Dates.now()-start)/1000)
             eta_in_s = (tmax-s.tn)/(speed)
             if eta_in_s > 356400
-                Logging.@info @sprintf("Progress: %.2f %%, ETA: XX:XX:XX, stepsize %.2e, err %.2f, repeated %d",
-                s.tn/tmax*100, s.dt, s.err, repeated_tot)
+                Logging.@info @sprintf("Progress: %.2f %%, ETA: XX:XX:XX, stepsize %.2e, err %.2f, repeated %d/%d steps",
+                s.tn/tmax*100, s.dt, s.err, repeated_tot, steps)
             else
                 eta_in_ms = Dates.Millisecond(ceil(eta_in_s*1000))
                 etad = Dates.DateTime(Dates.UTInstant(eta_in_ms))
-                Logging.@info @sprintf("Progress: %.2f %%, ETA: %s, stepsize %.2e, err %.2f, repeated %d",
-                    s.tn/tmax*100, Dates.format(etad, "HH:MM:SS"), s.dt, s.err, repeated_tot)
+                Logging.@info @sprintf("Progress: %.2f %%, ETA: %s, stepsize %.2e, err %.2f, repeated %d/%d steps",
+                    s.tn/tmax*100, Dates.format(etad, "HH:MM:SS"), s.dt, s.err, repeated_tot, steps)
             end
             flush(stderr)
             tic = Dates.now()