diff --git a/docs/src/sequence_optimisation.md b/docs/src/sequence_optimisation.md
index d42899b9a7471df541defe47ba5ec21267d73f51..4a1b91855b933c0e1e2db27d3af336b27aac8f37 100644
--- a/docs/src/sequence_optimisation.md
+++ b/docs/src/sequence_optimisation.md
@@ -12,9 +12,15 @@ Internally, MRIBuilder will then optimise the `BuildingBlock` free parameters to
 This optimisation uses the [Ipopt](https://github.com/coin-or/Ipopt) optimiser accessed through the [JuMP.jl](https://jump.dev/JuMP.jl/stable/) library.
 
 ## [Summary variables](@id variables)
+All variables are available as members of the [`variables`](@ref) structure.
 ```@meta
 CollapsedDocStrings = true
 ```
+```@docs
+variables.
+```
+
+## Variables interface
 ```@autodocs
 Modules = [
     MRIBuilder.Variables,
diff --git a/ext/MakieMRIBuilder/MakieMRIBuilder.jl b/ext/MakieMRIBuilder/MakieMRIBuilder.jl
index e6f4a0bc07fd07cf59c4740ff20ab7eff9c7db05..f84ecd3400aa45236703a2fb5ae2093debcdd21f 100644
--- a/ext/MakieMRIBuilder/MakieMRIBuilder.jl
+++ b/ext/MakieMRIBuilder/MakieMRIBuilder.jl
@@ -3,7 +3,7 @@ using Makie
 import MakieCore
 import MakieCore: @recipe, theme, generic_plot_attributes!, Attributes, automatic
 import MRIBuilder.Plot: SequenceDiagram, range_full, normalise, plot_sequence
-import MRIBuilder: BaseSequence, BaseBuildingBlock, duration
+import MRIBuilder: BaseSequence, BaseBuildingBlock, variables
 
 @recipe(Plot_Sequence, sequence) do scene
     attr = Attributes(
@@ -53,7 +53,7 @@ function Makie.plot!(scene:: Plot_Sequence)
                 current_y += (upper - lower) + 0.1
             end
         end
-        Makie.lines!(scene, [-duration(sequence) / 10., 0], [current_y/2, current_y/2], color=(:red, 0.))
+        Makie.lines!(scene, [-variables.duration(sequence) / 10., 0], [current_y/2, current_y/2], color=(:red, 0.))
     end
     
 end
diff --git a/src/MRIBuilder.jl b/src/MRIBuilder.jl
index b00d3db134e98b47bef5b311f5f138d02b0a62f3..591f342628e04b4b530b67d025ef04af49b861d6 100644
--- a/src/MRIBuilder.jl
+++ b/src/MRIBuilder.jl
@@ -22,17 +22,17 @@ export build_sequence, global_model, global_scanner, fixed
 import .Scanners: Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra, Default_Scanner
 export Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra, Default_Scanner
 
-import .Variables: variables, effective_time, make_generic, @defvar, duration, get_pulse, get_readout, get_pathway, get_gradient
-export variables, effective_time, make_generic, @defvar, duration, get_pulse, get_readout, get_pathway, get_gradient
+import .Variables: variables, make_generic, @defvar, get_pulse, get_readout, get_pathway, get_gradient
+export variables, make_generic, @defvar, get_pulse, get_readout, get_pathway, get_gradient
 
 import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times
 export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times
 
-import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, iter_blocks, iter_instant_gradients, iter_instant_pulses
-export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, iter_blocks, iter_instant_gradients, iter_instant_pulses
+import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, iter_blocks, iter_instant_gradients, iter_instant_pulses
+export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, iter_blocks, iter_instant_gradients, iter_instant_pulses
 
-import .Pathways: Pathway, duration_transverse, duration_dephase, bval, bmat
-export Pathway, duration_transverse, duration_dephase, bval, bmat
+import .Pathways: Pathway
+export Pathway
 
 import .Parts: dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size
 export dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size
diff --git a/src/components/abstract_types.jl b/src/components/abstract_types.jl
index ab311463cb2e523cf37cf7ed79207142ce78f4e6..b57d1611f3b55bd05ce9d2c88308ff402a53c4c2 100644
--- a/src/components/abstract_types.jl
+++ b/src/components/abstract_types.jl
@@ -1,5 +1,5 @@
 module AbstractTypes
-import ...Variables: AbstractBlock, duration, variables, adjustable, gradient_orientation
+import ...Variables: AbstractBlock, variables, adjustable, gradient_orientation
 
 """
 Super-type for all individual components that form an MRI sequence (i.e., RF pulse, gradient waveform, or readout event).
@@ -67,6 +67,6 @@ Edges are defined as any time, when:
 
 Edges that are within `tol` ms of each other are considered to be one edge (default: 1 ns).
 """
-edge_times(comp::BaseComponent) = [0., duration(comp)]
+edge_times(comp::BaseComponent) = [0., variables.duration(comp)]
 
 end
\ No newline at end of file
diff --git a/src/components/gradient_waveforms/changing_gradient_blocks.jl b/src/components/gradient_waveforms/changing_gradient_blocks.jl
index 57bd1590090cc553f5cd97b9ec3ea74f8864086a..f670fff3e5dfb65af46f0ed38e9530aafee8e7eb 100644
--- a/src/components/gradient_waveforms/changing_gradient_blocks.jl
+++ b/src/components/gradient_waveforms/changing_gradient_blocks.jl
@@ -63,22 +63,22 @@ _mult(g1::AbstractVector, g2::AbstractVector) = g1 .* permutedims(g2)
         #   g1^2 * duration^3 / 3 +
         #   g1 * (g2 - g1) * duration^3 / 4 +
         #   (g2 - g1)^2 * duration^3 / 10
-        grad_aver = 2 .* grad_start(cgb) .+ grad_end(cgb)
+        grad_aver = 2 .* variables.grad_start(cgb) .+ variables.grad_end(cgb)
         return (
-            _mult(qstart, qstart) .* duration(cgb) .+
-            duration(cgb)^2 .* _mult(qstart, grad_aver) .* 2Ï€ ./ 3 .+
+            _mult(qstart, qstart) .* variables.duration(cgb) .+
+            variables.duration(cgb)^2 .* _mult(qstart, grad_aver) .* 2Ï€ ./ 3 .+
             variables.bmat_gradient(cgb)
         )
     end
 
     function bmat_gradient(cgb::ChangingGradient)
-        gs = grad_start(cgb)
-        diff = slew_rate(cgb) .* duration(cgb)
+        gs = variables.grad_start(cgb)
+        diff = variables.slew_rate(cgb) .* variables.duration(cgb)
         return (2Ï€)^2 .* (
             _mult(gs, gs) ./ 3 .+
             _mult(gs, diff) ./ 4 .+
             _mult(diff, diff) ./ 10
-        ) .* duration(cgb)^3
+        ) .* variables.duration(cgb)^3
     end
 end
 
@@ -95,7 +95,7 @@ For N times this returns a vector with the N+1 replacement [`ConstantGradient`](
 """
 function split_gradient(cgb::ChangingGradient, times::VariableType...)
     all_times = [0., times...]
-    durations = [times[1], [t[2] - t[1] for t in zip(times[1:end-1], times[2:end])]..., duration(cgb) - times[end]]
+    durations = [times[1], [t[2] - t[1] for t in zip(times[1:end-1], times[2:end])]..., variables.duration(cgb) - times[end]]
     if cgb isa ChangingGradient1D
         return [ChangingGradient1D(cgb.gradient_strength .+ cgb.slew_rate .* t, cgb.slew_rate, cgb.orientation, d, cgb.group) for (t, d) in zip(all_times, durations)]
     else
diff --git a/src/components/gradient_waveforms/constant_gradient_blocks.jl b/src/components/gradient_waveforms/constant_gradient_blocks.jl
index 32e06ba63df3010fadb5d05e34bed65577931f3b..7e6252fcd3a20378b554b4e648f0e60d2e1b98a1 100644
--- a/src/components/gradient_waveforms/constant_gradient_blocks.jl
+++ b/src/components/gradient_waveforms/constant_gradient_blocks.jl
@@ -47,8 +47,8 @@ _mult(g1::AbstractVector, g2::AbstractVector) = g1 .* permutedims(g2)
 
 @defvar begin
     function bmat_gradient(cgb::ConstantGradient)
-        grad = 2Ï€ .* gradient_strength(cgb)
-        return _mult(grad, grad) .* duration(cgb)^3 ./3
+        grad = 2Ï€ .* variables.gradient_strength(cgb)
+        return _mult(grad, grad) .* variables.duration(cgb)^3 ./3
     end
 
     function bmat_gradient(cgb::ConstantGradient, qstart::AbstractVector)
@@ -56,10 +56,10 @@ _mult(g1::AbstractVector, g2::AbstractVector) = g1 .* permutedims(g2)
         #   qstart^2 * duration +
         #   qstart * grad * duration^2 +
         #   grad * grad * duration^3 / 3 +
-        grad = 2Ï€ .* gradient_strength(cgb)
+        grad = 2Ï€ .* variables.gradient_strength(cgb)
         return (
-            _mult(qstart, qstart) .* duration(cgb) .+
-            _mult(qstart, grad) .* duration(cgb)^2 .+
+            _mult(qstart, qstart) .* variables.duration(cgb) .+
+            _mult(qstart, grad) .* variables.duration(cgb)^2 .+
             variables.bmat_gradient(cgb)
         )
     end
diff --git a/src/components/gradient_waveforms/no_gradient_blocks.jl b/src/components/gradient_waveforms/no_gradient_blocks.jl
index 0dbe8461e5927ffc03d4b110e79e286dbee0b597..4ae4f6e4630aa65be52ff1da9bbc85f1158ee7ea 100644
--- a/src/components/gradient_waveforms/no_gradient_blocks.jl
+++ b/src/components/gradient_waveforms/no_gradient_blocks.jl
@@ -26,12 +26,12 @@ end
 
 @defvar gradient begin
     bmat_gradient(::NoGradient) = 0.
-    bmat_gradient(ngb::NoGradient, qstart::VariableType) = qstart^2 * duration(ngb)
-    bmat_gradient(ngb::NoGradient, qstart::AbstractVector{<:VariableType}) = @. qstart * permutedims(qstart) * duration(ngb)
+    bmat_gradient(ngb::NoGradient, qstart::VariableType) = qstart^2 * variables.duration(ngb)
+    bmat_gradient(ngb::NoGradient, qstart::AbstractVector{<:VariableType}) = @. qstart * permutedims(qstart) * variables.duration(ngb)
 end
 
 function split_gradient(ngb::NoGradient, times::VariableType...)
-    durations = [times[1], [t[2] - t[1] for t in zip(times[1:end-1], times[2:end])]..., duration(ngb) - times[end]]
+    durations = [times[1], [t[2] - t[1] for t in zip(times[1:end-1], times[2:end])]..., variables.duration(ngb) - times[end]]
     return [NoGradient(d) for d in durations]
 end
 
diff --git a/src/components/pulses/composite_pulses.jl b/src/components/pulses/composite_pulses.jl
index 459c77f0dc466417395461d942c42ac1c5f8405f..6a3c2a96c496fe4d3da825bb9559bd7834b68ba4 100644
--- a/src/components/pulses/composite_pulses.jl
+++ b/src/components/pulses/composite_pulses.jl
@@ -76,13 +76,13 @@ end
 
 
 function get_pulse_index(comp::CompositePulse, time::Number)
-    t_first_center = 0.5 * duration(comp.pulses[1])
+    t_first_center = 0.5 * variables.duration(comp.pulses[1])
     index = Int(div(time - t_first_center, comp.pulse_time, RoundNearest)) + 1
     if index < 1 || index > length(comp)
         return (index, 0.)
     else
         relative_to_center = time - t_first_center - (index - 1) * comp.pulse_time
-        relative_to_start = relative_to_center + 0.5 * duration(comp.pulses[index])
+        relative_to_start = relative_to_center + 0.5 * variables.duration(comp.pulses[index])
         return (index, relative_to_start)
     end
 end
@@ -96,7 +96,7 @@ for (fn, default) in [
         (index, rtime) = get_pulse_index(pulse, time)
         if (
             index < 1 || index > length(pulse) ||
-            rtime < 0 || rtime > duration(pulse.pulses[index])
+            rtime < 0 || rtime > variables.duration(pulse.pulses[index])
         )
             return $default
         end
@@ -106,9 +106,9 @@ end
 
 function edge_times(comp::CompositePulse)
     res = Float64[]
-    t1 = duration(comp.pulses[1]) * 0.5
+    t1 = variables.duration(comp.pulses[1]) * 0.5
     for (index, pulse) in enumerate(comp.pulses)
-        d = duration(pulse)
+        d = variables.duration(pulse)
         push!(res, t1 - d/2 + comp.pulse_time * (index - 1))
         push!(res, t1 + d/2 + comp.pulse_time * (index - 1))
     end
@@ -118,12 +118,12 @@ end
 
 function make_generic(comp::CompositePulse)
     sub_generics = make_generic.(comp.pulses)
-    t1 = duration(comp.pulses[1]) * 0.5
+    t1 = variables.duration(comp.pulses[1]) * 0.5
     times = Float64[]
     amplitude = Float64[]
     phase = Float64[]
     for (index, generic) in enumerate(sub_generics)
-        start_time = t1 + (index - 1) * comp.pulse_time - duration(comp.pulses[index]) / 2
+        start_time = t1 + (index - 1) * comp.pulse_time - variables.duration(comp.pulses[index]) / 2
         push!(times, start_time)
         push!(amplitude, 0.)
         push!(phase, generic.phase[1])
diff --git a/src/components/pulses/constant_pulses.jl b/src/components/pulses/constant_pulses.jl
index ae9b7a05671712b2659f8bad0201e0b6e8cba5b2..9ba22d7e8403646cb3681d33d59b981f6919571b 100644
--- a/src/components/pulses/constant_pulses.jl
+++ b/src/components/pulses/constant_pulses.jl
@@ -39,7 +39,7 @@ function ConstantPulse(; amplitude=nothing, duration=nothing, phase=nothing, fre
 end
 
 @defvar duration(pulse::ConstantPulse) = pulse.duration
-@defvar effective_time(pulse::ConstantPulse) = duration(pulse) / 2
+@defvar effective_time(pulse::ConstantPulse) = variables.duration(pulse) / 2
 
 @defvar pulse begin
     amplitude(pulse::ConstantPulse) = pulse.amplitude
@@ -49,21 +49,21 @@ end
 end
 
 @defvar pulse begin
-    flip_angle(pulse::ConstantPulse) = amplitude(pulse) * duration(pulse) * 360
-    inverse_bandwidth(pulse::ConstantPulse) = duration(pulse) * 4
+    flip_angle(pulse::ConstantPulse) = variables.amplitude(pulse) * variables.duration(pulse) * 360
+    inverse_bandwidth(pulse::ConstantPulse) = variables.duration(pulse) * 4
 
     phase(pulse::ConstantPulse, time::Number) = variables.phase(pulse) + variables.frequency(pulse) * (time - variables.effective_time(pulse)) * 360.
     frequency(pulse::ConstantPulse, time::Number) = variables.frequency(pulse)
 end
 
 function make_generic(block::ConstantPulse)
-    d = duration(block)
-    final_phase = phase(block) + d * frequency(block) * 360
+    d = variables.duration(block)
+    final_phase = variables.phase(block) + d * variables.frequency(block) * 360
     return GenericPulse(
         [0., d], 
-        [amplitude(block), amplitude(block)],
-        [phase(block), final_phase],
-        effective_time(block)
+        [variables.amplitude(block), variables.amplitude(block)],
+        [variables.phase(block), final_phase],
+        variables.effective_time(block)
     )
 end
 
diff --git a/src/components/pulses/sinc_pulses.jl b/src/components/pulses/sinc_pulses.jl
index 41e9514dbf5899c126e8cb6af2360690f385d511..74efad56b88ebb8746245314b667a599faa55621 100644
--- a/src/components/pulses/sinc_pulses.jl
+++ b/src/components/pulses/sinc_pulses.jl
@@ -92,28 +92,28 @@ end
 end
 
 @defvar begin
-    duration(pulse::SincPulse) = (N_left(pulse) + N_right(pulse)) * lobe_duration(pulse)
-    flip_angle(pulse::SincPulse) = (pulse.norm_flip_angle[1] + pulse.norm_flip_angle[2]) * amplitude(pulse) * lobe_duration(pulse) * 360
-    inverse_bandwidth(pulse::SincPulse) = lobe_duration(pulse)
-    effective_time(pulse::SincPulse) = N_left(pulse) * lobe_duration(pulse)
+    duration(pulse::SincPulse) = (variables.N_left(pulse) + variables.N_right(pulse)) * variables.lobe_duration(pulse)
+    flip_angle(pulse::SincPulse) = (pulse.norm_flip_angle[1] + pulse.norm_flip_angle[2]) * variables.amplitude(pulse) * variables.lobe_duration(pulse) * 360
+    inverse_bandwidth(pulse::SincPulse) = variables.lobe_duration(pulse)
+    effective_time(pulse::SincPulse) = variables.N_left(pulse) * variables.lobe_duration(pulse)
 end
 
 @defvar begin
-    amplitude(pulse::SincPulse, time::Number) = variables.amplitude(pulse) * normalised_function(abs((time - effective_time(pulse))) / lobe_duration(pulse), N_left(pulse), N_right(pulse); apodise=pulse.apodise)
-    phase(pulse::SincPulse, time::Number) = variables.phase(pulse) + variables.frequency(pulse) * (time - effective_time(pulse)) * 360.
+    amplitude(pulse::SincPulse, time::Number) = variables.amplitude(pulse) * normalised_function(abs((time - variables.effective_time(pulse))) / variables.lobe_duration(pulse), N_left(pulse), N_right(pulse); apodise=pulse.apodise)
+    phase(pulse::SincPulse, time::Number) = variables.phase(pulse) + variables.frequency(pulse) * (time - variables.effective_time(pulse)) * 360.
     frequency(pulse::SincPulse, time::Number) = variables.frequency(pulse)
 end
 
 function make_generic(block::SincPulse)
-    normed_times = -N_left(block):0.1:N_right(block) + 1e-5
-    times = max.(0., (normed_times .+ N_left(block))) .* lobe_duration(block)
-    amplitudes = amplitude(block) .* (normalised_function.(normed_times, N_left(block), N_right(block); apodise=block.apodise))
-    phases = [frequency(block) .* lobe_duration(block)] .* normed_times .* 360
-    return GenericPulse(times, amplitudes, phases, effective_time(block))
+    normed_times = -variables.N_left(block):0.1:variables.N_right(block) + 1e-5
+    times = max.(0., (normed_times .+ variables.N_left(block))) .* variables.lobe_duration(block)
+    amplitudes = variables.amplitude(block) .* (normalised_function.(normed_times, variables.N_left(block), variables.N_right(block); apodise=block.apodise))
+    phases = [variables.frequency(block) .* variables.lobe_duration(block)] .* normed_times .* 360
+    return GenericPulse(times, amplitudes, phases, variables.effective_time(block))
 end
 
 function split_timestep(block::SincPulse, precision)
-    max_second_derivative = π^2/3 * (amplitude(block) / lobe_duration(block))^2
+    max_second_derivative = π^2/3 * (variables.amplitude(block) / variables.lobe_duration(block))^2
     return sqrt(2 * precision / max_second_derivative)
 end
 
diff --git a/src/components/readouts/ADCs.jl b/src/components/readouts/ADCs.jl
index 430fc4999e014a1fa4bbd874f4693b01858cf80c..0cfd1e512dc047aadc391e1946a615b6f93e21ac 100644
--- a/src/components/readouts/ADCs.jl
+++ b/src/components/readouts/ADCs.jl
@@ -39,10 +39,10 @@ function ADC(; resolution=nothing, dwell_time=nothing, time_to_center=nothing, c
     @constraint global_model() res.dwell_time >= 0
     @constraint global_model() res.oversample >= 1
     if center_halfway
-        apply_simple_constraint!(duration(res), 2 * res.time_to_center)
+        apply_simple_constraint!(variables.duration(res), 2 * res.time_to_center)
     else
         @constraint global_model() res.time_to_center >= 0
-        @constraint global_model() res.time_to_center <= duration(res)
+        @constraint global_model() res.time_to_center <= variables.duration(res)
     end
     set_simple_constraints!(res, kwargs)
     return res
@@ -55,18 +55,18 @@ end
     resolution(adc::ADC) = adc.resolution
 end
 
-@defvar readout nsamples(adc::ADC) = resolution(adc) * oversample(adc)
+@defvar readout nsamples(adc::ADC) = variables.resolution(adc) * variables.oversample(adc)
 
-@defvar readout_times(adc::ADC) = ((1:Int(nsamples(adc))) .- 0.5) .* dwell_time(adc)
+@defvar readout_times(adc::ADC) = ((1:Int(variables.nsamples(adc))) .- 0.5) .* variables.dwell_time(adc)
 @defvar begin
-    duration(adc::ADC) = nsamples(adc) * dwell_time(adc)
-    effective_time(adc::ADC) = time_to_center(adc)
+    duration(adc::ADC) = variables.nsamples(adc) * variables.dwell_time(adc)
+    effective_time(adc::ADC) = variables.time_to_center(adc)
 end
 
 function fixed(adc::ADC)
     # round nsamples during fixing
-    r = Int(round(value(resolution(adc)), RoundNearest))
-    n = Int(round(value(nsamples(adc)), RoundNearest))
+    r = Int(round(value(variables.resolution(adc)), RoundNearest))
+    n = Int(round(value(variables.nsamples(adc)), RoundNearest))
     if iszero(n)
         return ADC(0, NaN, NaN, NaN)
     end
diff --git a/src/components/readouts/readouts.jl b/src/components/readouts/readouts.jl
index b229e4c32c4c2e9a93dabfec880efa67fd74a039..32f13a22549c28ab0b7efe99c54321a4a6a135e2 100644
--- a/src/components/readouts/readouts.jl
+++ b/src/components/readouts/readouts.jl
@@ -3,7 +3,7 @@ include("ADCs.jl")
 include("single_readouts.jl")
 
 import ..AbstractTypes: ReadoutComponent, split_timestep
-import .ADCs: ADC, readout_times
+import .ADCs: ADC
 import .SingleReadouts: SingleReadout
 
 split_timestep(rc::ReadoutComponent, precision) = Inf
diff --git a/src/containers/abstract.jl b/src/containers/abstract.jl
index bf18ab93a318d74a6bd94ec2ffb70f017e3c9263..16f66ff2bc7898e266657055777b0a6c6c8b2913 100644
--- a/src/containers/abstract.jl
+++ b/src/containers/abstract.jl
@@ -1,5 +1,5 @@
 module Abstract
-import ...Variables: AbstractBlock, variables, VariableType, get_pulse, get_gradient, @defvar, duration
+import ...Variables: AbstractBlock, variables, VariableType, get_pulse, get_gradient, @defvar
 import ...Components: BaseComponent, InstantPulse, InstantGradient, ReadoutComponent, NoGradient, RFPulseComponent, GradientWaveform
 
 """
@@ -15,7 +15,7 @@ abstract type ContainerBlock <: AbstractBlock end
 
 Returns the start time of component with given `indices` with respect to the start of the [`ContainerBlock`](@ref).
 
-Also see [`duration`](@ref), [`end_time`](@ref), and [`effective_time`](@ref)
+Also see [`variables.duration`](@ref), [`end_time`](@ref), and [`effective_time`](@ref)
 """
 function start_time(container::ContainerBlock, index1, index2, indices...)
     start_time(container, index1) + start_time(container[index1], index2, indices...)
@@ -27,11 +27,11 @@ end
 
 Returns the start time of component with given `indices` with respect to the start of the [`ContainerBlock`](@ref).
 
-Also see [`duration`](@ref), [`start_time`](@ref), and [`effective_time`](@ref)
+Also see [`variables.duration`](@ref), [`start_time`](@ref), and [`effective_time`](@ref)
 """
 end_time(container::ContainerBlock, index, indices...) = start_time(container, index) + end_time(container[index], indices...)
-end_time(block::AbstractBlock) = duration(block)
-end_time(block::Tuple{<:VariableType, <:AbstractBlock}) = duration(block[2])
+end_time(block::AbstractBlock) = variables.duration(block)
+end_time(block::Tuple{<:VariableType, <:AbstractBlock}) = variables.duration(block[2])
 
 @defvar begin
     effective_time(container::ContainerBlock, index, indices...) = start_time(container, index) + variables.effective_time(container[index], indices...)
@@ -44,7 +44,7 @@ Returns the start time of component with given `indices` with respect to the sta
 
 This will crash if the component does not have an [`effective_time`](@ref) (e.g., if it is (part of) a gradient waveform).
 
-Also see [`duration`](@ref), [`start_time`](@ref), and [`end_time`](@ref)
+Also see [`variables.duration`](@ref), [`start_time`](@ref), and [`end_time`](@ref)
 """
 effective_time
 
diff --git a/src/containers/base_sequences.jl b/src/containers/base_sequences.jl
index 49f24b272c69b423ad2dd76505e87c535b6a7ef4..e121764ce2db7590db4d9810d2fc5fda46a68dc7 100644
--- a/src/containers/base_sequences.jl
+++ b/src/containers/base_sequences.jl
@@ -8,7 +8,7 @@ import ...Variables: get_free_variable, VariableType, variables, set_simple_cons
 import ...BuildSequences: global_model, global_scanner
 import ...Components: EventComponent, NoGradient, edge_times
 import ...Scanners: Scanner, B0
-import ..Abstract: ContainerBlock, start_time, readout_times
+import ..Abstract: ContainerBlock, start_time
 import ..BuildingBlocks: Wait, BuildingBlock, BaseBuildingBlock
 
 """
@@ -35,15 +35,15 @@ function Base.getindex(bs::BaseSequence{N}, index::Integer) where {N}
 end
 
 function (sequence::BaseSequence{N})(time::AbstractFloat) where {N}
-    var_time = mod(time, duration(sequence))
+    var_time = mod(time, variables.duration(sequence))
     for block in sequence
-        var_time -= duration(block)
+        var_time -= variables.duration(block)
         if var_time < 0
-            return (var_time + duration(block), block)
+            return (var_time + variables.duration(block), block)
         end
     end
     if abs(var_time) <= 1e-6
-        return (duration(sequence[N]) + var_time, sequence[N])
+        return (variables.duration(sequence[N]) + var_time, sequence[N])
     end
     error("Total duration of building blocks does not match duration of the sequence.")
 end
@@ -60,11 +60,11 @@ function start_time(bs::BaseSequence{N}, index::Integer) where {N}
         throw(BoundsError(bs, index))
     end
     base_index = ((index - 1) % N) + 1
-    base_time = sum(i -> duration(bs[i]), 1:base_index-1; init=0.)
+    base_time = sum(i -> variables.duration(bs[i]), 1:base_index-1; init=0.)
     if iszero(nTR)
         return base_time
     else
-        return nTR * duration(bs) + base_time
+        return nTR * variables.duration(bs) + base_time
     end
 end
 
@@ -101,8 +101,8 @@ How often sequence should be repeated.
 nrepeat(bs::BaseSequence) = 1
 
 @defvar begin
-    duration(bs::BaseSequence{0}) = 0.
-    duration(bs::BaseSequence) = sum(variables.duration.(bs); init=0.)
+    variables.duration(bs::BaseSequence{0}) = 0.
+    variables.duration(bs::BaseSequence) = sum(variables.duration.(bs); init=0.)
 end
 
 function edge_times(seq::BaseSequence; tol=1e-6)
@@ -181,7 +181,7 @@ Converst object into something that can be included in the sequence:
 to_block(cb::ContainerBlock) = cb
 to_block(s::Symbol) = to_block(Val(s))
 to_block(s::Union{VariableType, Nothing, Val{:min}, Val{:max}}) = Wait(s)
-to_block(ec::EventComponent) = BuildingBlock([NoGradient(duration(ec)), (0., ec)])
+to_block(ec::EventComponent) = BuildingBlock([NoGradient(variables.duration(ec)), (0., ec)])
 
 
 function make_generic(seq::BaseSequence)
diff --git a/src/containers/building_blocks.jl b/src/containers/building_blocks.jl
index ee27f14ba47af650a1977ad294968feb05935ddb..7d2db3d4ce2e3bff0af6de95ac7f7d2df9a3aabb 100644
--- a/src/containers/building_blocks.jl
+++ b/src/containers/building_blocks.jl
@@ -5,7 +5,7 @@ module BuildingBlocks
 import LinearAlgebra: norm
 import JuMP: @constraint
 import StaticArrays: SVector
-import ..Abstract: ContainerBlock, start_time, readout_times, end_time, iter
+import ..Abstract: ContainerBlock, start_time, end_time, iter
 import ...BuildSequences: global_model
 import ...Components: BaseComponent, GradientWaveform, EventComponent, NoGradient, ChangingGradient, ConstantGradient, split_gradient, RFPulseComponent, ReadoutComponent, InstantGradient, edge_times
 import ...Variables: VariableType, make_generic, get_pulse, get_readout, scanner_constraints!, get_gradient, gradient_orientation, variables, @defvar, get_free_variable
@@ -101,7 +101,7 @@ function waveform(bb::BaseBuildingBlock)
     end
     tol = sqrt(eps(Float64))
     for (_, block) in waveform_sequence(bb)
-        new_time = result[end][1] + max(duration(block), 0)
+        new_time = result[end][1] + max(variables.duration(block), 0)
         prev_grad = result[end][2]
         if block isa NoGradient
             @assert all(abs.(prev_grad) .<= 1e-12) "$(typeof(bb)) inserts NoGradient before the gradient is zero. This is probably caused by an improper implementation of this BuildingBlock."
@@ -132,7 +132,7 @@ function start_time(building_block::BaseBuildingBlock, index)
         component = building_block[key]
         if component isa GradientWaveform
             prev_time = time
-            time += duration(component)
+            time += variables.duration(component)
         end
         if equal_key(key, index)
             delay = component isa Tuple ? component[1] : 0.
@@ -168,14 +168,14 @@ function waveform_sequence(bb::BaseBuildingBlock, first, last)
         if equal_key(key, first)
             @assert !started
             started = true
-            current_start = effective_time(bb[key])
+            current_start = variables.effective_time(bb[key])
         end
         if equal_key(key, last)
             @assert started
             if isnothing(current_start)
-                push!(parts, (current_grad_key, split_gradient(bb[current_grad_key], effective_time(bb[key]))[1]))
+                push!(parts, (current_grad_key, split_gradient(bb[current_grad_key], variables.effective_time(bb[key]))[1]))
             else
-                push!(parts, (current_grad_key, split_gradient(bb[current_grad_key], current_start, effective_time(bb[key]))[2]))
+                push!(parts, (current_grad_key, split_gradient(bb[current_grad_key], current_start, variables.effective_time(bb[key]))[2]))
             end
             return parts
         end
@@ -198,7 +198,7 @@ end
         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)
+        t2 = isnothing(index2) ? variables.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.qvec(event)
@@ -271,7 +271,7 @@ function get_gradient(bb::BaseBuildingBlock, time::Number)
             return (block, time - start_time(bb, key))
         end
     end
-    error("$bb with duration $(duration(bb)) does not define a gradient at time $time.")
+    error("$bb with duration $(variables.duration(bb)) does not define a gradient at time $time.")
 end
 
 @defvar function gradient_strength(bb::BaseBuildingBlock, time::Number)
diff --git a/src/containers/containers.jl b/src/containers/containers.jl
index 0c0162949d0a86e76e3a5c66a54f99506253e43f..96198e5d98ad698ee789af1dacbcc23a85f1555e 100644
--- a/src/containers/containers.jl
+++ b/src/containers/containers.jl
@@ -4,7 +4,7 @@ include("building_blocks.jl")
 include("base_sequences.jl")
 include("alternatives.jl")
 
-import .Abstract: ContainerBlock, start_time, end_time, effective_time, readout_times, iter_blocks, iter_instant_gradients, iter_instant_pulses
+import .Abstract: ContainerBlock, start_time, end_time, iter_blocks, iter_instant_gradients, iter_instant_pulses
 import .BuildingBlocks: BaseBuildingBlock, BuildingBlock, Wait, waveform, waveform_sequence, events, ndim_grad
 import .BaseSequences: BaseSequence, Sequence, nrepeat, get_index_single_TR
 import .Alternatives: AlternativeBlocks, match_blocks!
diff --git a/src/parts/spoilt_slice_selects.jl b/src/parts/spoilt_slice_selects.jl
index ea4d3b1876f70de971d9ddf1f42aab5396da412a..2d2fe4de1bd55228f10d95fc60193f77511f7e66 100644
--- a/src/parts/spoilt_slice_selects.jl
+++ b/src/parts/spoilt_slice_selects.jl
@@ -90,14 +90,14 @@ end
 end
 
 Base.keys(::SpoiltSliceSelect) = Val.((:rise1, :flat1, :fall1, :flat_pulse, :pulse, :rise2, :flat2, :fall2))
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:rise1}) = ChangingGradient(0., slew_rate(spoilt), gradient_orientation(spoilt), rise_time(spoilt)[1], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat1}) = ConstantGradient(slew_rate(spoilt) * rise_time(spoilt)[1], gradient_orientation(spoilt), flat_time(spoilt)[1], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall1}) = ChangingGradient(slew_rate(spoilt) * rise_time(spoilt)[1], -slew_rate(spoilt), gradient_orientation(spoilt), fall_time(spoilt)[1], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat_pulse}) = ConstantGradient(slew_rate(spoilt) * spoilt.diff_time, gradient_orientation(spoilt), duration(spoilt.pulse), spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:rise1}) = ChangingGradient(0., variables.slew_rate(spoilt), gradient_orientation(spoilt), variables.rise_time(spoilt)[1], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat1}) = ConstantGradient(variables.slew_rate(spoilt) * variables.rise_time(spoilt)[1], gradient_orientation(spoilt), variables.flat_time(spoilt)[1], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall1}) = ChangingGradient(variables.slew_rate(spoilt) * variables.rise_time(spoilt)[1], -variables.slew_rate(spoilt), gradient_orientation(spoilt), variables.fall_time(spoilt)[1], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat_pulse}) = ConstantGradient(variables.slew_rate(spoilt) * spoilt.diff_time, gradient_orientation(spoilt), variables.duration(spoilt.pulse), spoilt.group)
 Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:pulse}) = spoilt.pulse
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:rise2}) = ChangingGradient(slew_rate(spoilt) * spoilt.diff_time, slew_rate(spoilt), gradient_orientation(spoilt), rise_time(spoilt)[2], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat2}) = ConstantGradient(slew_rate(spoilt) * fall_time(spoilt)[2], gradient_orientation(spoilt), flat_time(spoilt)[2], spoilt.group)
-Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall2}) = ChangingGradient(slew_rate(spoilt) * fall_time(spoilt)[2], -slew_rate(spoilt), gradient_orientation(spoilt), fall_time(spoilt)[2], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:rise2}) = ChangingGradient(variables.slew_rate(spoilt) * spoilt.diff_time, variables.slew_rate(spoilt), gradient_orientation(spoilt), variables.rise_time(spoilt)[2], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:flat2}) = ConstantGradient(variables.slew_rate(spoilt) * variables.fall_time(spoilt)[2], gradient_orientation(spoilt), variables.flat_time(spoilt)[2], spoilt.group)
+Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall2}) = ChangingGradient(variables.slew_rate(spoilt) * variables.fall_time(spoilt)[2], -variables.slew_rate(spoilt), gradient_orientation(spoilt), variables.fall_time(spoilt)[2], spoilt.group)
 
 @defvar begin
     rise_time(spoilt::SpoiltSliceSelect) = (spoilt.rise_time1, spoilt.fall_time2 - spoilt.diff_time)
@@ -107,17 +107,17 @@ Base.getindex(spoilt::SpoiltSliceSelect, ::Val{:fall2}) = ChangingGradient(slew_
 end
 
 @defvar begin
-    duration(spoilt::SpoiltSliceSelect) = sum(rise_time(spoilt)) + sum(flat_time(spoilt)) + sum(flat_time(spoilt)) + variables.duration(spoilt.pulse)
-    inverse_slice_thickness(spoilt::SpoiltSliceSelect) = slew_rate(spoilt) * spoilt.diff_time * variables.duration(spoilt.pulse) * 1e3
-    gradient_strength(spoilt::SpoiltSliceSelect) = slew_rate(spoilt) * max(spoilt.rise_time1, spoilt.fall_time2)
+    duration(spoilt::SpoiltSliceSelect) = sum(variables.rise_time(spoilt)) + sum(variables.flat_time(spoilt)) + sum(variables.flat_time(spoilt)) + variables.duration(spoilt.pulse)
+    inverse_slice_thickness(spoilt::SpoiltSliceSelect) = variables.slew_rate(spoilt) * spoilt.diff_time * variables.duration(spoilt.pulse) * 1e3
+    gradient_strength(spoilt::SpoiltSliceSelect) = variables.slew_rate(spoilt) * max(spoilt.rise_time1, spoilt.fall_time2)
 end
 
 get_pulse(spoilt::SpoiltSliceSelect) = spoilt.pulse
 
 @defvar function all_gradient_strengths(spoilt::SpoiltSliceSelect)
-    grad1 = spoilt.slew_rate * rise_time(spoilt)[1]
-    grad2 = grad1 - spoilt.slew_rate * flat_time(spoilt)[1]
-    grad3 = spoilt.slew_rate * fall_time(spoilt)[2]
+    grad1 = spoilt.slew_rate * variables.rise_time(spoilt)[1]
+    grad2 = grad1 - spoilt.slew_rate * variables.flat_time(spoilt)[1]
+    grad3 = spoilt.slew_rate * variables.fall_time(spoilt)[2]
     return [grad1, grad2, grad3]
 end
 
diff --git a/src/parts/trapezoids.jl b/src/parts/trapezoids.jl
index 0f7ab8e1e82d99962060d21c9b823e35c8de2be2..9dd2b5d0a3086720f8ee9549301b30ca3214110d 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(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))
+Base.getindex(pg::BaseTrapezoid{N}, ::Val{:rise}) where {N} = ChangingGradient(zeros(3), variables.slew_rate(pg), variables.rise_time(pg), get_group(pg))
+Base.getindex(pg::BaseTrapezoid, ::Val{:flat}) = ConstantGradient(variables.gradient_strength(pg), variables.flat_time(pg), get_group(pg))
+Base.getindex(pg::BaseTrapezoid, ::Val{:fall}) = ChangingGradient(variables.gradient_strength(pg), -variables.slew_rate(pg), variables.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
@@ -114,11 +114,11 @@ get_group(pg::BaseTrapezoid) = get_group(get_gradient(pg))
 end
 
 @defvar gradient begin
-    gradient_strength(g::Trapezoid) = slew_rate(g) .* rise_time(g)
-    δ(g::Trapezoid) = rise_time(g) + flat_time(g)
+    gradient_strength(g::Trapezoid) = variables.slew_rate(g) .* variables.rise_time(g)
+    δ(g::Trapezoid) = variables.rise_time(g) + variables.flat_time(g)
 end
 
-@defvar duration(g::BaseTrapezoid) = 2 * rise_time(g) + flat_time(g)
+@defvar duration(g::BaseTrapezoid) = 2 * variables.rise_time(g) + variables.flat_time(g)
 
 @defvar qvec(g::BaseTrapezoid, ::Nothing, ::Nothing) = variables.δ(g) .* variables.gradient_strength(g) .* 2π
 
@@ -209,18 +209,18 @@ struct LineReadout{N} <: BaseTrapezoid{N}
     ramp_overlap :: VariableType
 end
 
-function LineReadout(adc::ADC; ramp_overlap=nothing, orientation=nothing, group=nothing, variables...)
+function LineReadout(adc::ADC; ramp_overlap=nothing, orientation=nothing, group=nothing, vars...)
     res = LineReadout(
         Trapezoid(; orientation=orientation, group=group),
         adc,
         get_free_variable(ramp_overlap)
     )
-    set_simple_constraints!(res, variables)
+    set_simple_constraints!(res, vars)
     if !(res.ramp_overlap isa Number)
         @constraint global_model() res.ramp_overlap >= 0
         @constraint global_model() res.ramp_overlap <= 1
     end
-    @constraint global_model() (res.ramp_overlap * rise_time(res.trapezoid) * 2 + flat_time(res.trapezoid)) == duration(res.adc)
+    @constraint global_model() (res.ramp_overlap * variables.rise_time(res.trapezoid) * 2 + variables.flat_time(res.trapezoid)) == variables.duration(res.adc)
     return res
 end
 
diff --git a/src/pathways.jl b/src/pathways.jl
index 23b5c36cc8efc9383c637c20979258952f55ad77..3283817b1700cfc0c9709dc8b03b91c687edc8f7 100644
--- a/src/pathways.jl
+++ b/src/pathways.jl
@@ -86,10 +86,10 @@ The requested state can be set using `transverse` and `positive` as follows:
 - `transverse=false`, `positive=false`: -longitudinal
 - `transverse=true`, `positive=false`: -transverse
 """
-duration_state
+variables.duration_state
 
 @defvar pathway function duration_transverse(pathway::Pathway)
-    return duration_state(pathway, true, true) + duration_state(pathway, true, false)
+    return variables.duration_state(pathway, true, true) + variables.duration_state(pathway, true, false)
 end
 """
     duration_transverse(pathway::Pathway)
@@ -99,10 +99,10 @@ This determines the amount of T2-weighting as ``e^{t/T_2}``, where ``t`` is the
 
 Also see [`duration_dephase`](@ref) for T2'-weighting.
 """
-duration_transverse
+variables.duration_transverse
 
 @defvar pathway function duration_dephase(pathway::Pathway)
-    return duration_state(pathway, true, true) - duration_state(pathway, true, false)
+    return variables.duration_state(pathway, true, true) - variables.duration_state(pathway, true, false)
 end
 """
     duration_dephase(pathway::Pathway)
@@ -112,7 +112,7 @@ This determines the amount of T2'-weighting as ``e^{t/T_2'}``, where ``t`` is th
 
 Also see [`duration_transverse`](@ref) for T2-weighting.
 """
-duration_dephase
+variables.duration_dephase
 
 
 @defvar pathway net_dephasing(pathway::Pathway) = pathway.qvec
@@ -125,10 +125,10 @@ Only gradients active while the spins are in the transverse plane are considered
 
 Returns a NamedTuple with the `qvec` for all gradient groups.
 """
-net_dephasing
+variables.net_dephasing
 
 
-@defvar pathway area_under_curve(pathway::Pathway) = norm(qvec(pathway))
+@defvar pathway area_under_curve(pathway::Pathway) = norm(variables.net_dephasing(pathway))
 """
     area_under_curve(pathway::Pathway)
 
@@ -138,10 +138,10 @@ Only gradients active while the spins are in the transverse plane are considered
 
 Returns a NamedTuple with the `area_under_curve` for all gradient groups.
 """
-area_under_curve
+variables.area_under_curve
 
 
-@defvar pathway bmat(pathway::Pathway)  = pathway.bmat
+@defvar pathway bmat(pathway::Pathway) = pathway.bmat
 """
     bmat(pathway::Pathway)
 
@@ -151,9 +151,9 @@ Only gradients active while the spins are in the transverse plane are considered
 
 Returns a NamedTuple with the `bmat` for all gradient groups.
 """
-bmat
+variables.bmat
 
-@defvar pathway bval(pathway::Pathway) = tr(bmat(pathway))
+@defvar pathway bval(pathway::Pathway) = tr(variables.bmat(pathway))
 """
     bval(pathway::Pathway)
 
@@ -163,7 +163,7 @@ Only gradients active while the spins are in the transverse plane will contribut
 
 Returns a NamedTuple with the `bval` for all gradient groups.
 """
-bval
+variables.bval
 
 
 """
diff --git a/src/printing.jl b/src/printing.jl
index ec6b7f4a4c16503cd8e80341345d031128a281e9..3846185386be84d4486286e92776c53168504599 100644
--- a/src/printing.jl
+++ b/src/printing.jl
@@ -1,7 +1,7 @@
 module Printing
 import JuMP: value
 import Printf: @sprintf
-import ..Variables: VariableType, variables, AbstractBlock 
+import ..Variables: VariableType, variables, AbstractBlock, Variable
 
 function _robust_value(possible_number::VariableType)
     try
@@ -45,13 +45,18 @@ function Base.show(io::IO, block::AbstractBlock)
         print(io, name, "=", repr(getproperty(block, name)), ", ")
     end
 
-    for fn in values(variables)
+    for fn_name in names(variables, all=true)
+        fn = getproperty(variables, fn_name)
+        if !(fn isa Variable)
+            continue
+        end
+
         try
             numeric_value = _robust_value(fn(block))
             if isnothing(numeric_value)
                 continue
             end
-            print(io, "$(nameof(fn))=$(numeric_value), ")
+            print(io, "$(fn_name)=$(numeric_value), ")
         catch e
             continue
         end
diff --git a/src/sequence_io/pulseq.jl b/src/sequence_io/pulseq.jl
index 5f1c5f6d864d38dee32eab015e560675ffe471df..774a9d69c978b8cab8110359f47043edb5ce1d1f 100644
--- a/src/sequence_io/pulseq.jl
+++ b/src/sequence_io/pulseq.jl
@@ -172,7 +172,7 @@ function PulseqBlock(block::BaseBuildingBlock; BlockDurationRaster, AdcRasterTim
     end
     
     return PulseqBlock(
-        Int(div(1e-3 * duration(block), BlockDurationRaster, RoundNearest)),
+        Int(div(1e-3 * variables.duration(block), BlockDurationRaster, RoundNearest)),
         rf,
         grads...,
         adc,
diff --git a/src/variables.jl b/src/variables.jl
index 7df980611f87fe8818835f1fcffbea76f583d201..29779eedca5c47947933bbe7bdbb4c52bd15927c 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -179,7 +179,7 @@ function _defvar(func_def, getter=nothing)
             fn_def = MacroTools.splitdef(ex)
             push!(func_names, fn_def[:name])
             new_def = Dict{Symbol, Any}()
-            new_def[:name] = Expr(:., fn_def[:name], QuoteNode(:f))
+            new_def[:name] = Expr(:., Expr(:., :variables, QuoteNode(fn_def[:name])), QuoteNode(:f))
             new_def[:args] = esc.(fn_def[:args])
             new_def[:kwargs] = esc.(fn_def[:kwargs])
             new_def[:body] = esc(fn_def[:body])
@@ -206,24 +206,19 @@ function _defvar(func_def, getter=nothing)
     expressions = Expr[]
     for func_name in func_names
         push!(expressions, quote
-            $(esc(func_name)) = try
-                variables.$(func_name)
-            catch e
-                if !(e isa UndefVarError)
-                    rethrow()
-                end
+            if !($(QuoteNode(func_name)) in names(variables; all=true))
                 function $(func_name) end
                 variables.$(func_name) = Variable($(QuoteNode(func_name)), $(func_name), $getter)
             end
-            if $(esc(func_name)) isa AlternateVariable
-                error("$($(esc(func_name)).name) is defined through $($(esc(func_name)).other_var). Please define that variable instead.")
+            if variables.$(func_name) isa AlternateVariable
+                error("$($(esc(func_name)).name) is defined through $(variables.$(func_name).other_var). Please define that variable instead.")
             end
-            if !isnothing($getter) && $(esc(func_name)).getter != $getter
-                if isnothing($(esc(func_name)).getter)
-                    $(esc(func_name)).getter = $getter
+            if !isnothing($getter) && variables.$(func_name).getter != $getter
+                if isnothing(variables.$(func_name).getter)
+                    variables.$(func_name).getter = $getter
                 else
-                    name = $(esc(func_name)).name
-                    error("$(name) is already defined as a variable for $($(esc(func_name)).getter). Cannot switch to $($getter).")
+                    name = variables.$(func_name).name
+                    error("$(name) is already defined as a variable for $(variables.$(func_name).getter). Cannot switch to $($getter).")
                 end
             end
         end