diff --git a/docs/EmbeddPlotlyJSSyncPlotLiterate.jl b/docs/EmbeddPlotlyJSSyncPlotLiterate.jl
index 7e22f26a3..950ca0061 100644
--- a/docs/EmbeddPlotlyJSSyncPlotLiterate.jl
+++ b/docs/EmbeddPlotlyJSSyncPlotLiterate.jl
@@ -21,6 +21,9 @@ function Base.show(io::IO, ::MIME"text/html", fig::PlotlyJS.SyncPlot)
     # Copy is required because we pop layout width/height for Plotly's inner JSON.
     # Mutating fig.plot directly would change size(fig) and user-visible behavior.
     plot = copy(fig.plot)
+    plot.frames = fig.plot.frames
+    plot.config = deepcopy(fig.plot.config)
+    plot.config.displayModeBar = false
     default_width = layout_to_html_default_size!(plot.layout.fields, :width, get(plot.layout.fields, :width, nothing))
     default_height = layout_to_html_default_size!(plot.layout.fields, :height, get(plot.layout.fields, :height, nothing))
 
@@ -28,6 +31,7 @@ function Base.show(io::IO, ::MIME"text/html", fig::PlotlyJS.SyncPlot)
     PlotlyJS.PlotlyBase.to_html(
         html_buffer,
         plot;
+        autoplay=false,
         full_html=true,
         include_plotlyjs="cdn",
         default_width=default_width,
diff --git a/examples/3.tutorials/adiabatic_pulses/BIR4.mat b/examples/3.tutorials/adiabatic_pulses/BIR4.mat
deleted file mode 100644
index 83ff7fcdd..000000000
Binary files a/examples/3.tutorials/adiabatic_pulses/BIR4.mat and /dev/null differ
diff --git a/examples/3.tutorials/adiabatic_pulses/BIR4InversionB0B1.jl b/examples/3.tutorials/adiabatic_pulses/BIR4InversionB0B1.jl
deleted file mode 100644
index 4264503f6..000000000
--- a/examples/3.tutorials/adiabatic_pulses/BIR4InversionB0B1.jl
+++ /dev/null
@@ -1,53 +0,0 @@
-# This document replicates the results of Figure 15g of the paper "The Return of the Frequency Sweep:
-# Designing Adiabatic Pulses for Contemporary NMR" by Michael Garwood and Lance DelaBarre.
-using KomaMRI, MAT, PlotlyJS, LinearAlgebra
-
-RF_wf = matread("./examples/4.adiabatic_pulses/BIR4.mat")
-for R = [28] #Product duration and bandwidth
-    B1_nom = 13.5e-6; w1_2pi_Hz = γ*B1_nom
-    Trfs = [7] * 1e-3  # (4:0.1:7.5) * 1e-3 #ms
-    f0s = R ./ (2Trfs) # R = (2 fmax) * Trf
-    ΔB1s = Array(0.3:0.01:1.5) # %
-    NB0s = 200
-    B0_max = 100 #+-100Hz
-    Gz = B0_max / γ
-    z = range(-1, 1, NB0s)
-    B0s = Array(γ*Gz*z)
-
-    #Init
-    score = zeros(length(f0s), length(Trfs), length(ΔB1s))
-    MagXY = zeros(ComplexF64, NB0s, length(ΔB1s))
-    MagZ = zeros(ComplexF64, NB0s, length(ΔB1s))
-    N = prod(size(score))
-    counter = 1
-    for (i, f0) = enumerate(f0s), (j, Trf) = enumerate(Trfs), (k, ΔB1) = enumerate(ΔB1s)
-        b1 = B1_nom * ΔB1
-        B1 = b1 * RF_wf["b1"][:]
-        Δf = RF_wf["df"][:] * f0
-        dt = Trf / length(B1)
-        T90 = .5e-3
-        B190 = 90 / (360 * γ * T90)
-        rf90 = Sequence()
-        @addblock rf90 += (RF(ΔB1 * B190, T90, 0, 0), z=Grad(Gz, T90, 0))
-        bir4 = Sequence()
-        @addblock bir4 += (RF(B1, Trf, Δf, 0), z=Grad(Gz, Trf, 0))
-        seq = Sequence()
-        # @addblock seq += rf90
-        @addblock seq += bir4
-        # @addblock seq += rf90
-        sim_params = Dict{String,Any}("Δt_rf"=>dt)
-
-        M = simulate_slice_profile(seq; sim_params, z)
-        # display(plot_seq(seq))
-
-        MagXY[:, counter] = M.xy
-        MagZ[:, counter] = M.z
-        println("################ $counter / $N = $(round(counter/N*100; digits=3)) % ################")
-        counter = counter + 1
-    end
-    # Heatmap
-    fmax = round(f0s[1] * 1e-3; digits=2)
-    p1 = plot(heatmap(y=B0s, x=ΔB1s, z=abs.(MagXY), colorscale="Jet"), Layout(title="MXY BIR-4 Trf=7 ms R=$R fmax=$fmax kHz"))
-    p2 = plot(heatmap(y=B0s, x=ΔB1s, z=real.(MagZ), colorscale="Jet"), Layout(title="MZ BIR-4 Trf=7 ms R=$R fmax=$fmax kHz"))
-    display([p1; p2])
-end
diff --git a/examples/3.tutorials/adiabatic_pulses/BIR4InversionProfile.jl b/examples/3.tutorials/adiabatic_pulses/BIR4InversionProfile.jl
deleted file mode 100644
index 5d55f41f7..000000000
--- a/examples/3.tutorials/adiabatic_pulses/BIR4InversionProfile.jl
+++ /dev/null
@@ -1,37 +0,0 @@
-# This document replicates the results of Figure 15g of the paper "The Return of the Frequency Sweep:
-# Designing Adiabatic Pulses for Contemporary NMR" by Michael Garwood and Lance DelaBarre.
-using KomaMRI, MAT, PlotlyJS, LinearAlgebra, ProgressMeter
-
-RF_wf = matread("./examples/4.adiabatic_pulses/BIR4.mat")
-R = 200 #Product duration and bandwidth
-B1 = 120e-6; w1_2pi_Hz = γ*B1
-Trfs = [2] * 1e-3 #(4:0.1:7.5) * 1e-3 #ms
-f0s = R ./ (2Trfs) # R = 2 fmax * Trf
-ΔB1s = 0.8:0.2:1.2 # %
-#Init
-score = zeros(length(f0s), length(Trfs), length(ΔB1s))
-p = [scatter() for i=1:length(f0s), j=1:length(Trfs), k=1:length(ΔB1s), l=1:2]
-N = prod(size(score))
-count = 1
-for (i, f0) = enumerate(f0s), (j, Trf) = enumerate(Trfs), (k, ΔB1) = enumerate(ΔB1s)
-    b1 = 120e-6 * ΔB1
-    B1 = b1 * RF_wf["b1"][:]
-    Δf = RF_wf["df"][:] * f0
-    dt = Trf / length(B1)
-    fmax_sim = 10e3
-    Gz = fmax_sim / γ
-    seq = Sequence()
-    @addblock seq += (RF(B1, Trf, Δf, 0), z=Grad(Gz, Trf, 0))
-    sim_params = Dict{String,Any}("Δt_rf"=>dt)
-
-    z = range(-1, 1, 200)
-    M = simulate_slice_profile(seq; sim_params, z)
-
-    f = γ*Gz*z
-    p[i,j,k,1] = scatter(x=f,y=(1 .- M.z)/2,name="B1/B1_nom=$ΔB1")
-    # p[i,j,k,2] = scatter(x=f,y=abs.(M.xy),name="Mxy, B1/B1_nom=$ΔB1")
-    println("################ $count / $N = $(round(count/N*100; digits=3)) % ################")
-    global count = count + 1
-end
-# Heatmap
-plot(p[:], Layout(title="Fraction Refocused (BIR-4, Trf=$(Trfs[1]*1e3)ms, TBP=Trf*BW=$(R))"))
diff --git a/examples/3.tutorials/adiabatic_pulses/HSInversionProfile.jl b/examples/3.tutorials/adiabatic_pulses/HSInversionProfile.jl
deleted file mode 100644
index 6e0c6d262..000000000
--- a/examples/3.tutorials/adiabatic_pulses/HSInversionProfile.jl
+++ /dev/null
@@ -1,30 +0,0 @@
-using KomaMRI, PlotlyJS
-# RF Pulse Paramters, https://onlinelibrary.wiley.com/doi/epdf/10.1002/jmri.26021?saml_referrer
-b1max = 13e-6      #Peak amplitude (uT)
-Trf = 18.3e-3      #Pulse duration (ms)
-β = 4e2            #frequency modulation param (rad/s)
-μ = 6              #phase modulation parameter (dimensionless)
-fmax = μ * β / (2π) # 2fmax = BW
-# Adiabatic condition b1max >> β*√μ/γ:
-b1max > β*sqrt(μ)/(2π*γ)
-# Pulse Shape
-t = range(-Trf/2, Trf/2, 201)
-B1 =  b1max .* sech.(β.*t)
-Δf = -fmax  .* tanh.(β.*t)
-# Sequence generation
-fmax_sim = 2e3
-Gz = fmax_sim / γ
-seq = Sequence()
-@addblock seq += (RF(B1, Trf, Δf, 0), z=Grad(Gz, Trf, 0))
-p1 = plot_seq(seq; max_rf_samples=Inf, slider=false)
-KomaMRI.get_flip_angles(seq)[1]
-# Simulation
-sim_params = Dict{String,Any}("Δt_rf"=>t[2]-t[1])
-z = range(-1, 1, 400)
-M = simulate_slice_profile(seq; sim_params, z)
-# Plot
-f = γ*Gz*z
-s1 = scatter(x=f,y=abs.(M.xy),name="|Mxy|")
-s2 = scatter(x=f,y=M.z,name="Mz")
-p2 = plot([s1,s2], Layout(title="Hyperbolic-Secant (HS) Adiabatic Inversion Pulse (μ=$μ, β=$β rad/s)", xaxis_title="Frequency [Hz]"))
-[p1; p2]
diff --git a/examples/3.tutorials/lit-03-ChemicalShiftEPI.jl b/examples/3.tutorials/lit-03-ChemicalShiftEPI.jl
index 8fdeeac84..2faad835d 100644
--- a/examples/3.tutorials/lit-03-ChemicalShiftEPI.jl
+++ b/examples/3.tutorials/lit-03-ChemicalShiftEPI.jl
@@ -8,8 +8,22 @@ sys = Scanner(); #hide
 obj = brain_phantom2D() # a slice of a brain
 p1 = plot_phantom_map(obj, :T2 ; height=400, width=400, view_2d=true)
 p2 = plot_phantom_map(obj, :Δw ; height=400, width=400, view_2d=true)
-#md [p1 p2] #hide
-#jl display([p1 p2])
+p = [p1 p2] #hide
+p.plot.layout.fields[:xaxis2][:scaleanchor] = "y2" #hide
+p.plot.layout.fields[:xaxis1][:domain] = [0.0, 0.40] #hide
+p.plot.layout.fields[:xaxis2][:domain] = [0.58, 0.98] #hide
+for axis in (:xaxis1, :yaxis1, :xaxis2, :yaxis2) #hide
+    p.plot.layout.fields[axis][:title] = "" #hide
+    p.plot.layout.fields[axis][:ticks] = "" #hide
+    p.plot.layout.fields[axis][:showticklabels] = false #hide
+end #hide
+for (trace, x) in zip(p.plot.data, (0.43, 1.02)) #hide
+    trace.fields[:marker][:colorbar][:x] = x #hide
+    trace.fields[:marker][:colorbar][:len] = 0.55 #hide
+    trace.fields[:marker][:colorbar][:thickness] = 14 #hide
+end #hide
+#md p #hide
+#jl display(p)
 
 # At the left, you can see the ``T_2`` map of the phantom,
 # and at the right, the off-resonance ``\Delta\omega``.
diff --git a/examples/3.tutorials/lit-04-3DSliceSelective.jl b/examples/3.tutorials/lit-04-3DSliceSelective.jl
index 14ad80602..4968f7330 100644
--- a/examples/3.tutorials/lit-04-3DSliceSelective.jl
+++ b/examples/3.tutorials/lit-04-3DSliceSelective.jl
@@ -47,6 +47,7 @@ p3 = plot_signal(raw; slider=false, height=300)
 # Finally, we reconstruct the acquiered images.
 
 ## Get the acquisition data
+raw.params["trajectory"] = "other"
 acq = AcquisitionData(raw)
 
 ## Setting up the reconstruction parameters and perform reconstruction
@@ -58,5 +59,15 @@ image = reconstruction(acq, reconParams)
 p4 = plot_image(abs.(image[:, :, 1]); height=360, title="Slice 1")
 p5 = plot_image(abs.(image[:, :, 2]); height=360, title="Slice 2")
 p6 = plot_image(abs.(image[:, :, 3]); height=360, title="Slice 3")
-#md [p4 p5 p6] #hide
-#jl display([p4 p5 p6])
+p = [p4 p5 p6] #hide
+foreach(t -> t.fields[:showscale] = false, p.plot.data) #hide
+for (i, xref) in enumerate(("x", "x2", "x3")) #hide
+    xaxis = Symbol("xaxis", i) #hide
+    yaxis = Symbol("yaxis", i) #hide
+    p.plot.layout.fields[yaxis][:scaleanchor] = xref #hide
+    p.plot.layout.fields[yaxis][:constrain] = "domain" #hide
+    p.plot.layout.fields[xaxis][:range] = [-0.5, Nx - 0.5] #hide
+    p.plot.layout.fields[yaxis][:range] = [-0.5, Ny - 0.5] #hide
+end #hide
+#md p #hide
+#jl display(p)
diff --git a/examples/3.tutorials/lit-04b-AdiabaticRFPulse.jl b/examples/3.tutorials/lit-04b-AdiabaticRFPulse.jl
new file mode 100644
index 000000000..5c5e88cef
--- /dev/null
+++ b/examples/3.tutorials/lit-04b-AdiabaticRFPulse.jl
@@ -0,0 +1,255 @@
+# # Adiabatic RF Pulse
+
+using KomaMRI #hide
+sys = Scanner(); #hide
+
+# In this tutorial, we will build a hyperbolic-secant (HS) adiabatic inversion
+# pulse. The key point is that KomaMRI can keep the RF frequency modulation
+# explicit, then simulate the RF-frame dynamics and the resulting
+# ``B_0``/``B_1`` robustness.
+
+# ## Defining a frequency-modulated pulse
+#
+# > Parameters follow the default ``\mathrm{HS}_{4,6}`` pulse in Figure 1(c) of this
+# > [JMRI adiabatic inversion example](https://doi.org/10.1002/jmri.26021).
+
+b1max = 13.5e-6
+Trf = 18.3e-3
+β̂ = 4
+μ = 6
+β = 2 * β̂ / Trf;
+
+# For an HS pulse, the hardest point for adiabatic following is near the
+# center of the sweep, where ``B_1`` is maximal and
+# ``|\mathrm{d}\Delta\omega/\mathrm{d}t|=\mu\beta^2``. Requiring the RF precession rate to
+# dominate that sweep rate gives the threshold used in the paper:
+#
+# ```math
+# \frac{(\gamma_\mathrm{rad} B_1)^2}{\mu\beta^2} \ge 1
+# \quad\Rightarrow\quad
+# B_1 \geq \frac{\sqrt{\mu}\beta}{\gamma_\mathrm{rad}}.
+# ```
+
+b1_threshold = β * sqrt(μ) / (2π * γ)
+b1max > b1_threshold
+
+# First, we define the RF amplitude as a hyperbolic secant and the frequency
+# sweep as a hyperbolic tangent.
+
+t = range(-Trf / 2, Trf / 2, 201)
+B1 = b1max .* sech.(β .* t)
+Δf = -μ * β .* tanh.(β .* t) ./ (2π);
+
+# The ``\Delta f`` argument of `RF` can be a scalar, for a constant RF offset
+# such as a slice offset, or a waveform, as in this frequency-modulated pulse.
+
+f = range(-2e3, 2e3, 161) |> collect
+seq = Sequence()
+@addblock seq += RF(B1, Trf, Δf, 0);
+
+# This representation is native to KomaMRI. Pulseq RF events do not store
+# frequency-modulation waveforms, so exporting this sequence with `write_seq`
+# requires first converting the sweep into RF phase samples. KomaMRI may do this
+# conversion automatically in the future.
+
+# ## Plotting an adiabatic pulse
+#
+# `plot_seq` can show both views. The default plot keeps ``\Delta f(t)``
+# explicit; the `freq_in_phase` keyword shows the same pulse after moving the
+# frequency sweep into the RF phase.
+
+using PlotlyJS #hide
+function show_only_traces!(p, names) #hide
+    for trace in p.plot.data #hide
+        name = get(trace.fields, :name, "") #hide
+        trace.fields[:showlegend] = false #hide
+        trace.fields[:visible] = name in names #hide
+    end #hide
+    return p #hide
+end #hide
+function scale_trace_y!(p, name, scale) #hide
+    for trace in p.plot.data #hide
+        get(trace.fields, :name, "") == name || continue #hide
+        y = trace.fields[:y] #hide
+        trace.fields[:customdata] = y #hide
+        trace.fields[:y] = scale .* y #hide
+        trace.fields[:hovertemplate] = "(%{x:.4f} ms, Δf_FM: %{customdata:.4f} kHz)" #hide
+    end #hide
+    return p #hide
+end #hide
+p_freq = plot_seq(seq; max_rf_samples=Inf, slider=false, height=360, title="Frequency-modulated RF", showlegend=false)
+p_phase = plot_seq(seq; freq_in_phase=true, max_rf_samples=Inf, slider=false, height=360, title="Phase-modulated RF", showlegend=false)
+show_only_traces!(p_freq, ("|B1|_AM", "Δf_FM")) #hide
+show_only_traces!(p_phase, ("|B1|_AM", "∠B1_AM")) #hide
+scale_trace_y!(p_freq, "Δf_FM", 9) #hide
+p_rf = [p_freq p_phase] #hide
+relayout!(p_rf, height=380, margin=attr(t=56)) #hide
+p_rf #hide
+#jl display(p_rf)
+
+# Keeping the frequency sweep explicit is useful because KomaMRI simulates a
+# frequency-modulated RF pulse in the RF rotating frame. If the same pulse is
+# first converted to sampled phase modulation, the simulation is more sensitive
+# to the RF time sampling.
+
+# ## Comparing RF and rotating frames
+#
+# Using a callback, we can record the magnetization during the RF pulse.
+
+trajectory = NamedTuple[]
+call_every_N_blocks = 1
+
+record_traj = Callback(
+    call_every_N_blocks,
+    (progress_info, sim_blocks_info, device_data, sim_params) -> begin
+        j = last(sim_blocks_info.parts[progress_info.block])
+        push!(trajectory, (;
+            Mxy=device_data.Xt.xy[1], Mz=device_data.Xt.z[1],
+            ψ=device_data.seqd.ψ[j], B1=device_data.seqd.B1[j], Δf=device_data.seqd.Δf[j],
+        ))
+    end,
+)
+
+sim_params = KomaMRICore.default_sim_params()
+sim_params["return_type"] = "state"
+sim_params["max_rf_block_length"] = 1; # very inefficient; just for plots
+obj0 = Phantom(; x=[0.0], Δw=[0.0]);
+simulate(obj0, seq, sys; sim_params, callbacks=(record_traj,), verbose=false);
+
+function adiabatic_frame(p) #hide
+    ωeff = (-real(p.B1), -imag(p.B1), p.Δf / γ) #hide
+    ω̂rf = ωeff ./ sqrt(sum(abs2, ωeff)) #hide
+    Mxy_rf = p.Mxy * cis(-p.ψ) #hide
+    ω̂xy_rot = complex(ω̂rf[1], ω̂rf[2]) * cis(p.ψ) #hide
+    return (; #hide
+        Mrf=(real(Mxy_rf), imag(Mxy_rf), p.Mz), #hide
+        Mrot=(real(p.Mxy), imag(p.Mxy), p.Mz), #hide
+        ω̂rf, #hide
+        ω̂rot=(real(ω̂xy_rot), imag(ω̂xy_rot), ω̂rf[3]), #hide
+    ) #hide
+end #hide
+trajectory_frames = adiabatic_frame.(trajectory) #hide
+xyz(points) = ntuple(i -> getindex.(points, i), 3) #hide
+Mrf_xyz = xyz(getproperty.(trajectory_frames, :Mrf)) #hide
+Mrot_xyz = xyz(getproperty.(trajectory_frames, :Mrot)) #hide
+ω̂rf_xyz = xyz(getproperty.(trajectory_frames, :ω̂rf)) #hide
+ω̂rot_xyz = xyz(getproperty.(trajectory_frames, :ω̂rot)) #hide
+anim_idx = unique(round.(Int, range(1, length(trajectory_frames), 36))) #hide
+function sphere_mesh(; nθ=36, nφ=18) #hide
+    θ = range(0, 2π, nθ) #hide
+    φ = range(0, π, nφ) #hide
+    return ( #hide
+        x=[sin(p) * cos(t) for p in φ, t in θ], #hide
+        y=[sin(p) * sin(t) for p in φ, t in θ], #hide
+        z=[cos(p) for p in φ, t in θ], #hide
+    ) #hide
+end #hide
+sphere = sphere_mesh() #hide
+function bloch_traces(i, M, ω̂; scene="scene", showlegend=true) #hide
+    path_trace = scatter3d(; #hide
+        x=M[1][1:i], y=M[2][1:i], z=M[3][1:i], #hide
+        mode="lines", name="M path", scene=scene, showlegend=showlegend, #hide
+        line=attr(color="#111827", width=6), #hide
+    ) #hide
+    magnetization = scatter3d(; #hide
+        x=[0, M[1][i]], y=[0, M[2][i]], z=[0, M[3][i]], #hide
+        mode="lines+markers", name="M", scene=scene, showlegend=showlegend, #hide
+        line=attr(color="#dc2626", width=7), #hide
+        marker=attr(color="#dc2626", size=4), #hide
+    ) #hide
+    effective_field = scatter3d(; #hide
+        x=[0, ω̂[1][i]], y=[0, ω̂[2][i]], z=[0, ω̂[3][i]], #hide
+        mode="lines+markers", name="ωeff", scene=scene, showlegend=showlegend, #hide
+        line=attr(color="#2563eb", width=7), #hide
+        marker=attr(color="#2563eb", size=4), #hide
+    ) #hide
+    return [ #hide
+        path_trace, #hide
+        magnetization, #hide
+        effective_field, #hide
+    ] #hide
+end #hide
+bloch_sphere(scene) = [surface(; x=sphere.x, y=sphere.y, z=sphere.z, opacity=0.14, showscale=false, name="Bloch sphere", scene=scene, colorscale=[[0, "#dbe4ef"], [1, "#dbe4ef"]])] #hide
+function animation_frame(i) #hide
+    traces = [ #hide
+        bloch_traces(i, Mrf_xyz, ω̂rf_xyz)..., #hide
+        bloch_traces(i, Mrot_xyz, ω̂rot_xyz; scene="scene2", showlegend=false)..., #hide
+    ] #hide
+    return frame(; name=string(i), data=traces, traces=[1, 2, 3, 5, 6, 7]) #hide
+end #hide
+function bloch_scene(domain) #hide
+    axis(range) = attr(; #hide
+        range=range, title=attr(text=""), ticks="", #hide
+        showgrid=false, showbackground=false, showticklabels=false, zeroline=false, #hide
+    ) #hide
+    return attr(; #hide
+        domain=domain, aspectmode="cube", #hide
+        xaxis=axis([-1, 1]), yaxis=axis([-1, 1]), zaxis=axis([-1, 1]), #hide
+        camera=attr(eye=attr(x=1.31, y=1.31, z=0.75), up=attr(x=0, y=0, z=1)), #hide
+    ) #hide
+end #hide
+base_traces = [ #hide
+    bloch_sphere("scene")..., #hide
+    bloch_traces(first(anim_idx), Mrf_xyz, ω̂rf_xyz)..., #hide
+    bloch_sphere("scene2")..., #hide
+    bloch_traces(first(anim_idx), Mrot_xyz, ω̂rot_xyz; scene="scene2", showlegend=false)..., #hide
+] #hide
+frames = animation_frame.(anim_idx) #hide
+p_bloch = plot( #hide
+    base_traces, #hide
+    Layout(; #hide
+        height=340, margin=attr(l=0, r=0, t=32, b=0), #hide
+        scene=bloch_scene(attr(x=[0.0, 0.48], y=[0.0, 0.93])), #hide
+        scene2=bloch_scene(attr(x=[0.52, 1.0], y=[0.0, 0.93])), #hide
+        annotations=[ #hide
+            attr(text="RF frame", x=0.24, y=0.98, xref="paper", yref="paper", xanchor="center", showarrow=false, font=attr(size=18)), #hide
+            attr(text="Rotating frame", x=0.76, y=0.98, xref="paper", yref="paper", xanchor="center", showarrow=false, font=attr(size=18)), #hide
+        ], #hide
+        updatemenus=[ #hide
+            attr(; #hide
+                type="buttons", direction="right", x=0, y=1.0, #hide
+                xanchor="left", yanchor="top", #hide
+                buttons=[attr(label="Play", method="animate", args=[nothing, attr(frame=attr(duration=80, redraw=true), fromcurrent=true)])], #hide
+            ), #hide
+        ], #hide
+    ), #hide
+    frames, #hide
+) #hide
+#jl display(p_bloch)
+
+# The blue vector is the normalized rotation axis
+# ``\hat{\boldsymbol{\omega}}_\mathrm{eff}``, where
+# ``\boldsymbol{\omega}_\mathrm{eff} = -\gamma \mathbf{B}_\mathrm{eff}``. With this sign convention,
+# the magnetization precesses right-handed around the plotted axis.
+
+# ## HS adiabatic pulse B0 and B1 robustness
+#
+# To showcase the off-resonance and ``B_1`` robustness of this type of pulse, we
+# can show its effect in a heatmap with ``B_0 \in [-2, 2]\,\mathrm{kHz}`` and
+# ``B_1 \in [0, 16]\,\mu\mathrm{T}``:
+
+obj = Phantom(; x=zeros(length(f)), Δw=2π .* f);
+b1_scales = range(0.05, 1.2, 47) |> collect;
+sim_params = Dict{String, Any}("return_type" => "state");
+
+Mz = map(b1_scales) do scale
+    seq_scale = Sequence()
+    @addblock seq_scale += RF(scale .* B1, Trf, Δf, 0)
+    simulate(obj, seq_scale, sys; sim_params, verbose=false).z
+end;
+
+b1_axis = 1e6 .* b1max .* b1_scales #hide
+threshold = 1e6 * b1_threshold #hide
+Mz_map = round.(reduce(hcat, Mz); digits=3) #hide
+p_hs = plot( #hide
+    [ #hide
+        heatmap(; x=f, y=b1_axis, z=permutedims(Mz_map), zmin=-1, zmax=1, colorscale="RdBu", colorbar=attr(title="Mz")), #hide
+        scatter(; x=f, y=fill(threshold, length(f)), name="threshold", mode="lines", hoverinfo="skip", line=attr(color="#a62023", dash="dash", width=3)), #hide
+    ], #hide
+    Layout(; title="HS<sub>4,6</sub> inversion profile", xaxis_title="Off-resonance [Hz]", yaxis_title="B1,max [μT]", height=420, margin=attr(t=64)), #hide
+) #hide
+p_hs #hide
+#jl display(p_hs)
+
+# The dashed line is the analytic adiabatic threshold for this HS pulse, where we
+# can see that the inversion is achieved after the threshold is surpassed.
diff --git a/examples/3.tutorials/lit-06-DiffusionMotion.jl b/examples/3.tutorials/lit-06-DiffusionMotion.jl
index 0021e848b..0fa4c1f62 100644
--- a/examples/3.tutorials/lit-06-DiffusionMotion.jl
+++ b/examples/3.tutorials/lit-06-DiffusionMotion.jl
@@ -5,7 +5,7 @@ using PlotlyJS #hide
 using Random, Suppressor #hide
 
 # The purpose of this tutorial is to showcase the simulation of diffusion-related effects. 
-# For this, we are going to define a [`path`](@ref) motion to simulate the Brownian motion of spins.
+# For this, we are going to define a [`KomaMRI.path`](@ref) motion to simulate the Brownian motion of spins.
 # This is not the most efficient way of simulating diffusion, but it is a good way to understand the phenomenon.
 # In particular, we will going to simulate isotropic diffusion, characterized by the Apparent Diffusion Coefficient (ADC).
 
@@ -21,7 +21,7 @@ obj = Phantom(;
     T2 = ones(Nspins) * 100e-3,
 );
 
-# Now we will define the Brownian motion of the spins using the [`path`](@ref) motion definition.
+# Now we will define the Brownian motion of the spins using the [`KomaMRI.path`](@ref) motion definition.
 # The motion will be defined by the displacements in the x, y, and z directions (`dx`, `dy`, and `dz`) 
 # of the spins. The displacements will be generated by a random walk with mean square displacement 
 #