diff --git a/src/components/gradient_waveforms/changing_gradient_blocks.jl b/src/components/gradient_waveforms/changing_gradient_blocks.jl
index db276161f8db14ac95f3fb4abb79ebd4ac61d66f..57bd1590090cc553f5cd97b9ec3ea74f8864086a 100644
--- a/src/components/gradient_waveforms/changing_gradient_blocks.jl
+++ b/src/components/gradient_waveforms/changing_gradient_blocks.jl
@@ -34,16 +34,16 @@ struct ChangingGradient3D <: ChangingGradient{3}
     group :: Union{Nothing, Symbol}
 end
 
-@defvar begin
-    duration(cgb::ChangingGradient) = cgb.duration
-end
+@defvar duration(cgb::ChangingGradient) = cgb.duration
 
 @defvar gradient begin
-    grad_start(cgb::ChangingGradient) = cgb.gradient_strength_start
-    slew_rate(cgb::ChangingGradient) = cgb.slew_rate
+    grad_start(cgb::ChangingGradient1D) = cgb.gradient_strength_start .* cgb.orientation
+    grad_start(cgb::ChangingGradient3D) = cgb.gradient_strength_start
+    slew_rate(cgb::ChangingGradient1D) = cgb.slew_rate .* cgb.orientation
+    slew_rate(cgb::ChangingGradient3D) = cgb.slew_rate
     grad_end(cgb::ChangingGradient) = variables.grad_start(cgb) .+ variables.slew_rate(cgb) .* variables.duration(cgb)
     gradient_strength(cgb::ChangingGradient) = max.(variables.grad_start(cgb), variables.grad_end(cgb))
-    qval(cgb::ChangingGradient) = (variables.grad_start(cgb) .+ variables.grad_end(cgb)) .* (variables.duration(cgb) * π)
+    qvec(cgb::ChangingGradient) = (variables.grad_start(cgb) .+ variables.grad_end(cgb)) .* (variables.duration(cgb) * π)
 
     gradient_strength(cgb::ChangingGradient, time::Number) = variables.slew_rate(cgb) .* time .+ variables.grad_start(cgb)
 end
@@ -51,9 +51,6 @@ end
 _mult(g1::VariableType, g2::VariableType) = g1 * g2
 _mult(g1::AbstractVector, g2::AbstractVector) = g1 .* permutedims(g2)
 
-to_vec(cgb::ChangingGradient1D, g::VariableType) = cgb.orientation .* g
-to_vec(::ChangingGradient3D, g::AbstractVector) = g
-
 @defvar gradient begin
     function bmat_gradient(cgb::ChangingGradient, qstart::AbstractVector)
         # grad = (g1 * (duration - t) + g2 * t) / duration
@@ -66,7 +63,7 @@ to_vec(::ChangingGradient3D, g::AbstractVector) = g
         #   g1^2 * duration^3 / 3 +
         #   g1 * (g2 - g1) * duration^3 / 4 +
         #   (g2 - g1)^2 * duration^3 / 10
