diff --git a/src/components/gradient_waveforms/no_gradient_blocks.jl b/src/components/gradient_waveforms/no_gradient_blocks.jl
index 6beb655a33680d39610b3f85213c738d9faf7095..ce7667aed5a957e46f0ca8491129068edd22aef1 100644
--- a/src/components/gradient_waveforms/no_gradient_blocks.jl
+++ b/src/components/gradient_waveforms/no_gradient_blocks.jl
@@ -1,6 +1,6 @@
 module NoGradientBlocks
 import StaticArrays: SVector, SMatrix
-import ....Variables: VariableType, duration, qval, qvec, bmat_gradient, gradient_strength, slew_rate, get_free_variable
+import ....Variables: VariableType, duration, qval, bmat_gradient, gradient_strength, slew_rate, get_free_variable
 import ...AbstractTypes: GradientWaveform
 import ..ChangingGradientBlocks: split_gradient
 
@@ -31,9 +31,6 @@ end
 
 NoGradient(duration) = NoGradient{1}(duration)
 
-qvec(::NoGradient, index1, index2) = zero(SVector{3, Float64})
-qvec(::NoGradient) = zero(SVector{3, Float64})
-
 bmat_gradient(::NoGradient{1}) = 0.
 bmat_gradient(::NoGradient{3}) = zero(SMatrix{3, 3, Float64, 9})
 bmat_gradient(ngb::NoGradient{1}, qstart::VariableType) = qstart^2 * duration(ngb)
diff --git a/src/containers/building_blocks.jl b/src/containers/building_blocks.jl
index 437e99980f4a5a1ec110fc7afd96b471ab623c09..5076f1303edc5cb5b5a2e6e5c31007255912fc1e 100644
--- a/src/containers/building_blocks.jl
+++ b/src/containers/building_blocks.jl
@@ -7,7 +7,7 @@ import JuMP: @constraint
 import ..Abstract: ContainerBlock, start_time
 import ...BuildSequences: global_model
 import ...Components: BaseComponent, GradientWaveform, EventComponent, NoGradient, ChangingGradient, ConstantGradient, split_gradient, DelayedEvent, RFPulseComponent, ReadoutComponent
-import ...Variables: qval, bmat_gradient, effective_time, get_free_variable
+import ...Variables: qval, bmat_gradient, effective_time, get_free_variable, qval3
 import ...Variables: VariableType, duration, make_generic, get_pulse, get_readout, scanner_constraints!
 
 """
@@ -19,7 +19,7 @@ Main interface:
 - [`waveform_sequence`](@ref) returns just the gradient waveform as a sequence of [`GradientWaveform`](@ref) objects.
 - [`waveform`](@ref) returns just the gradient waveform as a sequence of (time, gradient_strength) tuples.
 - [`events`](@ref) returns the RF pulses and readouts.
-- [`qvec`](@ref) returns area under curve for (part of) the gradient waveform.
+- [`qval3`](@ref) returns area under curve for (part of) the gradient waveform.
 
 Sub-types need to implement:
 - `Base.keys`: returns sequence of keys to all the components.
@@ -183,7 +183,7 @@ Similarly, if `last_event` is set to something else than `nothing`, only the gra
 """
 function qval(bb::BaseBuildingBlock, index1, index2)
     if (!isnothing(index1)) && (index1 == index2)
-        return zeros(3)
+        return 0.
     end
     sum([qval(wv) for (_, wv) in waveform_sequence(bb, index1, index2)]; init=0.)
 end
@@ -198,7 +198,7 @@ function bmat_gradient(bb::BaseBuildingBlock, qstart, index1, index2)
 
     for (_, part) in waveform_sequence(bb, index1, index2)
         result = result .+ bmat_gradient(part, qcurrent)
-        qcurrent = qcurrent .+ qvec(part, qcurrent)
+        qcurrent = qcurrent .+ qval3(part, qcurrent)
     end
     return result
 end
diff --git a/src/parts/epi_readouts.jl b/src/parts/epi_readouts.jl
index 960e427d9ee1f624b89d790bc69d7672bdda96fa..acb4a32112c642c9060ba55b031dd4401603a9a3 100644
--- a/src/parts/epi_readouts.jl
+++ b/src/parts/epi_readouts.jl
@@ -2,7 +2,7 @@ module EPIReadouts
 import ...Containers: BaseSequence, get_index_single_TR
 import ..Trapezoids: Trapezoid, opposite_kspace_lines, LineReadout
 import ...Components: ADC
-import ...Variables: get_free_variable, VariableType, qval, qvec, set_simple_constraints!, resolution, inverse_voxel_size, inverse_fov, resolution, get_readout, apply_simple_constraint!, effective_time
+import ...Variables: get_free_variable, VariableType, qval, qval3, set_simple_constraints!, resolution, inverse_voxel_size, inverse_fov, resolution, get_readout, apply_simple_constraint!, effective_time
 import ...Pathways: PathwayWalker, update_walker_till_time!, walk_pathway!
 
 """
@@ -49,10 +49,10 @@ function EPIReadout(; resolution::AbstractVector{<:Integer}, recenter=false, gro
         get_free_variable(nothing),
         ky_lines
     )
-    apply_simple_constraint!(qvec(res.start_gradient), VariableType[-qval(pos)/2, ky_lines[1] * res.ky_step, 0.])
+    apply_simple_constraint!(qval3(res.start_gradient), VariableType[-qval(pos)/2, ky_lines[1] * res.ky_step, 0.])
     if recenter
         sign = isodd(length(ky_lines)) ? -1 : 1
-        apply_simple_constraint!(qvec(res.recenter_gradient), VariableType[sign * qval(pos)/2, -ky_lines[end] * res.ky_step, 0.])
+        apply_simple_constraint!(qval3(res.recenter_gradient), VariableType[sign * 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/pathways.jl b/src/pathways.jl
index 6fc3a77b69ae3463c5db9cccf6657b71557ee501..cc46223637f558f89c7c365337052f3e982dd469 100644
--- a/src/pathways.jl
+++ b/src/pathways.jl
@@ -3,7 +3,7 @@ import LinearAlgebra: norm, tr
 import StaticArrays: SVector, SMatrix
 import ..Components: NoGradient, RFPulseComponent, ReadoutComponent, InstantGradient, GradientWaveform, DelayedEvent
 import ..Containers: BaseSequence, Sequence, BaseBuildingBlock, waveform, events, waveform_sequence, start_time, AlternativeBlocks
-import ..Variables: qvec, qval, bmat_gradient, VariableType, effective_time, duration, TR, bmat, bval, area_under_curve, duration_dephase, duration_transverse
+import ..Variables: qvec, qval, qval3, bmat_gradient, VariableType, effective_time, duration, TR, bmat, bval, area_under_curve, duration_dephase, duration_transverse
 
 
 """
@@ -481,7 +481,7 @@ function update_walker_gradient!(gradient::GradientWaveform, walker::PathwayWalk
     # update qvec/bmat during gradient
     tracker = walker.gradient_trackers[key]
     tracker.bmat = tracker.bmat .+ bmat_gradient(gradient, tracker.qvec)
-    tracker.qvec = tracker.qvec .+ qvec(gradient)
+    tracker.qvec = tracker.qvec .+ qval3(gradient)
     tracker.last_gradient_time = gradient_start_time + duration(gradient)
 end
 
diff --git a/src/variables.jl b/src/variables.jl
index 1ce233bb762b6bc32cdeb5750303a356f93e3c87..3ebb8fe5773fe035d755e7b9efe45dd386621319 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -56,7 +56,7 @@ all_variables_symbols = [
         :inverse_slice_thickness => "Inverse of slice thickness of an RF pulse that is active during a gradient in 1/mm. Also, see [`slice_thickness`](@ref).",
     ],
     :gradient => [
-        :qvec => "The spatial range with orientation on which the displacements can be detected due to this gradient in rad/um (also see [`qval`](@ref)).",
+        :qval3 => "The spatial range with orientation on which the displacements can be detected due to this gradient in rad/um (also see [`qval`](@ref)).",
         :qval => "The spatial range on which the displacements can be detected due to this gradient in rad/um. This will be a scalar if the orientation is fixed and a scalar otherwise. If you always want a vector, use [`qvec`](@ref).",
         :qval_square => "Square of [`qval`](@ref) in rad^2/um^2.",
         :δ => "Effective duration of a gradient pulse ([`rise_time`](@ref) + [`flat_time`](@ref)) in ms.",
@@ -228,7 +228,7 @@ end
 
 for (target_name, all_vars) in all_variables_symbols
     for (variable_func, _) in all_vars
-        if variable_func in [:qvec, ]
+        if variable_func == :qval3
             continue
         end
         get_func = Symbol("get_" * string(target_name))
@@ -276,7 +276,7 @@ for (target_name, all_vars) in all_variables_symbols
 end
 
 
-function qvec(bb::AbstractBlock, args...; kwargs...)
+function qval3(bb::AbstractBlock, args...; kwargs...)
     value = qval(bb, args...; kwargs...)
     if value isa Number && iszero(value)
         return value
@@ -340,6 +340,7 @@ apply_simple_constraint!(variable, ::Val{:min}) = @objective global_model() Min
 apply_simple_constraint!(variable, ::Val{:max}) = @objective global_model() Min objective_function(global_model()) - variable
 apply_simple_constraint!(variable, value::VariableType) = @constraint global_model() variable == value
 apply_simple_constraint!(variable::AbstractVector, value::AbstractVector) = [apply_simple_constraint!(v1, v2) for (v1, v2) in zip(variable, value)]
+apply_simple_constraint!(variable::AbstractVector, value::VariableType) = [apply_simple_constraint!(v1, value) for v1 in variable]
 apply_simple_constraint!(variable::Number, value::Number) = @assert variable ≈ value "Variable set to multiple incompatible values."
 function apply_simple_constraint!(variable::NamedTuple, value::NamedTuple)
     for key in keys(value)