-        grad_aver = to_vec(cgb, 2 .* grad_start(cgb) .+ grad_end(cgb))
+        grad_aver = 2 .* grad_start(cgb) .+ grad_end(cgb)
         return (
             _mult(qstart, qstart) .* duration(cgb) .+
             duration(cgb)^2 .* _mult(qstart, grad_aver) .* 2Ï€ ./ 3 .+
@@ -75,8 +72,8 @@ to_vec(::ChangingGradient3D, g::AbstractVector) = g
     end
 
     function bmat_gradient(cgb::ChangingGradient)
-        gs = to_vec(cgb, grad_start(cgb))
-        diff = to_vec(cgb, slew_rate(cgb) .* duration(cgb))
+        gs = grad_start(cgb)
+        diff = slew_rate(cgb) .* duration(cgb)
         return (2Ï€)^2 .* (
             _mult(gs, gs) ./ 3 .+
             _mult(gs, diff) ./ 4 .+
diff --git a/src/components/gradient_waveforms/constant_gradient_blocks.jl b/src/components/gradient_waveforms/constant_gradient_blocks.jl
index 397b70368b42cb09f1f811030669f48adae2b4bd..32e06ba63df3010fadb5d05e34bed65577931f3b 100644
--- a/src/components/gradient_waveforms/constant_gradient_blocks.jl
+++ b/src/components/gradient_waveforms/constant_gradient_blocks.jl
@@ -34,11 +34,10 @@ end
 @defvar duration(cgb::ConstantGradient) = cgb.duration
 
 @defvar gradient begin
-    gradient_strength(cgb::ConstantGradient) = cgb.gradient_strength
-    slew_rate(::ConstantGradient1D) = 0.
-    slew_rate(::ConstantGradient3D) = zero(SVector{3, Float64})
-    qval(cgb::ConstantGradient1D) = variables.duration(cgb) * variables.gradient_strength(cgb) * 2Ï€
-    qval(cgb::ConstantGradient3D) = variables.duration(cgb) .* variables.gradient_strength(cgb) .* 2Ï€
+    gradient_strength(cgb::ConstantGradient1D) = cgb.gradient_strength * cgb.orientation
+    gradient_strength(cgb::ConstantGradient3D) = cgb.gradient_strength
+    slew_rate(::ConstantGradient) = zero(SVector{3, Float64})
+    qvec(cgb::ConstantGradient) = variables.duration(cgb) * variables.gradient_strength(cgb) * 2Ï€
 
     gradient_strength(cgb::ConstantGradient, time::Number) = variables.gradient_strength(cgb)
 end
@@ -46,12 +45,9 @@ end
 _mult(g1::VariableType, g2::VariableType) = g1 * g2
 _mult(g1::AbstractVector, g2::AbstractVector) = g1 .* permutedims(g2)
 
-to_vec(cgb::ConstantGradient1D, g::VariableType) = cgb.orientation .* g
-to_vec(::ConstantGradient3D, g::AbstractVector) = g
-
 @defvar begin
     function bmat_gradient(cgb::ConstantGradient)
-        grad = to_vec(cgb, 2Ï€ .* gradient_strength(cgb))
+        grad = 2Ï€ .* gradient_strength(cgb)
         return _mult(grad, grad) .* duration(cgb)^3 ./3
     end
 
@@ -60,7 +56,7 @@ to_vec(::ConstantGradient3D, g::AbstractVector) = g
         #   qstart^2 * duration +
         #   qstart * grad * duration^2 +
         #   grad * grad * duration^3 / 3 +
-        grad = to_vec(cgb, 2Ï€ .* gradient_strength(cgb))
+        grad = 2Ï€ .* gradient_strength(cgb)
         return (
             _mult(qstart, qstart) .* duration(cgb) .+
             _mult(qstart, grad) .* duration(cgb)^2 .+
diff --git a/src/components/gradient_waveforms/gradient_waveforms.jl b/src/components/gradient_waveforms/gradient_waveforms.jl
index 5e60f1d808ddea2da46b39b1a8f0817deb61a39a..2b9f53e33cc9cb3251438d7c016c3505574cb6a1 100644
--- a/src/components/gradient_waveforms/gradient_waveforms.jl
+++ b/src/components/gradient_waveforms/gradient_waveforms.jl
@@ -11,7 +11,7 @@ These parts are combined into a full gradient waveform in a `BuildingBlock`.
 Each part of this gradient waveform can compute:
 - `gradient_strength`: maximum gradient strength in each dimension.
 - `slew_rate`: maximum slew rate in each dimension.
-- `qval`/`qval3`: area under curve
+- `qval`/`qvec`: area under curve
 - `bmat_gradient`: diffusion weighting (scalar in 1D or matrix in 3D).
 """
 module GradientWaveforms
diff --git a/src/components/gradient_waveforms/no_gradient_blocks.jl b/src/components/gradient_waveforms/no_gradient_blocks.jl
index cd0dbcc8b7bc081ecdec547c4508efd094117c35..0dbe8461e5927ffc03d4b110e79e286dbee0b597 100644
--- a/src/components/gradient_waveforms/no_gradient_blocks.jl
+++ b/src/components/gradient_waveforms/no_gradient_blocks.jl
@@ -15,8 +15,8 @@ struct NoGradient <: GradientWaveform{0}
     duration :: VariableType
 end
 
-for func in (:qval, :gradient_strength, :slew_rate)
-    @eval variables.$(func).f(::NoGradient) = 0.
+for func in (:qvec, :gradient_strength, :slew_rate)
+    @eval variables.$(func).f(::NoGradient) = zero(SVector{3, Float64})
 end
 
 @defvar begin
diff --git a/src/components/instant_gradients.jl b/src/components/instant_gradients.jl
index 441a6cb0d79e5a6d0b784e18321aa545465a5dcd..68f6b745985c855ed0abbb32ff3a575ba9c3a851 100644
--- a/src/components/instant_gradients.jl
+++ b/src/components/instant_gradients.jl
@@ -42,7 +42,7 @@ struct InstantGradient1D <: InstantGradient{1}
     group :: Union{Nothing, Symbol}
 end
 
-@defvar gradient qval(ig::InstantGradient1D) = ig.qval
+@defvar gradient qvec(ig::InstantGradient1D) = ig.qval .* ig.orientation
 
 """
 An [`InstantGradient`](@ref) with a variable orientation.
@@ -61,7 +61,7 @@ function InstantGradient3D(; qval=[nothing, nothing, nothing], group=nothing, va
     return res
 end
 
-@defvar gradient qval(ig::InstantGradient3D) = ig.qvec
+@defvar gradient qvec(ig::InstantGradient3D) = ig.qvec
 
 
 variables.duration.f(::InstantGradient) = 0.
diff --git a/src/containers/building_blocks.jl b/src/containers/building_blocks.jl
index 3b67ff9104b29ef34ffc1fd4425d7454e071905d..c827705c33de6bf3239d208f1215e7667237a0a3 100644
--- a/src/containers/building_blocks.jl
+++ b/src/containers/building_blocks.jl
@@ -191,22 +191,22 @@ function waveform_sequence(bb::BaseBuildingBlock, first, last)
 end
 
 @defvar begin
-    function qval(bb::BaseBuildingBlock, index1, index2)
+    function qvec(bb::BaseBuildingBlock, index1, index2)
         if (!isnothing(index1)) && (index1 == index2)
             return 0.
         end
-        res = sum([variables.qval(wv) for (_, wv) in waveform_sequence(bb, index1, index2)])
+        res = sum([variables.qvec(wv) for (_, wv) in waveform_sequence(bb, index1, index2)])
 
         t1 = isnothing(index1) ? 0. : start_time(bb, index1)
         t2 = isnothing(index2) ? duration(bb) : start_time(bb, index2)
         for (key, event) in events(bb)
             if event isa InstantGradient && (t1 <= start_time(bb, key) <= t2)
-                res = res .+ variables.qval(event)
+                res = res .+ variables.qvec(event)
             end
         end
         return res
     end
-    qval(bb::BaseBuildingBlock) = variables.qval(bb, nothing, nothing)
+    qvec(bb::BaseBuildingBlock) = variables.qvec(bb, nothing, nothing)
 
     function bmat_gradient(bb::BaseBuildingBlock, qstart, index1, index2)
         if (!isnothing(index1)) && (index1 == index2)
@@ -217,7 +217,7 @@ end
 
         for (_, part) in waveform_sequence(bb, index1, index2)
             result = result .+ variables.bmat_gradient(part, qcurrent)
-            qcurrent = qcurrent .+ qval3(part, qcurrent)
+            qcurrent = qcurrent .+ variables.qvec(part, qcurrent)
         end
         return result
     end
@@ -225,14 +225,14 @@ end
 end
 
 """
-    qval(overlapping[, first_event, last_event])
+    qvec(overlapping[, first_event, last_event])
 
 Computes the area under the curve for the gradient waveform in [`BaseBuildingBlock`](@ref).
 
 If `first_event` is set to something else than `nothing`, only the gradient waveform after this RF pulse/Readout will be considered.
 Similarly, if `last_event` is set to something else than `nothing`, only the gradient waveform up to this RF pulse/Readout will be considered.
 """
-qval
+qvec
 
 function edge_times(bb::BaseBuildingBlock)
     res = Float64[]
diff --git a/src/containers/linearise.jl b/src/containers/linearise.jl
index ed8b5b92f7a1ebce1c3f0541c1d7f90d63f7a2be..3b06f6184980f487660689fab717099223141471 100644
--- a/src/containers/linearise.jl
+++ b/src/containers/linearise.jl
@@ -182,7 +182,7 @@ function LinearSequence(container::BaseSequence, times::AbstractVector{<:Number}
                     continue
                 end
                 if pulse isa InstantGradient1D
-                    pulse = InstantGradient3D(qval(pulse) .* pulse.orientation, pulse.group)
+                    pulse = InstantGradient3D(qvec(pulse) .* pulse.orientation, pulse.group)
                 end
                 index = findmin(t -> abs(t - real_time), times)[2] 
                 to_store[index] = pulse
diff --git a/src/parts/epi_readouts.jl b/src/parts/epi_readouts.jl
index de8d2e767d2f43f3c91823dd3e1db4b2b288ad85..e74d20d3b4690c37fe8cc6062999160f5c922feb 100644
--- a/src/parts/epi_readouts.jl
+++ b/src/parts/epi_readouts.jl
@@ -49,10 +49,10 @@ function EPIReadout(; resolution::AbstractVector{<:Integer}, recenter=false, gro
         get_free_variable(nothing),
         ky_lines
     )
-    apply_simple_constraint!(variables.qval(res.start_gradient), VariableType[-variables.qval(pos)/2, ky_lines[1] * res.ky_step, 0.])
+    apply_simple_constraint!(variables.qvec(res.start_gradient), VariableType[-variables.qval(pos)/2, ky_lines[1] * res.ky_step, 0.])
     if recenter
         sign = isodd(length(ky_lines)) ? -1 : 1
-        apply_simple_constraint!(variables.qval(res.recenter_gradient), VariableType[sign * variables.qval(pos)/2, -ky_lines[end] * res.ky_step, 0.])
+        apply_simple_constraint!(variables.qvec(res.recenter_gradient), VariableType[sign * variables.qval(pos)/2, -ky_lines[end] * res.ky_step, 0.])
     end
     for shift in unique(ky_lines[2:end] - ky_lines[1:end-1])
         res.blips[shift] = Trapezoid(orientation=[0, shift > 0 ? 1 : -1, 0], group=group, qval=abs(shift) * res.ky_step)
diff --git a/src/parts/helper_functions.jl b/src/parts/helper_functions.jl
index aac6a8903a43bd189633264001aa726030bbe08f..038e3a39ed756353e62427894a175830ba2b8ddd 100644
--- a/src/parts/helper_functions.jl
+++ b/src/parts/helper_functions.jl
@@ -215,10 +215,10 @@ function dwi_gradients(; type=:trapezoid, optimise=false, scanner=nothing, refoc
         (g1, g2) = [func_dict[type](;
             group=group, orientation=o, Dict(key => get_index(value, i) for (key, value) in pairs(real_variables))...
             ) for (i, o) in enumerate((orientation, other_orientation))]
-        if !isnothing(orientation) || refocus
-            apply_simple_constraint!(variables.qval(g1), variables.qval(g2))
+        if refocus
+            apply_simple_constraint!(variables.qvec(g1), variables.qvec(g2))
         else
-            apply_simple_constraint!(variables.qval(g1), -variables.qval(g2))
+            apply_simple_constraint!(variables.qvec(g1), -variables.qvec(g2))
         end
         for var_func in match
             if var_func isa Symbol
diff --git a/src/parts/slice_select_rephases.jl b/src/parts/slice_select_rephases.jl
index 94a7bfa19187491e136ee9cc110bee2b9afccf98..1def64ee05e9a263cebdf8a79729ba035ca6b07d 100644
--- a/src/parts/slice_select_rephases.jl
+++ b/src/parts/slice_select_rephases.jl
@@ -23,11 +23,7 @@ function SliceSelectRephase(pulse::RFPulseComponent; orientation=nothing, group=
         SliceSelect(pulse; orientation=orientation, group=group, slice_thickness=slice_thickness, kwargs...),
         Trapezoid(; orientation=rephase_orientation, group=group, kwargs...)
     )
-    if N == 1
-        apply_simple_constraint!(variables.qval(res.slice_select, :pulse, nothing), variables.qval(res.rephase))
-    else
-        apply_simple_constraint!(variables.qval(res.slice_select, :pulse, nothing), -variables.qval(res.rephase))
-    end
+    apply_simple_constraint!(variables.qvec(res.slice_select, :pulse, nothing), -variables.qvec(res.rephase))
     return res
 end
 
diff --git a/src/parts/spoilt_slice_selects.jl b/src/parts/spoilt_slice_selects.jl
index fd70b365ab916d80f7d34e289630c921a2fe796b..ea4d3b1876f70de971d9ddf1f42aab5396da412a 100644
--- a/src/parts/spoilt_slice_selects.jl
+++ b/src/parts/spoilt_slice_selects.jl
@@ -35,7 +35,7 @@ struct SpoiltSliceSelect <: BaseBuildingBlock
     slew_rate :: Number
 end
 
-function SpoiltSliceSelect(pulse::RFPulseComponent; orientation=[0, 0, 1], group=nothing, spoiler_scale=nothing, slice_thickness=nothing, kwargs...)
+function SpoiltSliceSelect(pulse::RFPulseComponent; orientation=[0, 0, 1], group=nothing, slice_thickness=nothing, kwargs...)
     model = global_model()
     res = nothing
     if slice_thickness isa Number && isinf(slice_thickness)
@@ -72,18 +72,9 @@ function SpoiltSliceSelect(pulse::RFPulseComponent; orientation=[0, 0, 1], group
         end
         @constraint model res.diff_time <= res.rise_time1
         @constraint model res.diff_time <= res.fall_time2
-        @constraint model qval(res, nothing, :pulse) == qval(res, :pulse, nothing)
+        @constraint model qvec(res, nothing, :pulse) == qvec(res, :pulse, nothing)
     end
 
-    if !isnothing(spoiler_scale)
-        if spoiler_scale == :min # note that spoiler_scale is inverse of qval
-            @objective model Min objective_function(model) - qval(res, nothing, 1)
-        elseif spoiler_scale == :max
-            @objective model Min objective_function(model) + qval(res, nothing, 1)
-        else
-            @constraint model qval(res, nothing, :pulse) >= 2Ï€ * 1e-3 / spoiler_scale
-        end
-    end
     set_simple_constraints!(res, kwargs)
 
     max_time = gradient_strength(global_scanner()) / res.slew_rate
diff --git a/src/parts/trapezoids.jl b/src/parts/trapezoids.jl
index aa32bf9e3a0fe58b9e10844a4ee1e3d8d18b795a..ffe6045d7afda1c21e6fd852a1c9d4a32f808353 100644
--- a/src/parts/trapezoids.jl
+++ b/src/parts/trapezoids.jl
@@ -96,9 +96,9 @@ end
 
 Base.keys(::Trapezoid) = (Val(:rise), Val(:flat), Val(:fall))
 
-Base.getindex(pg::BaseTrapezoid{N}, ::Val{:rise}) where {N} = ChangingGradient(N == 3 ? zeros(3) : 0., slew_rate(pg), gradient_orientation(pg), rise_time(pg), get_group(pg))
-Base.getindex(pg::BaseTrapezoid, ::Val{:flat}) = ConstantGradient(gradient_strength(pg), gradient_orientation(pg), flat_time(pg), get_group(pg))
-Base.getindex(pg::BaseTrapezoid, ::Val{:fall}) = ChangingGradient(gradient_strength(pg), -slew_rate(pg), gradient_orientation(pg), rise_time(pg), get_group(pg))
+Base.getindex(pg::BaseTrapezoid{N}, ::Val{:rise}) where {N} = ChangingGradient(zeros(3), slew_rate(pg), rise_time(pg), get_group(pg))
+Base.getindex(pg::BaseTrapezoid, ::Val{:flat}) = ConstantGradient(gradient_strength(pg), flat_time(pg), get_group(pg))
+Base.getindex(pg::BaseTrapezoid, ::Val{:fall}) = ChangingGradient(gradient_strength(pg), -slew_rate(pg), rise_time(pg), get_group(pg))
 gradient_orientation(::BaseTrapezoid{3}) = nothing
 gradient_orientation(pg::BaseTrapezoid{1}) = gradient_orientation(get_gradient(pg))
 gradient_orientation(pg::Trapezoid{1}) = pg.orientation
@@ -109,7 +109,8 @@ get_group(pg::BaseTrapezoid) = get_group(get_gradient(pg))
 @defvar gradient begin
     rise_time(pg::Trapezoid) = pg.rise_time
     flat_time(pg::Trapezoid) = pg.flat_time
-    slew_rate(g::Trapezoid) = g.slew_rate
+    slew_rate(g::Trapezoid1D) = g.slew_rate .* g.orientation
+    slew_rate(g::Trapezoid3D) = g.slew_rate
 end
 
 @defvar gradient begin
@@ -119,7 +120,7 @@ end
 
 @defvar duration(g::BaseTrapezoid) = 2 * rise_time(g) + flat_time(g)
 
-@defvar qval(g::BaseTrapezoid, ::Nothing, ::Nothing) = variables.δ(g) .* gradient_strength(g) .* 2π
+@defvar qvec(g::BaseTrapezoid, ::Nothing, ::Nothing) = variables.δ(g) .* variables.gradient_strength(g) .* 2π
 
 adjustable(::BaseTrapezoid) = :gradient
 
diff --git a/src/pathways.jl b/src/pathways.jl
index 7fd46b549031e204d6208635bb986b2c5d180f79..24cd3e838b0e36b2fa5c9246dd25780c16d58935 100644
--- a/src/pathways.jl
+++ b/src/pathways.jl
@@ -473,7 +473,7 @@ function update_walker_gradient!(gradient::GradientWaveform, walker::PathwayWalk
         # update qvec/bmat during gradient
         tracker = walker.gradient_trackers[key]
         tracker.bmat = tracker.bmat .+ variables.bmat_gradient(gradient, tracker.qvec)
-        tracker.qvec = tracker.qvec .+ variables.qval3(gradient)
+        tracker.qvec = tracker.qvec .+ variables.qvec(gradient)
         tracker.last_gradient_time = gradient_start_time + variables.duration(gradient)
     end
 end
@@ -482,15 +482,11 @@ end
     update_walker_instant_gradient!(gradient, walker)
 """
 function update_walker_instant_gradient!(gradient::InstantGradient{N}, walker::PathwayWalker, gradient_start_time::VariableType) where {N}
-    if N == 1
-        qvec3 = qval(gradient) .* gradient.orientation
-    else
-        qvec3 = qval(gradient)
-    end
+    qvec3 = variables.qvec(gradient)
     for key in (isnothing(gradient.group) ? [nothing] : [nothing, gradient.group])
         update_gradient_tracker_till_time!(walker, key, gradient_start_time)
         tracker = walker.gradient_trackers[key]
-        tracker.qvec = tracker.qvec  .+ qvec3
+        tracker.qvec = tracker.qvec .+ qvec3
     end
 end
 
diff --git a/src/variables.jl b/src/variables.jl
index b324470e76e49e47fbe36b47ab51e92eb7b301c2..c5d9c2e8be1e63e342e5bf269c1bd78077efd037 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -74,8 +74,8 @@ end
 struct AlternateVariable <: AnyVariable
     name :: Symbol
     other_var :: Symbol
-    to_other :: Function
     from_other :: Function
+    to_other :: Union{Nothing, Function}
     inverse :: Bool
 end
 
@@ -203,11 +203,20 @@ Duration of the sequence or building block in ms.
 duration
 
 
-function def_alternate_variable!(name::Symbol, other_var::Symbol, to_other::Function, from_other::Function, inverse::Bool)
-    getfield(variables, :variables)[name] = AlternateVariable(name, other_var, to_other, from_other, inverse)
+function def_alternate_variable!(name::Symbol, other_var::Symbol, from_other::Function, to_other::Union{Nothing, Function}, inverse::Bool)
+    getfield(variables, :variables)[name] = AlternateVariable(name, other_var, from_other, to_other, inverse)
 end
 
-def_alternate_variable!(:spoiler_scale, :spoiler_scale, q->1e-3 * 2Ï€/q, l->1e-3 * 2Ï€/l, true)
+def_alternate_variable!(:spoiler_scale, :qval, q->1e-3 * 2Ï€/q, l->1e-3 * 2Ï€/l, true)
+def_alternate_variable!(:qval, :qval_square, sqrt, q -> q * q, false)
+def_alternate_variable!(:qval_square, :qvec, qv -> sum(q -> q * q, qv), nothing, false)
+
+for vec_variable in [:gradient_strength, :slew_rate]
+    vec_square = Symbol(string(vec_variable) * "_square")
+    vec_norm = Symbol(string(vec_variable) * "_norm")
+    def_alternate_variable!(:vec_norm, :vec_square, sqrt, v -> v * v, false)
+    def_alternate_variable!(:vec_square, :vec_variable, v -> v[1] * v[1] + v[2] * v[2] + v[3] * v[3], nothing, false)
+end
 
 for name in [:slice_thickness, :bandwidth, :fov, :voxel_size]
     inv_name = Symbol("inverse_" * string(name))
@@ -333,8 +342,8 @@ end
 
 function (var::AlternateVariable)(args...; kwargs...)
     other_var = variables[var.other_var]
-    apply_from_other(res::Number) = var.from_other(res)
-    apply_from_other(res::AbstractArray{<:Number}) = var.from_other.(res)
+    apply_from_other(res::VariableType) = var.from_other(res)
+    apply_from_other(res::AbstractArray{<:VariableType}) = var.from_other.(res)
     apply_from_other(res::NamedTuple) = NamedTuple(k => apply_from_other(v) for (k, v) in pairs(res))
     return apply_from_other(other_var(args...; kwargs...))
 end
diff --git a/test/test_post_hoc.jl b/test/test_post_hoc.jl
index 8233ea28fabbccd01baaee6b7865f4337d963a02..a396dbb31826e032b10ae280d292b03185716650 100644
--- a/test/test_post_hoc.jl
+++ b/test/test_post_hoc.jl
@@ -6,37 +6,37 @@
             dwi = DWI(bval=1., TE=:min)
             @test bval(dwi) ≈ 1.
             qval_orig = qval(dwi[:gradient])
-            @test all(qval3(dwi[:gradient]) .≈ [qval_orig, 0., 0.])
+            @test all(qvec(dwi[:gradient]) .≈ [qval_orig, 0., 0.])
 
             @testset "scale and change orientation" begin
                 new_dwi = adjust(dwi, diffusion=(scale=0.5, orientation=[0., 1., 0.]))
                 @test bval(new_dwi) ≈ 0.25
-                @test all(qval3(new_dwi[:gradient]) .≈ [0., qval_orig/2, 0.])
+                @test all(qvec(new_dwi[:gradient]) .≈ [0., qval_orig/2, 0.])
             end
             @testset "Rotate gradient" begin
                 new_dwi = adjust(dwi, gradient=(rotation=RotationVec(0., 0., π/4), ))
                 @test bval(new_dwi) ≈ 1.
-                @test all(qval3(new_dwi[:gradient]) .≈ [qval_orig/√2, qval_orig/√2, 0.])
+                @test all(qvec(new_dwi[:gradient]) .≈ [qval_orig/√2, qval_orig/√2, 0.])
             end
         end
         @testset "instant gradients" begin
             dwi = DWI(bval=1., TE=80, Δ=40, gradient=(type=:instant, ))
             @test bval(dwi) ≈ 1.
             qval_orig = qval(dwi[:gradient])
-            @test all(qval3(dwi[:gradient]) .≈ [qval_orig, 0., 0.])
+            @test all(qvec(dwi[:gradient]) .≈ [qval_orig, 0., 0.])
 
             @testset "scale and change orientation" begin
                 new_dwi = adjust(dwi, diffusion=(scale=0.5, orientation=[0., 1., 0.]))
                 @test bval(new_dwi) ≈ 0.25
-                @test all(qval3(new_dwi[:gradient]) .≈ [0., qval_orig/2, 0.])
+                @test all(qvec(new_dwi[:gradient]) .≈ [0., qval_orig/2, 0.])
 
                 @testset "multiple adjustments" begin
                     new_dwi = adjust(dwi, diffusion=(scale=[0.5, 1.], orientation=[0., 1., 0.]), merge=false)
                     @test length(new_dwi) == 2
                     @test bval(new_dwi[1]) ≈ 0.25
                     @test bval(new_dwi[2]) ≈ 1.
-                    @test all(qval3(new_dwi[1][:gradient]) .≈ [0., qval_orig/2, 0.])
-                    @test all(qval3(new_dwi[2][:gradient]) .≈ [0., qval_orig, 0.])
+                    @test all(qvec(new_dwi[1][:gradient]) .≈ [0., qval_orig/2, 0.])
+                    @test all(qvec(new_dwi[2][:gradient]) .≈ [0., qval_orig, 0.])
                     
                     new_dwi = adjust(dwi, diffusion=(scale=[0.5, 1.], orientation=[0., 1., 0.]))
                     @test duration(new_dwi) ≈ 160
@@ -50,7 +50,7 @@
             @testset "Rotate gradient" begin
                 new_dwi = adjust(dwi, gradient=(rotation=RotationVec(0., 0., π/4), ))
                 @test bval(new_dwi) ≈ 1.
-                @test all(qval3(new_dwi[:gradient]) .≈ [qval_orig/√2, qval_orig/√2, 0.])
+                @test all(qvec(new_dwi[:gradient]) .≈ [qval_orig/√2, qval_orig/√2, 0.])
 
             end
         end