diff --git a/KomaMRICore/src/datatypes/Phantom.jl b/KomaMRICore/src/datatypes/Phantom.jl
index 0b945056e..6fc62d0ca 100644
--- a/KomaMRICore/src/datatypes/Phantom.jl
+++ b/KomaMRICore/src/datatypes/Phantom.jl
@@ -1,10 +1,10 @@
 """
-    phantom = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, ux, uy, uz)
+    obj = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, ux, uy, uz)
 
 The Phantom struct.
 
 # Arguments
-- `name`: (`::String`) name of the Phantom
+- `name`: (`::String`) name of the phantom
 - `x`: (`::AbstractVector{T}`, `[m]`) vector of x-positions of the spins
 - `y`: (`::AbstractVector{T}`, `[m]`) vector of y-positions of the spins
 - `z`: (`::AbstractVector{T}`, `[m]`) vector of z-positions of the spins
@@ -21,10 +21,17 @@ The Phantom struct.
 - `uz`: (`::Function`) displacement field in the z-axis
 
 # Returns
-- `phantom`: (`::Phantom`) Phantom struct
+- `obj`: (`::Phantom`) Phantom struct
+
+# Examples
+```julia-repl
+julia> obj = Phantom()
+
+julia> obj.ρ
+```
 """
  @with_kw mutable struct Phantom{T<:Real}
-	name::String = "spins"
+    name::String = "spins"
 	x::AbstractVector{T}
 	y::AbstractVector{T} =   zeros(size(x))
 	z::AbstractVector{T} =   zeros(size(x))
@@ -45,10 +52,26 @@ The Phantom struct.
 	uy::Function = (x,y,z,t)->0
 	uz::Function = (x,y,z,t)->0
 end
-# Phantom() = Phantom(name="spin",x=zeros(1,1))
+
+"""Size and length of a phantom"""
 size(x::Phantom) = size(x.ρ)
 Base.length(x::Phantom) = length(x.ρ)
 
+"""Compare two phantoms"""
+Base.isapprox(obj1::Phantom, obj2::Phantom)  = begin
+    obj1.x     ≈ obj2.x    &&
+    obj1.y     ≈ obj2.y    &&
+    obj1.z     ≈ obj2.z    &&
+    obj1.ρ     ≈ obj2.ρ    &&
+    obj1.T1    ≈ obj2.T1   &&
+    obj1.T2    ≈ obj2.T2   &&
+    obj1.T2s   ≈ obj2.T2s  &&
+    obj1.Δw    ≈ obj2.Δw   &&
+    obj1.Dλ1   ≈ obj2.Dλ1  &&
+    obj1.Dλ2   ≈ obj2.Dλ2  &&
+    obj1.Dθ    ≈ obj2.Dθ
+end
+
 """Separate object spins in a sub-group."""
 Base.getindex(obj::Phantom, p::AbstractRange) = begin
 	Phantom(name=obj.name,
@@ -70,7 +93,8 @@ Base.getindex(obj::Phantom, p::AbstractRange) = begin
 			uz=obj.uz
 			)
 end
-"""Separate object spins in a sub-group."""
+
+"""Separate object spins in a sub-group (lightweigth)."""
 Base.view(obj::Phantom, p::AbstractRange) = begin
 	@views Phantom(name=obj.name,
 			x=obj.x[p],
@@ -92,9 +116,8 @@ Base.view(obj::Phantom, p::AbstractRange) = begin
 			)
 end
 
-# Compartment enabling routines:
-# Addition of compartments
-+(s1::Phantom,s2::Phantom) =begin
+"""Addition of phantoms"""
++(s1::Phantom,s2::Phantom) = begin
 	Phantom(name=s1.name*"+"*s2.name,
 		x=[s1.x;s2.x],
 		y=[s1.y;s2.y],
@@ -114,7 +137,8 @@ end
 		uz=s1.uz
 	)
 end
-#Fraction of compartments
+
+"""Scalar multiplication of a phantom"""
 *(α::Real,obj::Phantom) = begin
 	Phantom(name=obj.name,
 		x=obj.x,
@@ -136,92 +160,6 @@ end
 	)
 end
 
-# Movement related commands
-# StartAt(s::Phantom,t0::Float64) = Phantom(s.name,s.x,s.y,s.ρ,s.T2,s.Δw,s.Dλ1,s.Dλ2,s.Dθ,(x,y,t)->s.ux(x,y,t.+t0),(x,y,t)->s.uy(x,y,t.+t0))
-# FreezeAt(s::Phantom,t0::Float64) = Phantom(s.name*"STILL",s.x.+s.ux(s.x,s.y,t0),s.y.+s.uy(s.x,s.y,t0),s.ρ,s.T2,s.Δw,s.Dλ1,s.Dλ2,s.Dθ,(x,y,t)->0,(x,y,t)->0)
-
-#TODO: jaw-pitch-roll, expand, contract, functions
-
-# Getting maps
-# get_DxDy2D(obj::Phantom) = begin
-# 	P(i) = rotz(obj.Dθ[i])[1:2,1:2];
-# 	D(i) = [obj.Dλ1[i] 0;0 obj.Dλ2[i]]
-# 	nx = [1;0]; ny = [0;1]
-# 	Dx = [nx'*P(i)'*D(i)*P(i)*nx for i=1:prod(size(obj.Dλ1))]
-# 	Dy = [ny'*P(i)'*D(i)*P(i)*ny for i=1:prod(size(obj.Dλ1))]
-# 	Dx, Dy
-# end
-
-"""
-    phantom = heart_phantom(α=1, β=1, γ=1, fat_bool::Bool=false)
-
-Heart-like LV phantom. The variable `α` is for streching, `β` for contraction, and `γ` for
-rotation.
-
-# Arguments
-- `α`: (`::Real`, `=1`) streching parameter
-- `β`: (`::Real`, `=1`) contraction parameter
-- `γ`: (`::Real`, `=1`) rotation parameter
-- `fat_bool`: (`::Bool`, `=false`) fat boolean parameter
-
-# Returns
-- `phantom`: (`::Phantom`) Heart-like LV phantom struct
-"""
-heart_phantom(α=1, β=1, γ=1, fat_bool::Bool=false) = begin
-	#PARAMETERS
-	FOV = 10e-2 #m Diameter ventricule
-	N = 21
-	Δxr = FOV/(N-1) #Aprox rec resolution, use Δx_pix and Δy_pix
-	Ns = 50 #number of spins per voxel
-	Δx = Δxr/sqrt(Ns) #spin separation
-	#POSITIONS
-	x = y = -FOV/2:Δx:FOV/2-Δx #spin coordinates
-	x, y = x .+ y'*0, x*0 .+ y' #grid points
-	#PHANTOM
-	⚪(R) =  (x.^2 .+ y.^2 .<= R^2)*1. #Circle of radius R
-	v = FOV/4 #m/s 1/16 th of the FOV during acquisition
-	ωHR = 2π/1 #One heart-beat in one second
-
-	# θ(t) = -π/4*γ*(sin.(ωHR*t).*(sin.(ωHR*t).>0)+0.25.*sin.(ωHR*t).*(sin.(ωHR*t).<0) )
-	ux(x,y,t) = begin
-		strech = 0 #α * v * (x.^2 .+ y.^2) / (FOV/2)^2 .* sign.(x)
-		contract = - β * v * x / (FOV/2)  #expand
-		rotate = - γ * v * y / (FOV/2)
-		def = (strech .+ contract .+ rotate) .* sin.(ωHR*t)
-	end
-	uy(x,y,t) = begin
-		strech = 0 #α * v * (x.^2 .+ y.^2) / (FOV/2)^2 .* sign.(y)
-		contract = - β * v * y / (FOV/2)
-		rotate = γ * v * x / (FOV/2)
-		def = (strech .+ contract .+ rotate) .* sin.(ωHR*t)
-	end
-	# Water spins
-	R = 9/10*FOV/2
-	r = 6/11*FOV/2
-	ring = ⚪(R) .- ⚪(r)
-	ρ = ⚪(r) .+ 0.9*ring #proton density
-	# Diffusion tensor model
-	D = 2e-9 #Diffusion of free water m2/s
-	D1, D2 = D, D/20
-	Dλ1 = D1*⚪(R) #Diffusion map
-	Dλ2 = D1*⚪(r) .+ D2*ring #Diffusion map
-	Dθ =  atan.(x,-y) .* ring #Diffusion map
-	T1 = (1400*⚪(r) .+ 1026*ring)*1e-3 #Myocardial T1
-	T2 = ( 308*⚪(r) .+ 42*ring  )*1e-3 #T2 map [s]
-	# Generating Phantoms
-	heart = Phantom("LeftVentricle",x,y,ρ,T2,Dλ1,Dλ2,Dθ,ux,uy)
-	# Fat spins
-	ring2 = ⚪(FOV/2) .- ⚪(R) #outside fat layer
-	ρ_fat = .5*ρ.*ring2
-	Δw_fat = 2π*220*ring2 #fat should be dependant on B0
-	T1_fat = 800*ring2*1e-3
-	T2_fat = 120*ring2*1e-3 #T2 map [s]
-	fat = Phantom("fat",x,y,ρ_fat,T2_fat,Δw_fat)
-	#Resulting phantom
-	obj = fat_bool ? heart + fat : heart #concatenating spins
-end
-
-
 """
     phantom = brain_phantom2D(;axis="axial", ss=4)
 
@@ -313,14 +251,14 @@ function brain_phantom2D(;axis="axial", ss=4)
 	T2 = T2*1e-3
 	T2s = T2s*1e-3
     phantom = Phantom{Float64}(name="brain2D_"*axis,
-					  x=y[ρ.!=0],
-					  y=x[ρ.!=0],
-					  z=0*x[ρ.!=0],
-					  ρ=ρ[ρ.!=0],
-					  T1=T1[ρ.!=0],
-					  T2=T2[ρ.!=0],
-					  T2s=T2s[ρ.!=0],
-					  Δw=Δw[ρ.!=0])
+				x=y[ρ.!=0],
+				y=x[ρ.!=0],
+				z=0*x[ρ.!=0],
+				ρ=ρ[ρ.!=0],
+				T1=T1[ρ.!=0],
+				T2=T2[ρ.!=0],
+				T2s=T2s[ρ.!=0],
+				Δw=Δw[ρ.!=0])
 	phantom
 end
 
@@ -418,13 +356,13 @@ function brain_phantom3D(;ss=4,start_end=[160, 200])
 	T2 = T2*1e-3
 	T2s = T2s*1e-3
     phantom = Phantom{Float64}(name ="brain3D",
-							   x = y[ρ.!=0],
-							   y = x[ρ.!=0],
-							   z = z[ρ.!=0],
-							   ρ = ρ[ρ.!=0],
-							   T1 = T1[ρ.!=0],
-							   T2 = T2[ρ.!=0],
-							   T2s = T2s[ρ.!=0],
-							   Δw = Δw[ρ.!=0])
+				x = y[ρ.!=0],
+				y = x[ρ.!=0],
+				z = z[ρ.!=0],
+				ρ = ρ[ρ.!=0],
+				T1 = T1[ρ.!=0],
+				T2 = T2[ρ.!=0],
+				T2s = T2s[ρ.!=0],
+				Δw = Δw[ρ.!=0])
 	phantom
 end
diff --git a/KomaMRICore/src/datatypes/Scanner.jl b/KomaMRICore/src/datatypes/Scanner.jl
index 634f53832..5f24e625e 100644
--- a/KomaMRICore/src/datatypes/Scanner.jl
+++ b/KomaMRICore/src/datatypes/Scanner.jl
@@ -1,14 +1,14 @@
 """
-    scanner = Scanner(B0, B1, Gmax, Smax, ADC_Δt, seq_Δt, GR_Δt, RF_Δt,
+    sys = Scanner(B0, B1, Gmax, Smax, ADC_Δt, seq_Δt, GR_Δt, RF_Δt,
         RF_ring_down_T, RF_dead_time_T, ADC_dead_time_T)
 
 The Scanner struct.
 
 # Arguments
-- `B0`: (`::Real`, `=1.5`, `[T]`) main magnetic field
+- `B0`: (`::Real`, `=1.5`, `[T]`) main magnetic field strength
 - `B1`: (`::Real`, `=10e-6`, `[T]`) maximum RF amplitude
-- `Gmax`: (`::Real`, `=60e-3`, `[T/m]`) maximum Gradient
-- `Smax`: (`::Real`, `=500`, `[mT/m/ms]`) maximum slew-rate
+- `Gmax`: (`::Real`, `=60e-3`, `[T/m]`) maximum gradient amplitude
+- `Smax`: (`::Real`, `=500`, `[mT/m/ms]`) gradient maximum slew-rate
 - `ADC_Δt`: (`::Real`, `=2e-6`, `[s]`) ADC raster time
 - `seq_Δt`: (`::Real`, `=1e-5`, `[s]`) sequence-block raster time
 - `GR_Δt`: (`::Real`, `=1e-5`, `[s]`) gradient raster time
@@ -18,7 +18,14 @@ The Scanner struct.
 - `ADC_dead_time_T`: (`::Real`, `=10e-6`, `[s]`) ADC dead time
 
 # Returns
-- `scanner`: (`::Scanner`) Scanner struct
+- `sys`: (`::Scanner`) Scanner struct
+
+# Examples
+```julia-repl
+julia> sys = Scanner()
+
+julia> sys.B0
+```
 """
 @with_kw mutable struct Scanner
     #Main
@@ -36,6 +43,3 @@ The Scanner struct.
     RF_dead_time_T::Real=100e-6
     ADC_dead_time_T::Real=10e-6
 end
-
-#Functions that check that Sequence satisfy hardware requirements:
-#check_sys_req(seq::Sequence,sys::Scanner)
diff --git a/KomaMRICore/src/datatypes/Sequence.jl b/KomaMRICore/src/datatypes/Sequence.jl
index eacecfe29..be3debcb8 100644
--- a/KomaMRICore/src/datatypes/Sequence.jl
+++ b/KomaMRICore/src/datatypes/Sequence.jl
@@ -436,204 +436,6 @@ Returns all the flip angles of the RF pulses in the sequence `x`.
 """
 get_flip_angles(x::Sequence) = get_flip_angle.(x.RF)[:]
 
-"""
-    y = get_ADC_on(seq::Sequence, t::Array{Float64,1})
-
-Get the ADC struct that are active.
-
-!!! note
-    This function is not being used.
-
-# Arguments
-- `seq`: (`::Sequence`) Sequence struct
-- `t`: (`::Vector{Float64}`, `[s]`) the time vector
-
-# Returns
-- `y`: (`::Vector{Bool}`) ADC struct that are active
-"""
-get_ADC_on(seq::Sequence, t::Array{Float64,1}) = begin
-	Δt = t[2]-t[1]
-	M, N = size(seq.ADC)
-	ADC = seq.ADC.N .> 0
-	T = [seq.ADC[i].T for i=1:N]; T = [0; T]
-	⊓(t) = (0 .<= t .< 1)*1.;
-	ADC_bool = sum([ADC[i]*⊓((t.-sum(T[1:i]))/(T[i+1]+Δt)) for i=1:N])
-	ADC_bool = (ADC_bool.>0)
-	circshift(convert.(Bool,	[:]),-1)
-end
-
-"""
-    bvalue = get_bvalue(DIF::Sequence)
-
-Calculates the `b`-value, in the PGSE sequnce is b = (2πγGδ)²(Δ-δ/3) and the normalized
-diffusion vector `n` from a *diffusion sequence*.
-
-    bvalue = get_bvalue(G, Δ, δ; null::Bool=false)
-
-Calculates the `b` value for PGSE and a moment nulled sequence.
-
-!!! note
-    This function is not being used.
-
-# Arguments
-- `DIF`: (`:Sequence`)
-- `G`
-- `Δ`
-- `δ`
-
-# Keywords
-- `null`: (`::Bool`, `=false`)
-
-# Returns
-- `bvalue`: the b value
-"""
-get_bvalue(DIF::Sequence) = begin
-	M, N = size(DIF.GR)
-	G = getproperty.(DIF.GR,:A) #[DIF.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
-	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
-	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
-	τ = dur(DIF) #End of sequence
-	# B-value
-	b = 0
-	for k=1:M,i=1:N,j=1:N
-		ij = max(i,j)
-		α = (i==j) ? 2/3 : 1/2
-		b += G[k,i]*G[k,j]*δ[i]*δ[j]*(τ-T[ij]-α*δ[ij])
-	end
-	b_value = (2π*γ)^2*b # Trace of B tensor
-end
-
-get_bvalue(G, Δ, δ; null::Bool=false) = (null ? 1/9 : 1)*(2π*γ*G*δ)^2*(Δ-δ/3)
-
-"""
-    bmatrix = get_Bmatrix(DIF::Sequence)
-
-Calculates the `b`-matrix, such as `b`-value = g' B g [s/mm2] with g [T/m].
-
-!!! note
-    This function is not being used.
-
-# Arguments
-- `DIF`: (`::Sequence`) diffusion Sequence struct
-
-# Returns
-- `bmatrix`: b matrix
-"""
-get_Bmatrix(DIF::Sequence) = begin
-	M, N = size(DIF.GR)
-	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
-	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
-	τ = dur(DIF) #End of sequence
-	# B-value, slower way of doing it
-	# b = zeros(M,N,N)
-	# for k=1,i=1:N,j=1:N
-	# 	ij = max(i,j)
-	# 	α = (i==j) ? 2/3 : 1/2
-	# 	b[k,i,j] = δ[i]*δ[j]*(τ-T[ij]-α*δ[ij])
-	# end
-	ij = [max(i,j) for i=1:N, j=1:N]
-	α = [(i==j) ? 2/3 : 1/2 for i=1:N, j=1:N]
-	b = (δ' .* δ) .* (τ .- T[ij] .- α .* δ[ij])
-	b_value = (2π*γ)^2*1e-6*b # Trace of B tensor
-	b_value
-end
-
-"""
-    qvector = get_qvector(DIF::Sequence; type::String="val")
-
-Calculates `q` = γδG diffusion vector from a *diffusion sequence*.
-
-!!! note
-    This function is not being used.
-
-# Arguments
-- `DIF`: (`::Sequence`) diffusion Sequence struct
-- `type`: (`::String`, `="val"`) option to get a value or plot
-
-# Returns
-- `qvector`: q diffusion vector as a value or a plot
-"""
-get_qvector(DIF::Sequence; type::String="val") = begin
-	M, N = size(DIF.GR)
-	G = getproperty.(DIF.GR,:A) #[DIF.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
-	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
-	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
-	τ = dur(DIF) #End of sequence
-	qτ = γ*sum([G[i,j]*δ[j] for i=1:M,j=1:N],dims=2) #Final q-value
-	norm(qτ) ≥ 1e-10 ? error("This is not a diffusion sequence, M0{G(t)}=∫G(t)dt!=0.") : 0
-	# q-trajectory
-	int_rect(t,A,δ) = (t.≥0) ? A*(t.-(t.≥δ).*(t.-δ)) : 0
-	q(t) = γ*sum([int_rect(t.-T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
- 	# Output q-vector and q-trajectory
-	if type=="val"
-		q(τ/2) #1/m
-	elseif type=="traj"
-		t = range(0,stop=τ,length=10^3)
-		qt = [q(t)[i] for i=1:M, t=t]
-		l = PlotlyJS.Layout(;title="q(t) trajectory", yaxis_title="q [1/mm]",
-	    xaxis_title="t [ms]",height=300)
-		p = [PlotlyJS.scatter() for j=1:M]
-		idx = ["qx" "qy" "qz"]
-		for j=1:M
-			p[j] = PlotlyJS.scatter(x=t*1e3, y=qt[j,:]*1e-3,name=idx[j])
-		end
-		PlotlyJS.plot(p, l)
-	end
-end
-
-"""
-    M0, M1, M2 = get_M0_M1_M2(seq::Sequence)
-
-Get the momentums `M0`, `M1` and `M2` from the sequence `seq`.
-
-!!! note
-    This function is not being used.
-
-# Arguments
-- `seq`: (`::Sequence`) Sequence struct
-
-# Returns
-- `M0`: M0 momentum
-- `M1`: M1 momentum
-- `M2`: M2 momentum
-"""
-get_M0_M1_M2(SEQ::Sequence) = begin
-	M, N = size(SEQ.GR)
-	G = [SEQ.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
-	δ = [SEQ.GR[1,j].T for j=1:N]; #Duration of pulse
-	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
-	τ = dur(SEQ) #End of sequence
-	#M0
-	M0i(t,T,A,δ) = (t.≥T) ? A*(t.-T.-(t.≥δ+T).*(t.-δ.-T)) : 0
-	M0(t) = sum([M0i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
-	#M1
-	M1i(t,T,A,δ) = (t.≥T) ? A/2*(t.^2 .-T^2 .-(t.≥δ+T).*(t.^2 .-(δ+T).^2))./2 : 0
-	M1(t) = sum([M1i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
-	#M1
-	M2i(t,T,A,δ) = (t.≥T) ? A/3*(t.^3 .-T^3 .-(t.≥δ+T).*(t.^3 .-(δ+T).^3))./3 : 0
-	M2(t) = sum([M2i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
-	M0, M1, M2
-end
-
-"""
-    Gmax = get_max_grad(x::Sequence)
-
-The maximum gradient of the sequence `x`.
-
-!!! note
-    This function is not being used.
-
-# Arguments
-- `seq`: (`::Sequence`) Sequence struct
-
-# Returns
-- `Gmax`: (`::Real`) maximum value of the gradient
-"""
-get_max_grad(x::Sequence) = begin
-	G = [x.GR[i,j].A for i=1:size(x.GR,1),j=1:size(x.GR,2)]
-	Gmax = maximum(sqrt.(sum(G.^2,dims=1)))
-end
-
 """
     rf_idx, rf_type = get_RF_types(seq, t)
 
@@ -907,5 +709,5 @@ get_eddy_currents(seq::Sequence; Δt=1, λ=80e-3) = begin
 	M2z_adc = linear_interpolation(ts,M2[:,3],extrapolation_bc=0)(t_adc)
 	M2_adc = [M2x_adc M2y_adc M2z_adc]
 	#Final
-	M2*1_000, M2_adc*1_000
+	M2, M2_adc
 end
diff --git a/KomaMRICore/src/datatypes/sequence/ADC.jl b/KomaMRICore/src/datatypes/sequence/ADC.jl
index d1709a6b9..85b148499 100644
--- a/KomaMRICore/src/datatypes/sequence/ADC.jl
+++ b/KomaMRICore/src/datatypes/sequence/ADC.jl
@@ -24,7 +24,7 @@ mutable struct ADC
     Δf::Float64
     ϕ::Float64
     function ADC(N, T, delay, Δf, ϕ)
-      T < 0 || delay < 0 ? error("ADC timings must be positive.") : new(N, T, delay, Δf, ϕ)
+        T < 0 || delay < 0 ? error("ADC timings must be positive.") : new(N, T, delay, Δf, ϕ)
     end
     function ADC(N, T, delay)
 		T < 0 || delay < 0 ? error("ADC timings must be positive.") : new(N, T, delay, 0, 0)
@@ -34,6 +34,12 @@ mutable struct ADC
     end
 end
 
+# ADC comparison
+Base.isapprox(adc1::ADC, adc2::ADC) = begin
+    return all(length(getfield(adc1, k)) ≈ length(getfield(adc2, k)) for k ∈ fieldnames(ADC))
+        all(getfield(adc1, k) ≈ getfield(adc2, k) for k ∈ fieldnames(ADC))
+end
+
 """
     y = getproperty(x::Vector{ADC}, f::Symbol)
 
@@ -43,7 +49,7 @@ directly without the need to iterate elementwise.
 # Arguments
 - `x`: (`::Vector{ADC}`) vector of ADC structs
 - `f`: (`::Symbol`, opts: [`:N`, `:T`, `:delay`, `:Δf`, `:ϕ`, `:dur`]) input symbol that
-    represents a property of the ACD structs
+    represents a property of the ADC structs
 
 # Returns
 - `y`: (`::Vector{Any}`) vector with the property defined by the `f` for all elements of
diff --git a/KomaMRICore/src/datatypes/sequence/Delay.jl b/KomaMRICore/src/datatypes/sequence/Delay.jl
index a6c2d25f6..4a21e982e 100644
--- a/KomaMRICore/src/datatypes/sequence/Delay.jl
+++ b/KomaMRICore/src/datatypes/sequence/Delay.jl
@@ -1,10 +1,8 @@
 """
     delay = Delay(T)
 
-The Delay struct. The input delay time `T` must be non-negative.
-
-!!! note
-    This struct is meant to add delays to a sequence struct.
+The Delay struct. It is a special "object" meant to add a delay to a sequence by using a
+sum operator.
 
 # Arguments
 - `T`: (`::Real`, `[s]`) time delay value
diff --git a/KomaMRICore/src/datatypes/sequence/Grad.jl b/KomaMRICore/src/datatypes/sequence/Grad.jl
index e20cf721c..e4e771df8 100644
--- a/KomaMRICore/src/datatypes/sequence/Grad.jl
+++ b/KomaMRICore/src/datatypes/sequence/Grad.jl
@@ -131,13 +131,13 @@ Base.show(io::IO, x::Grad) = begin
 	if !compact
 		wave = length(x.A) == 1 ? r(x.A*1e3) : "∿"
 		if x.rise == x.fall == 0.
-			print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"Grad($(wave) mT, $(r(x.T*1e3)) ms)")
+			print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"Grad($(wave) mT, $(r(sum(x.T)*1e3)) ms)")
 		else
-			print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"Grad($(wave) mT, $(r(x.T*1e3)) ms, ↑$(r(x.rise*1e3)) ms, ↓$(r(x.fall*1e3)) ms)")
+			print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"Grad($(wave) mT, $(r(sum(x.T)*1e3)) ms, ↑$(r(x.rise*1e3)) ms, ↓$(r(x.fall*1e3)) ms)")
 		end
 	else
 		wave = length(x.A) == 1 ? "⊓" : "∿"
-		print(io, (sum(abs.(x.A)) > 0 ? wave : "⇿")*"($(r((x.delay+x.rise+x.fall+x.T)*1e3)) ms)")
+		print(io, (sum(abs.(x.A)) > 0 ? wave : "⇿")*"($(r((x.delay+x.rise+x.fall+sum(x.T))*1e3)) ms)")
 	end
 end
 
@@ -177,6 +177,12 @@ getproperty(x::Matrix{Grad}, f::Symbol) = begin
 	end
 end
 
+# Gradient comparison
+Base.isapprox(gr1::Grad, gr2::Grad) = begin
+    return all(length(getfield(gr1, k)) ≈ length(getfield(gr2, k)) for k ∈ fieldnames(Grad)) &&
+        all(getfield(gr1, k) ≈ getfield(gr2, k) for k ∈ fieldnames(Grad))
+end
+
 # Gradient operations
 *(x::Grad,α::Real) = Grad(α*x.A,x.T,x.rise,x.fall,x.delay)
 *(α::Real,x::Grad) = Grad(α*x.A,x.T,x.rise,x.fall,x.delay)
@@ -184,7 +190,7 @@ end
 	N, M = size(GR)
 	[sum(A[i,1:N] .* GR[:,j]) for i=1:N, j=1:M]
 end
-zero(::Grad) = Grad(0,0)
+Base.zero(::Grad) = Grad(0,0)
 size(g::Grad, i::Int64) = 1 #To fix [g;g;;] concatenation of Julia 1.7.3
 /(x::Grad,α::Real) = Grad(x.A/α,x.T,x.rise,x.fall,x.delay)
 +(x::Grad,y::Grad) = Grad(x.A.+y.A,x.T,x.rise,x.fall,x.delay)
diff --git a/KomaMRICore/src/datatypes/sequence/RF.jl b/KomaMRICore/src/datatypes/sequence/RF.jl
index 394eecf9d..499676d08 100644
--- a/KomaMRICore/src/datatypes/sequence/RF.jl
+++ b/KomaMRICore/src/datatypes/sequence/RF.jl
@@ -88,6 +88,12 @@ getproperty(x::Matrix{RF}, f::Symbol) = begin
 	end
 end
 
+# RF comparison
+Base.isapprox(rf1::RF, rf2::RF) = begin
+    return all(length(getfield(rf1, k)) == length(getfield(rf2, k)) for k ∈ fieldnames(RF))
+        all(≈(getfield(rf1, k), getfield(rf2, k), atol=1e-9) for k ∈ fieldnames(RF))
+end
+
 # Properties
 size(r::RF, i::Int64) = 1 #To fix [r;r;;] concatenation of Julia 1.7.3
 *(α::Complex{T}, x::RF) where {T<:Real} = RF(α*x.A,x.T,x.Δf,x.delay)
diff --git a/KomaMRICore/src/io/MRiLab.jl b/KomaMRICore/src/io/MRiLab.jl
index f38e56010..2e64eec1c 100644
--- a/KomaMRICore/src/io/MRiLab.jl
+++ b/KomaMRICore/src/io/MRiLab.jl
@@ -48,9 +48,9 @@ function read_phantom_MRiLab(filename; B0=1.5, offset=[0,0,0], FRange_filename="
             # [((col+1)/2)*VOex.XDimRes;
             # ((row+1)/2)*VOex.YDimRes ;
             # ((layer+1)/2)*VOex.ZDimRes]; % Set matrix center as Object origin for motion simulation
-	xx = [(-FOVx/2:Δx[1]:FOVx/2)...;]
-	yy = [(-FOVy/2:Δx[2]:FOVy/2)...;;]
-	zz = [(-FOVz/2:Δx[3]:FOVz/2)...;;;]
+    xx = reshape((-FOVx/2:Δx[1]:FOVx/2),:,1,1) #[(-FOVx/2:Δx[1]:FOVx/2)...;]
+	yy = reshape((-FOVy/2:Δx[2]:FOVy/2),1,:,1) #[(-FOVy/2:Δx[2]:FOVy/2)...;;]
+	zz = reshape((-FOVz/2:Δx[3]:FOVz/2),1,1,:) #[(-FOVz/2:Δx[3]:FOVz/2)...;;;]
 	x = xx*1 .+ yy*0 .+ zz*0 .+ offset[1]	#spin x coordinates
 	y = xx*0 .+ yy*1 .+ zz*0 .+ offset[2]	#spin y coordinates
 	z = xx*0 .+ yy*0 .+ zz*1 .+ offset[3]	#spin z coordinates
diff --git a/KomaMRICore/src/sequences/PulseDesigner.jl b/KomaMRICore/src/sequences/PulseDesigner.jl
index 486428c92..14a12adbf 100644
--- a/KomaMRICore/src/sequences/PulseDesigner.jl
+++ b/KomaMRICore/src/sequences/PulseDesigner.jl
@@ -37,11 +37,9 @@ julia> plot_seq(ex)
 """
 RF_hard(B1, T, sys::Scanner; G=[0,0,0], Δf=0) = begin
 	ζ = sum(G) / sys.Smax
-	EX = Sequence([	Grad(G[1],T,ζ);	 #Gx
-					Grad(G[2],T,ζ);  #Gy
-					Grad(G[3],T,ζ);;], #Gz
-					 [RF(B1,T,Δf,ζ);;]	 #RF
-					)
+    gr = reshape([Grad(G[1],T,ζ); Grad(G[2],T,ζ); Grad(G[3],T,ζ)],:,1)
+    rf = reshape([RF(B1,T,Δf,ζ)],:,1)
+	EX = Sequence(gr, rf)
 	EX
 end
 
@@ -161,8 +159,8 @@ EPI(FOV::Float64, N::Int, sys::Scanner) = begin
 	# println("Pixel Δf in phase direction $(round(Δfx_pix_phase,digits=2)) Hz")
 	#Pre-wind and wind gradients
 	ϵ2 = Ta/(Ta+ζ)
-    PHASE =   Sequence(1/2*[Grad(-Ga, Ta, ζ)      ; ϵ2*Grad(-Ga, Ta, ζ);;]) #This needs to be calculated differently
-	DEPHASE = Sequence(1/2*[Grad((-1)^N*Ga, Ta, ζ); ϵ2*Grad(-Ga, Ta, ζ);;]) #for even N
+    PHASE =   Sequence(reshape(1/2*[Grad(      -Ga, Ta, ζ); ϵ2*Grad(-Ga, Ta, ζ)],:,1)) #This needs to be calculated differently
+	DEPHASE = Sequence(reshape(1/2*[Grad((-1)^N*Ga, Ta, ζ); ϵ2*Grad(-Ga, Ta, ζ)],:,1)) #for even N
 	seq = PHASE+EPI+DEPHASE
 	#Saving parameters
 	seq.DEF = Dict("Nx"=>Nx,"Ny"=>Ny,"Nz"=>1,"Name"=>"epi")
diff --git a/KomaMRICore/src/simulation/other/DiffusionModel.jl b/KomaMRICore/src/simulation/other/DiffusionModel.jl
index d71321904..e4fbb2b74 100644
--- a/KomaMRICore/src/simulation/other/DiffusionModel.jl
+++ b/KomaMRICore/src/simulation/other/DiffusionModel.jl
@@ -4,93 +4,7 @@ cross(n) = begin
    nx,ny,nz = n
     [0 -nz ny;
      nz 0 -nx;
-    -ny nx 0] 
+    -ny nx 0]
 end
 """Rodrigues' formula: Rotation matrix that when applied rotates with respect to "n" in an angle θ anti clock-wise"""
 Un(θ,n) = [1 0 0; 0 1 0; 0 0 1] * cos(θ) + sin(θ) * cross(n) + (1-cos(θ)) * (n * n')
-"""
-Slab oriented along the x axis.
-
-Bar, L., & Sochen, N. (2015). A spectral framework for NMR signal with restricted diffusion. Concepts in Magnetic Resonance Part A, 44(1), 16–53. doi:10.1002/cmr.a.21326 
-Grebenkov, D.S. (2008), Laplacian eigenfunctions in NMR. I. A numerical tool. Concepts Magn. Reson., 32A: 277-301. doi:10.1002/cmr.a.20117
-"""
-function Planes(L,D=2e-9,M=30)
-    #Normalization
-    ϵ(m) = m==0 ? 1 : sqrt(2)
-    #Eigen values
-    λ(m) = π^2*m^2
-    Λ =  D/L^2 * diagm([λ(m) for m=0:M])
-    #Matrix A
-    Ax = L * [m!=n ? ((-1)^(m+n)-1)*ϵ(m)*ϵ(n)*(λ(m)+λ(n))/(λ(m)-λ(n))^2 : 1/2 for m=0:M, n=0:M]
-    A = cat(Ax,0*Ax,0*Ax; dims=3) 
-    n = [1 0 0;
-         0 1 0;
-         0 0 1]
-    μ = (Λ, A, n)
-end
-"""
-Infinite cylinder oriented along the z axis.
-
-Bar, L., & Sochen, N. (2015). A spectral framework for NMR signal with restricted diffusion. Concepts in Magnetic Resonance Part A, 44(1), 16–53. doi:10.1002/cmr.a.21326 
-Grebenkov, D.S. (2008), Laplacian eigenfunctions in NMR. I. A numerical tool. Concepts Magn. Reson., 32A: 277-301. doi:10.1002/cmr.a.20117
-"""
-function Cylinder(R,D=2e-9,M=20)
-    #J'n(αnk) = 0
-    α = [1.841184, 3.054237, 3.831706, 4.201189, 5.317553,
-         5.331443, 6.415616, 6.706133, 7.015587, 7.501266,
-         8.015237, 8.536316, 8.577836, 9.282396, 9.647422, 
-         9.969468, 10.17347, 10.51986, 10.71143, 11.34592]
-    n = [1,2,0,3,4, 1,5,2,0,6, 3,1,7,4,8, 2,0,5,9,3]
-    #Eigen values
-    λ(nk) = nk==0 ? 0 : α[nk]^2
-    Λ =  D/R^2 * diagm([λ(nk) for nk=0:M])
-    #βnk
-    β(nk) = nk==0 ? 1 : sqrt(λ(nk)/(λ(nk)-n[nk]^2))    
-    #Some definitions
-    cond1(mk,nk) = abs(n[mk] - n[nk]) == 1
-    ϵ1(mk,nk) = sqrt( 1 + (n[mk] == 0) + (n[nk] == 0) )
-    cond2(mk,nk) = (n[mk] == n[nk]-1) - (n[mk] == n[nk]+1)
-    #Matrix A
-    Ax = R * [cond1(i,j)*ϵ1(i,j)*β(i)*β(j)*(λ(i)+λ(j)-2*n[i]*n[j])/(λ(i)-λ(j))^2 for i=0:M, j=0:M]
-    Ay = 1im*R * [cond2(i,j)*β(i)*β(j)*(λ(i)+λ(j)-2*n[i]*n[j])/(λ(i)-λ(j))^2 for i=0:M, j=0:M]
-    A = (Ax,Ay,0*Ax)
-    μ = (Λ, A)
-end
-
-"""
-Sphere of radius R.
-
-Bar, L., & Sochen, N. (2015). A spectral framework for NMR signal with restricted diffusion. Concepts in Magnetic Resonance Part A, 44(1), 16–53. doi:10.1002/cmr.a.21326 
-Grebenkov, D.S. (2008), Laplacian eigenfunctions in NMR. I. A numerical tool. Concepts Magn. Reson., 32A: 277-301. doi:10.1002/cmr.a.20117
-"""
-function Sphere(R,D=2e-9,M=20)
-    #j'n(αlk) = 0
-    α = [2.081576, 3.342094, 4.493409, 4.514100, 5.646704, 
-         5.940370, 6.756456, 7.289932, 7.725252, 7.851078, 
-         8.583755, 8.934839, 9.205840, 9.840446, 10.01037, 
-         10.61386, 10.90412, 11.07021, 11.07942, 11.97273]  
-    l = [1,2,0,3,4, 1,5,2,0,6, 3,7,1,4,8, 2,0,5,9,3]
-    m = [0,0,1,0,0, 1,0,1,2,0, 1,0,2,1,0, 2,3,1,0,2]
-    #Eigen values
-    λ(lm) = lm==0 ? 0 : α[lm]^2
-    Λ =  D/R^2 * diagm([λ(lm) for lm=0:M])
-    #βnk
-    β(lm) = lm==0 ? sqrt(3/2) : sqrt((2*l[lm]+1)*λ(lm)/(λ(lm)-l[lm]*(l[lm]+1)))    
-    # TODO: Finish matrix Ax
-    δlδm(i,j) = (abs(l[i]-l[j])==1)*(abs(m[i]-m[j])==1)
-    ϵ1(mk,nk) = (1+n[mk]+n[nk])/((2*n[mk]+1)*(2*n[nk]+1))
-    ϵ2(mk,nk) = n[mk]*(n[nk]+1) + n[nk]*(n[mk]+1) + 1
-    #Matrix A
-    Ax = R * [cond1(mk,nk)*ϵ1(mk,nk)*β(mk)*β(nk)*(λ(mk)+λ(nk)-ϵ2(mk,nk))/(λ(mk)-λ(nk))^2 for mk=0:M, nk=0:M]
-    A = (Ax,Ay,Az)
-    μ = (Λ, A)
-end
-
-function SignalE(μ, seq)
-    𝒊 = 1im;
-    M, N = size(seq.GR)
-    G, δ = seq.GR.A, seq.GR.T 
-    # E = [ Π exp( -(Λ + iγ Gn⋅A) ⋅ δn ) ]_{0,0}
-    A = [exp(-( μ[1] .+ 𝒊*2π*γ* .+([μ[2][:,:,m]*G[m,n] for m=1:M]...) )*δ[n]) for n=1:N]
-    E = *(A...)[1,1]
-end
diff --git a/KomaMRICore/src/simulation/other/OffResonance.jl b/KomaMRICore/src/simulation/other/OffResonance.jl
deleted file mode 100644
index 80487159c..000000000
--- a/KomaMRICore/src/simulation/other/OffResonance.jl
+++ /dev/null
@@ -1,50 +0,0 @@
-# B = μ (I-σI+1/3Χ) H
-#FT^-1{1/}
-
-using Plots
-using LinearAlgebra
-pyplot()
-
-Dipole(x,y) = begin
-    r = sqrt.(x.^2 .+ y.^2)
-    θ = atan.(x,y) #angle with respect to the main magnetic field
-    r < 1e-16 ? 0 : 1/(4*π)*(1 .+ 3*cos.(2*θ))./r.^3
-end
-
-B0 = 1 #From scanner object
-
-"""Schenck JF. The role of magnetic susceptibility in magnetic resonance
-imaging: MRI magnetic compatibility of the first and second kinds.
-Med Phys. 1996;23(6):815-850. doi:10.1118/1.597854"""
-Sphere(x,R,Δχ) = begin
-    if norm(x) <= R #Inside sphere
-        2Δχ/3*B0 #2/3
-    else #Outside
-        Δχ/3*B0*R^3*(2x[3]^2-x[1]^2-x[2]^2)/norm(x)^5
-    end
-end
-
-"""Schenck JF. The role of magnetic susceptibility in magnetic resonance
-imaging: MRI magnetic compatibility of the first and second kinds.
-Med Phys. 1996;23(6):815-850. doi:10.1118/1.597854"""
-CylinderZ(x,R,Δχ) = begin
-    if norm(x) <= R #Inside sphere
-        Δχ*B0 #1
-    else #Outside
-        0
-    end
-end
-
-"""Schenck JF. The role of magnetic susceptibility in magnetic resonance
-imaging: MRI magnetic compatibility of the first and second kinds.
-Med Phys. 1996;23(6):81 5-850. doi:10.1118/1.597854"""
-CylinderY(x,R,Δχ) = begin
-    if norm(x) <= R #Inside sphere
-        Δχ/2*B0 #1/2
-    else #Outside
-        Δχ/2*B0*R^2*(x[3]^2-x[1]^2)/(x[3]^2+x[1]^2)^2
-    end
-end
-
-# contourf([CylinderY([x,0,z],.1,.1) for z=-1:.001:1, x=-1:.001:1],aspect_ratio=:equal)
-# contourf([Sphere([x,0,z],.1,.1) for z=-1:.001:1, x=-1:.001:1],aspect_ratio=:equal)
diff --git a/KomaMRICore/src/simulation/other/TestCUDAKernel.jl b/KomaMRICore/src/simulation/other/TestCUDAKernel.jl
deleted file mode 100644
index a96cec5df..000000000
--- a/KomaMRICore/src/simulation/other/TestCUDAKernel.jl
+++ /dev/null
@@ -1,84 +0,0 @@
-using KomaMRI, CUDA
-
-#Test parameters
-ns = 50_000 #20630
-nt = 1_000
-#Phantom
-T1s = 1000.
-T2s = 50.
-M0s = 1.
-x = zeros(ns)
-T1 = T1s * ones(ns)
-T2 = T2s * ones(ns)
-M0 = M0s * ones(ns)
-obj = Phantom(x=x,ρ=M0,T1=T1,T2=T2)
-#DiscreteSequnce
-T = 50.
-Δt = T * ones(nt) / nt
-t = [0; cumsum(Δt)[1:end-1]]
-O = zeros(nt)
-A = zeros(Bool, nt)
-B1 = 1im*ones(nt)
-seqd = KomaMRI.DiscreteSequence(O,O,O,B1,O,A,t,Δt)
-#Magnetization
-Mxy = 1im*ones(ns)
-Mz = zeros(ns)
-M = Mag(Mxy, Mz)
-
-function BlochExcitationKernel(obj, M, seq)
-    φ = 0f0
-    nz = 1f0
-    nxy = 0f0im
-    for s ∈ seq
-        #Spinor Rotation
-        NVTX.@range "Rotation" begin
-        #Effective field
-        ΔB0 = obj.Δw./(2π*γ) .- s.Δf./γ
-        B = (s.Gx.*obj.x .+ s.Gy.*obj.y .+ s.Gz.*obj.z) .+ ΔB0
-        φ = -2π*γ*B.*s.Δt
-        #Spinor
-        α = cos.(φ./2) .- 1im.*nz.*sin.(φ./2)   
-        β = -1im.*nxy.*sin.(φ./2)
-        #Spinor rot
-        M.xy .= 2 .*conj.(α).*β.*M.z .+ conj.(α).^2 .*M.xy .- β.^2 .*conj.(M.xy)
-        M.z  .= (abs.(α).^2 .- abs.(β).^2).*M.z .- 2 .*real.(α.*β.*conj.(M.xy))
-        end
-        #Relaxation
-        NVTX.@range "Relaxation" begin
-        @. M.xy = M.xy * exp(-s.Δt / obj.T2)
-        @. M.z  = M.z  * exp(-s.Δt / obj.T1) + obj.ρ * (1 - exp(-s.Δt / obj.T1))
-        end
-    end
-    return M
-end
-
-function print_mxy()
-    Mxy_mean = sum(M.xy) / ns; Mz_mean = sum(M.z) / ns
-    println(exp(-sum(Δt)/T2s), " ", M0s*(1 - exp(-sum(Δt)/T1s)))
-    println(Mxy_mean, " ", Mz_mean)
-end
-println("##################################")
-println("------------ CPU Time ------------")
-M.xy = ones(ns); M.z = zeros(ns)
-M = BlochExcitationKernel(obj, M, seqd); print_mxy()
-M.xy = ones(ns); M.z = zeros(ns)
-stat_cpu = Base.@timed BlochExcitationKernel(obj, M, seqd);
-println("CPU time = $(stat_cpu.time) s")
-println("------------ GPU Time ------------")
-M.xy = ones(ns); M.z = zeros(ns)
-M = BlochExcitationKernel(KomaMRI.gpu(obj), KomaMRI.gpu(M), KomaMRI.gpu(seqd)); print_mxy()
-M.xy = ones(ns); M.z = zeros(ns)
-stat_gpu = CUDA.@timed BlochExcitationKernel(KomaMRI.gpu(obj), KomaMRI.gpu(M), KomaMRI.gpu(seqd));
-println("GPU time = $(stat_gpu.time) s")
-println("SPEEDUP = $(round(stat_cpu.time/stat_gpu.time))")
-# Serious profiling
-CUDA.reclaim()
-CUDA.@profile begin 
-    NVTX.@range "Objects to GPU" begin
-        obj_gpu  = obj  |> KomaMRI.gpu
-        M_gpu    = M    |> KomaMRI.gpu
-        seqd_gpu = seqd |> KomaMRI.gpu
-    end
-    NVTX.@range "BlochExcitationKernel" BlochExcitationKernel(obj_gpu, M_gpu, seqd_gpu)
-end
-CUDA.reclaim()
\ No newline at end of file
diff --git a/KomaMRICore/src/simulation/other/TestInterpolationGPU.jl b/KomaMRICore/src/simulation/other/TestInterpolationGPU.jl
deleted file mode 100644
index a54fe11a9..000000000
--- a/KomaMRICore/src/simulation/other/TestInterpolationGPU.jl
+++ /dev/null
@@ -1,15 +0,0 @@
-using Interpolations, Adapt, PlotlyJS
-
-b1 = [1, 1.5, 2, 1.5, 1, 0]
-dt = [1, 1, 2, 1, 1]
-t = [0; cumsum(dt)]
-itpp = extrapolate(interpolate((t, ), b1, Gridded(Constant{Previous}())), 0)
-itp = interpolate((t, ), b1, Gridded(Constant{Previous}()))
-cuitp = adapt(CuArray{Float32}, itp)
-x = 0:.001:sum(dt)
-y = @time itpp.(x);
-y = @time cuitp.(x);
-plot(x,y)
-##
-# using KomaMRI
-# (gx,gy,gz) = get_grads(seq, t)
\ No newline at end of file
diff --git a/KomaMRICore/test/runtests.jl b/KomaMRICore/test/runtests.jl
index daf9859b7..93b8068a9 100644
--- a/KomaMRICore/test/runtests.jl
+++ b/KomaMRICore/test/runtests.jl
@@ -3,6 +3,7 @@ using TestItems, TestItemRunner
 @run_package_tests filter=ti->!(:skipci in ti.tags)&&(:core in ti.tags) #verbose=true
 
 @testitem "Sequence" tags=[:core] begin
+    using Suppressor
     @testset "Init" begin
         sys = Scanner()
         B1 = sys.B1; durRF = π/2/(2π*γ*B1) #90-degree hard excitation pulse
@@ -72,29 +73,68 @@ using TestItems, TestItemRunner
         @test GR.A ≈ GR2.A
 
         #Sanity checks of contructors (A [T], T[s], rise[s], fall[s], delay[s])
+        A, T = 0.1, 1e-3
+        grad = Grad(A, T)
+
         A, T = rand(2)
         g1, g2 = Grad(A,T), Grad(A,T,0,0,0)
-        @test g1.A ≈ g2.A
-        @test g1.T ≈ g2.T
-        @test g1.rise ≈ g2.rise
-        @test g1.fall ≈ g2.fall
-        @test g1.delay ≈ g2.delay
+        @test g1 ≈ g2
 
         A, T, ζ = rand(3)
         g1, g2 = Grad(A,T,ζ), Grad(A,T,ζ,ζ,0)
-        @test g1.A ≈ g2.A
-        @test g1.T ≈ g2.T
-        @test g1.rise ≈ g2.rise
-        @test g1.fall ≈ g2.fall
-        @test g1.delay ≈ g2.delay
+        @test g1 ≈ g2
 
         A, T, delay, ζ = rand(4)
         g1, g2 = Grad(A,T,ζ,delay), Grad(A,T,ζ,ζ,delay)
-        @test g1.A ≈ g2.A
-        @test g1.T ≈ g2.T
-        @test g1.rise ≈ g2.rise
-        @test g1.fall ≈ g2.fall
-        @test g1.delay ≈ g2.delay
+        @test g1 ≈ g2
+
+        # Test construction with shape function
+        T, N = 1e-3, 100
+        f = t -> sin(π*t / T)
+        gradw = Grad(f, T, N)
+        @test gradw.A ≈ f.(range(0, T; length=N))
+
+        # Test Grad operations
+        α = 3
+        gradt = α * grad
+        @test size(grad, 1) == 1
+        @test gradt.A ≈ α * grad.A
+        gradt = grad * α
+        @test gradt.A ≈ α * grad.A
+        gradt = grad / α
+        @test gradt.A ≈ grad.A / α
+        grads = grad + gradt
+        @test grads.A ≈ grad.A + gradt.A
+        A1, A2, A3 = 0.1, 0.2, 0.3
+        v1 = [Grad(A1,T); Grad(A2,T); Grad(A3,T)]
+        v2 = [Grad(A2,T); Grad(A3,T); Grad(A1,T)]
+        v3 = v1 + v2
+        @test [v3[i].A for i=1:length(v3)] ≈ [v1[i].A + v2[i].A for i=1:length(v1)]
+        gradr = grad - gradt
+        @test gradr.A ≈ grad.A - gradt.A
+        gradt = -grad
+        @test gradt.A ≈ -grad.A
+        vc = vcat(v1, v2)
+        @test [vc[1,j].A for j=1:length(v1)] ≈ [v1[i].A for i=1:length(v1)]
+        @test [vc[2,j].A for j=1:length(v2)] ≈ [v2[i].A for i=1:length(v2)]
+        vc = vcat(v1, v2, v3)
+        @test [vc[1,j].A for j=1:length(v1)] ≈ [v1[i].A for i=1:length(v1)]
+        @test [vc[2,j].A for j=1:length(v2)] ≈ [v2[i].A for i=1:length(v2)]
+        @test [vc[3,j].A for j=1:length(v3)] ≈ [v3[i].A for i=1:length(v3)]
+        delay, rise, T, fall = 1e-6, 2e-6, 10e-3, 3e-6
+        gr = Grad(A, T, rise, fall, delay)
+        @test dur(gr) ≈ delay + rise + T + fall
+        T1, T2, T3 = 1e-3, 2e-3, 3e-3
+        vt = [Grad(A1,T1); Grad(A2,T2); Grad(A3,T3)]
+        @test dur(vt) ≈ [maximum([T1, T2, T3])]
+
+        # Test Grad output message
+        io = IOBuffer()
+        show(io, "text/plain", grad)
+        @test occursin("Grad(", String(take!(io)))
+        show(io, "text/plain", gr)
+        @test occursin("Grad(", String(take!(io)))
+
     end
 
     @testset "RF" begin
@@ -108,18 +148,177 @@ using TestItems, TestItemRunner
         #Sanity checks of constructors (A [T], T [s], Δf[Hz], delay [s])
         A, T = rand(2)
         r1, r2 = RF(A,T), RF(A,T,0,0)
-        @test r1.A ≈ r2.A
-        @test r1.T ≈ r2.T
-        @test r1.Δf ≈ r2.Δf
-        @test r1.delay ≈ r2.delay
+        @test r1 ≈ r2
 
         A, T, Δf = rand(3)
         r1, r2 = RF(A,T,Δf), RF(A,T,Δf,0)
-        @test r1.A ≈ r2.A
-        @test r1.T ≈ r2.T
-        @test r1.Δf ≈ r2.Δf
-        @test r1.delay ≈ r2.delay
+        @test r1 ≈ r2
+
+        # Test RF output message
+        io = IOBuffer()
+        show(io, "text/plain", r1)
+        @test occursin("RF(", String(take!(io)))
+
+        # Test Grad operations
+        B1x, B1y, T = rand(3)
+        A = B1x + im*B1y
+        α = Complex(rand())
+        rf = RF(A, T)
+        rft = α * rf
+        @test size(rf, 1) == 1
+        @test rft.A ≈ α * rf.A
+        @test dur(rf) ≈ rf.T
+        B1x, B1y, B2x, B2y, B3x, B3y, T1, T2, T3 = rand(9)
+        rf1, rf2, rf3 = RF(B1x + im*B1y, T1), RF(B1x + im*B1y, T2), RF(B3x + im*B3y, T3)
+        rv = [rf1; rf2; rf3]
+        @test dur(rv) ≈ sum(dur.(rv))
+
+    end
+
+    @testset "Delay" begin
+
+        # Test delay construction
+        T = 1e-3
+        delay = Delay(T)
+        @test delay.T ≈ T
+
+        # Test delay construction error for negative values
+        err = Nothing
+        try Delay(-T) catch err end
+        @test err isa ErrorException
+
+        # Test delay output message
+        io = IOBuffer()
+        show(io, "text/plain", delay)
+        @test occursin("Delay(", String(take!(io)))
+
+        # Test addition of a delay to a sequence
+        seq = Sequence()
+        ds = delay + seq
+        @test dur(ds[1]) ≈ delay.T && dur(ds[2]) ≈ .0
+        sd = seq + delay
+        @test dur(sd[1]) ≈ .0 && dur(sd[2]) ≈ delay.T
+
+    end
+    @testset "ADC" begin
+
+        # Test ADC construction
+        N, T, delay, Δf, ϕ  = 64, 1e-3, 2e-3, 1e-6, .25*π
+        adc = ADC(N, T, delay, Δf, ϕ)
+
+        adc1, adc2 = ADC(N, T), ADC(N,T,0,0,0)
+        @test adc1 ≈ adc2
+
+        adc1, adc2 = ADC(N, T, delay), ADC(N, T, delay, 0, 0)
+        @test adc1 ≈ adc2
+
+        adc1, adc2 = ADC(N, T, delay, Δf, ϕ), ADC(N, T, delay, Δf, ϕ)
+        @test adc1 ≈ adc2
+
+        # Test ADC construction errors for negative values
+        err = Nothing
+        try ADC(N, -T) catch err end
+        @test err isa ErrorException
+        try ADC(N, -T,  delay) catch err end
+        @test err isa ErrorException
+        try ADC(N,  T, -delay) catch err end
+        @test err isa ErrorException
+        try ADC(N, -T, -delay) catch err end
+        @test err isa ErrorException
+        try ADC(N, -T,  delay, Δf, ϕ) catch err end
+        @test err isa ErrorException
+        try ADC(N,  T, -delay, Δf, ϕ) catch err end
+        @test err isa ErrorException
+        try ADC(N, -T, -delay, Δf, ϕ) catch err end
+        @test err isa ErrorException
+
+        # Test ADC getproperties
+        Nb, Tb, delayb, Δfb, ϕb  = 128, 2e-3, 4e-3, 2e-6, .125*π
+        adb = ADC(Nb, Tb, delayb, Δfb, ϕb)
+        adcs = [adc, adb]
+        @test adcs.N ≈ [adc.N, adb.N] && adcs.T ≈ [adc.T, adb.T] && adcs.delay ≈ [adc.delay, adb.delay]
+        @test adcs.Δf ≈ [adc.Δf, adb.Δf] && adcs.ϕ ≈ [adc.ϕ, adb.ϕ] && adcs.dur ≈ [adc.T + adc.delay, adb.T + adb.delay]
+
+    end
+
+    @testset "DiscreteSequence" begin
+        path = @__DIR__
+        seq = @suppress read_seq(path*"/test_files/spiral.seq") #Pulseq v1.4.0, RF arbitrary
+        simParams = KomaMRICore.default_sim_params()
+        t, Δt = KomaMRICore.get_uniform_times(seq, simParams["Δt"]; Δt_rf=simParams["Δt_rf"])
+        seqd = KomaMRICore.discretize(seq)
+        i1, i2 = rand(1:Int(floor(0.5*length(seqd)))), rand(Int(ceil(0.5*length(seqd))):length(seqd))
+        @test seqd[i1].t ≈ [t[i1]]
+        @test seqd[i1:i2-1].t ≈ t[i1:i2]
+
+        T, N = 1, 4
+        seq = RF(1, 1)
+        seq += Sequence([Grad(1, 1)])
+        seq += ADC(N, 1)
+        simParams = KomaMRICore.default_sim_params()
+        simParams["Δt"], simParams["Δt_rf"] = T/N, T/N
+        seqd = KomaMRICore.discretize(seq; simParams)
+        i = Int(floor(length(seqd) / 3))
+        @test is_RF_on(seq[1]) == is_RF_on(seqd[1*i:1*i+1]) && is_GR_on(seq[1]) == is_GR_on(seqd[1*i:1*i+1]) && is_ADC_on(seq[1]) == is_ADC_on(seqd[1*i:1*i+1])
+        @test is_RF_on(seq[2]) == is_RF_on(seqd[2*i:2*i+1]) && is_GR_on(seq[2]) == is_GR_on(seqd[2*i:2*i+1]) && is_ADC_on(seq[2]) == is_ADC_on(seqd[2*i:2*i+1])
+        @test is_RF_on(seq[3]) == is_RF_on(seqd[3*i:3*i+1]) && is_GR_on(seq[3]) == is_GR_on(seqd[3*i:3*i+1]) && is_ADC_on(seq[3]) == is_ADC_on(seqd[3*i:3*i+1])
+        @test KomaMRICore.is_GR_off(seqd) ==  !KomaMRICore.is_GR_on(seqd)
+        @test KomaMRICore.is_RF_off(seqd) ==  !KomaMRICore.is_RF_on(seqd)
+        @test KomaMRICore.is_ADC_off(seqd) == !KomaMRICore.is_ADC_on(seqd)
+    end
+
+    @testset "SequenceFunctions" begin
+        path = @__DIR__
+        seq = @suppress read_seq(path*"/test_files/spiral.seq") #Pulseq v1.4.0, RF arbitrary
+        t, Δt = KomaMRICore.get_uniform_times(seq, 1)
+        t_adc =  KomaMRICore.get_adc_sampling_times(seq)
+        M2, M2_adc = KomaMRICore.get_slew_rate(seq)
+        Gx, Gy, Gz = KomaMRICore.get_grads(seq, t)
+        Gmx, Gmy, Gmz = KomaMRICore.get_grads(seq, reshape(t, 1, :))
+        @test reshape(Gmx, :, 1) ≈ Gx && reshape(Gmy, :, 1) ≈ Gy && reshape(Gmz, :, 1) ≈ Gz
+        @test is_ADC_on(seq) == is_ADC_on(seq, t)
+        @test is_RF_on(seq) == is_RF_on(seq, t)
+        @test KomaMRICore.is_Delay(seq) == !(is_GR_on(seq) || is_RF_on(seq) || is_ADC_on(seq))
+        @test size(M2, 1) == length(Δt) && size(M2_adc, 1) == length(t_adc)
+        io = IOBuffer()
+        show(io, "text/plain", seq)
+        @test occursin("Sequence[", String(take!(io)))
+
+        α = rand()
+        c = α + im*rand()
+        x = seq
+        y = @suppress read_seq(path*"/test_files/epi.seq") #Pulseq v1.4.0, RF arbitrary
+        z = x + y
+        @test z.GR.A ≈ [x.GR  y.GR].A && z.RF.A ≈ [x.RF y.RF].A && z.ADC.N ≈ [x.ADC; y.ADC].N
+        z = x - y
+        @test z.GR.A ≈ [x.GR  -y.GR].A
+        z = -x
+        @test z.GR.A ≈ -x.GR.A
+        z = x * α
+        @test z.GR.A ≈ α*x.GR.A
+        z = α * x
+        @test z.GR.A ≈ α*x.GR.A
+        z = x * c
+        @test z.RF.A ≈ c*x.RF.A
+        z = c * x
+        @test z.RF.A ≈ c*x.RF.A
+        z = x / α
+        @test z.GR.A ≈ x.GR.A/α
+        @test size(y) == size(y.GR[1,:])
+        z = x + x.GR[3,1]
+        @test z.GR.A[1, end] ≈ x.GR[3,1].A
+        z = x.GR[3,1] + x
+        @test z.GR.A[1, 1] ≈ x.GR[3,1].A
+        z = x + x.RF[1,1]
+        @test z.RF.A[1, end] ≈ x.RF[1,1].A
+        z = x.RF[1,1] + x
+        @test z.RF.A[1, 1] ≈ x.RF[1,1].A
+        z = x + x.ADC[3,1]
+        @test z.ADC.N[end] ≈ x.ADC[3,1].N
+        z = x.ADC[3,1] + x
+        @test z.ADC.N[1] ≈ x.ADC[3,1].N
     end
+
 end
 
 @testitem "PulseDesigner" tags=[:core] begin
@@ -154,18 +353,69 @@ end
 end
 
 @testitem "Phantom" tags=[:core] begin
+
+    # Test phantom struct creation
+    name = "Bulks"
+    x = [-2e-3; -1e-3; 0.; 1e-3; 2e-3]
+    y = [-4e-3; -2e-3; 0.; 2e-3; 4e-3]
+    z = [-6e-3; -3e-3; 0.; 3e-3; 6e-3]
+    ρ = [.2; .4; .6; .8; 1.]
+    T1 = [.9; .9; .5; .25; .4]
+    T2 = [.09; .05; .04; .07; .005]
+    T2s = [.1; .06; .05; .08; .015]
+    Δw = [-2e-6; -1e-6; 0.; 1e-6; 2e-6]
+    Dλ1 = [-4e-6; -2e-6; 0.; 2e-6; 4e-6]
+    Dλ2 = [-6e-6; -3e-6; 0.; 3e-6; 6e-6]
+    Dθ = [-8e-6; -4e-6; 0.; 4e-6; 8e-6]
+    u = (x,y,z,t)->0
+    obj = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ, ux=u, uy=u, uz=u)
+    obj2 = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ, ux=u, uy=u, uz=u)
+    @test obj ≈ obj2
+
+    # Test size and length definitions of a phantom
+    @test size(obj) == size(ρ)
+    #@test length(obj) == length(ρ)
+
+    # Test phantom comparison
+    ue = (x,y,z,t)->1
+    obe = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, ue, ue, ue)
+    @test obj ≈ obe
+
+    # Test phantom subset
+    rng = 1:2:5
+    ug = (x,y,z,t)->-1
+    obg = Phantom(name, x[rng], y[rng], z[rng], ρ[rng], T1[rng], T2[rng], T2s[rng], Δw[rng],
+                    Dλ1[rng], Dλ2[rng], Dθ[rng], ug, ug, ug)
+    @test obj[rng] ≈ obg
+    @test @view(obj[rng]) ≈ obg
+
+    # Test addition of phantoms
+    ua = (x,y,z,t)->2
+    oba = Phantom(name, [x; x[rng]], [y; y[rng]], [z; z[rng]], [ρ; ρ[rng]],
+                    [T1; T1[rng]], [T2; T2[rng]], [T2s; T2s[rng]], [Δw; Δw[rng]],
+                    [Dλ1; Dλ1[rng]], [Dλ2; Dλ2[rng]], [Dθ; Dθ[rng]], ua, ua, ua)
+    @test obj + obg ≈ oba
+
+    # Test scalar multiplication of a phantom
+    c = 7.
+    obc = Phantom(name, x, y, z, c*ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, u, u, u)
+    @test c*obj ≈ obc
+
     #Test brain phantom 2D
-    ph = brain_phantom2D()    #2D phantom
+    ph = brain_phantom2D()
     @test ph.name=="brain2D_axial"
 
     #Test brain phantom 3D
-    ph = brain_phantom3D()    #3D phantom
+    ph = brain_phantom3D()
     @test ph.name=="brain3D"
 end
 
 @testitem "Scanner" tags=[:core] begin
-    #Test somehow
-    @test true
+    B0, B1, Gmax, Smax = 1.5, 10e-6, 60e-3, 500
+    ADC_Δt, seq_Δt, GR_Δt, RF_Δt = 2e-6, 1e-5, 1e-5, 1e-6
+    RF_ring_down_T, RF_dead_time_T, ADC_dead_time_T = 20e-6, 100e-6, 10e-6
+    sys = Scanner(B0, B1, Gmax, Smax, ADC_Δt, seq_Δt, GR_Δt, RF_Δt, RF_ring_down_T, RF_dead_time_T, ADC_dead_time_T)
+    @test sys.B0 ≈ B0 && sys.B1 ≈ B1 && sys.Gmax ≈ Gmax && sys.Smax ≈ Smax
 end
 
 @testitem "Spinors×Mag" tags=[:core] begin
@@ -201,6 +451,43 @@ end
     xx2 = Un(θ,nx)*x; #3D rot matrix
     xx1 = [real(xx1.xy[1]), imag(xx1.xy[1]), xx1.z[1]]
     @test xx1 ≈ xx2
+
+    # Test Spinor struct
+    α, β = rand(2)
+    s = Spinor(α, β)
+    @test s[1].α ≈ [Complex(α)] && s[1].β ≈ [Complex(β)]
+    io = IOBuffer()
+    show(io, "text/plain", s)
+    @test occursin("Spinor(", String(take!(io)))
+    α2, β2 = rand(2)
+    s2 = Spinor(α2, β2)
+    sp = s * s2
+    @test sp.α ≈ s.α.*s2.α .- conj.(s2.β).*s.β && sp.β ≈ s.α.*s2.β .+ conj.(s2.α).*s.β
+    φ, φ1, θ, φ2 = rand(4)
+    Rm = KomaMRICore.Rg(φ1, θ, φ2)
+    @test Rm.α ≈ [cos(θ/2)*exp(-1im*(φ1+φ2)/2)] && Rm.β ≈ [sin(θ/2)*exp(-1im*(φ1-φ2)/2)]
+    Rn = KomaMRICore.Rφ(φ, θ)
+    @test Rn.α ≈ [cos(θ/2)+0im] && Rn.β ≈ [exp(1im*φ)*sin(θ/2)]
+    @test abs(s) ≈ [α^2 + β^2]
+end
+
+@testitem "TimeStepCalculation" tags=[:core] begin
+    ampRF = 1e-6
+    durRF = 1e-3
+    index_offset, number_rf_points = 2, 3
+    rf_key_points = []
+    rf = RF(ampRF, durRF)
+    seq = Sequence()
+    seq += rf
+    append!(rf_key_points, [index_offset; index_offset + number_rf_points + 1])
+    seq += Delay(durRF)
+    seq += rf
+    append!(rf_key_points, [rf_key_points[end] + 1; rf_key_points[end] + 1 + number_rf_points + 1])
+    seq += Delay(durRF)
+    seq += rf
+    append!(rf_key_points, [rf_key_points[end] + 1; rf_key_points[end] + 1 + number_rf_points + 1])
+    t, Δt = KomaMRICore.get_variable_times(seq; dt_rf=durRF)
+    @test KomaMRICore.get_breaks_in_RF_key_points(seq, t) == rf_key_points
 end
 
 @testitem "TrapezoidalIntegration" tags=[:core] begin
@@ -417,6 +704,25 @@ end
     @test  error0 + error1 + error2 < 0.1 #NMRSE < 0.1%
 end
 
+@testitem "BlochDict_CPU_single_thread" tags=[:important, :core] begin
+    using Suppressor
+    path = joinpath(@__DIR__, "test_files")
+    seq = @suppress read_seq(joinpath(path, "epi_100x100_TE100_FOV230.seq"))
+    obj = Phantom{Float64}(x=[0.], T1=[1000e-3], T2=[100e-3])
+    sys = Scanner()
+    simParams = Dict("gpu"=>false, "Nthreads"=>1, "sim_method"=>KomaMRICore.Bloch(), "return_type"=>"mat")
+    sig = @suppress simulate(obj, seq, sys; simParams)
+    sig = sig / prod(size(obj))
+    simParams["sim_method"] = KomaMRICore.BlochDict()
+    sig2 = simulate(obj, seq, sys; simParams)
+    sig2 = sig2 / prod(size(obj))
+    @test sig ≈ sig2
+
+    io = IOBuffer()
+    show(io, "text/plain", KomaMRICore.BlochDict())
+    @test occursin("BlochDict(", String(take!(io)))
+end
+
 @testitem "IO" tags=[:core] begin
     using Suppressor
     #Test Pulseq
@@ -448,6 +754,18 @@ end
         @test raw.params["patientName"] == "brain2D_axial"
         @test raw.params["trajectory"] == "other"
         @test raw.params["systemVendor"] == "KomaMRI.jl"
+
+        # Test signal_to_raw_data
+        signal1 = Vector()
+        for i=1:length(raw.profiles)
+            signal1 = [signal1; raw.profiles[i].data]
+        end
+        rawmrd = signal_to_raw_data(signal1, seq)
+        @test rawmrd.params["institutionName"] == raw.params["institutionName"]
+        io = IOBuffer()
+        show(io, "text/plain", rawmrd)
+        @test occursin("RawAcquisitionData[", String(take!(io)))
+
     end
     #Test JEMRIS
     @testset "JEMRIS" begin
@@ -455,4 +773,12 @@ end
         obj = read_phantom_jemris(path*"/test_files/column1d.h5")
         @test obj.name == "column1d.h5"
     end
-end
\ No newline at end of file
+    #Test JEMRIS
+    @testset "MRiLab" begin
+        path = @__DIR__
+        filename = path * "/test_files/brain_mrilab.mat"
+        FRange_filename = path * "/test_files/FRange.mat" #Slab within slice thickness
+        obj = read_phantom_MRiLab(filename; FRange_filename)
+        @test obj.name == "brain_mrilab.mat"
+    end
+end
diff --git a/KomaMRICore/test/test_files/FRange.mat b/KomaMRICore/test/test_files/FRange.mat
new file mode 100644
index 000000000..44e28cc2a
Binary files /dev/null and b/KomaMRICore/test/test_files/FRange.mat differ
diff --git a/KomaMRICore/test/test_files/brain_mrilab.mat b/KomaMRICore/test/test_files/brain_mrilab.mat
new file mode 100644
index 000000000..50bf44e1a
Binary files /dev/null and b/KomaMRICore/test/test_files/brain_mrilab.mat differ
diff --git a/KomaMRIPlots/Project.toml b/KomaMRIPlots/Project.toml
index 72c833b27..703c7748a 100644
--- a/KomaMRIPlots/Project.toml
+++ b/KomaMRIPlots/Project.toml
@@ -1,7 +1,7 @@
 name = "KomaMRIPlots"
 uuid = "76db0263-63f3-4d26-bb9a-5dba378db904"
 authors = ["Carlos Castillo Passi <cncastillo@uc.cl>"]
-version = "0.7.5"
+version = "0.7.6"
 
 [deps]
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
diff --git a/KomaMRIPlots/src/KomaMRIPlots.jl b/KomaMRIPlots/src/KomaMRIPlots.jl
index e7892864e..6ee62bc71 100644
--- a/KomaMRIPlots/src/KomaMRIPlots.jl
+++ b/KomaMRIPlots/src/KomaMRIPlots.jl
@@ -8,6 +8,7 @@ include("ui/DisplayFunctions.jl")
 using Reexport
 @reexport using PlotlyJS: savefig
 
-export plot_seq, plot_M0, plot_M1, plot_M2, plot_eddy_currents, plot_seqd, plot_kspace, plot_phantom_map, plot_signal, plot_image, plot_dict
+export plot_seq, plot_M0, plot_M1, plot_M2, plot_eddy_currents, plot_seqd,
+        plot_slew_rate, plot_kspace, plot_phantom_map, plot_signal, plot_image, plot_dict
 
 end
diff --git a/KomaMRIPlots/src/ui/DisplayFunctions.jl b/KomaMRIPlots/src/ui/DisplayFunctions.jl
index 0054c081a..574d2d1ac 100644
--- a/KomaMRIPlots/src/ui/DisplayFunctions.jl
+++ b/KomaMRIPlots/src/ui/DisplayFunctions.jl
@@ -373,6 +373,57 @@ julia> plot_eddy_currents(seq, 80e-3)
 ```
 """
 function plot_eddy_currents(seq, λ; α=ones(size(λ)), width=nothing, height=nothing, slider=true, show_seq_blocks=false, darkmode=false, range=[], title="")
+	#Times
+	dt = 1
+	t, Δt = KomaMRICore.get_uniform_times(seq + ADC(100, 100e-3), dt)
+	t = t[2:end]
+	ΔT = KomaMRICore.durs(seq)
+	T0 = cumsum([0; ΔT],dims=1)
+	Gx, Gy, Gz = KomaMRICore.get_grads(seq, t)
+	#Eddy currents per lambda
+	Gec = zeros(length(t), 3)
+	for (i, l) in enumerate(λ)
+		aux, _ =  KomaMRICore.get_eddy_currents(seq + ADC(100, 100e-3); Δt=dt, λ=l)
+		Gec .+= α[i] .* aux
+	end
+	#Plot eddy currents
+	p = [scatter() for j=1:4]
+	p[1] = scatter(x=t*1e3, y=(Gx*0 .+ Gec[:,1])*1e3, hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", name="ECx")
+	p[2] = scatter(x=t*1e3, y=(Gy*0 .+ Gec[:,2])*1e3, hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", name="ECy")
+	p[3] = scatter(x=t*1e3, y=(Gz*0 .+ Gec[:,3])*1e3, hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", name="ECz")
+	#Layout and config
+	l, config = generate_seq_time_layout_config(title, width, height, range, slider, show_seq_blocks, darkmode; T0)
+	plot(p, l; config)
+end
+
+"""
+    p = plot_slew_rate(seq; height=nothing, width=nothing, slider=true, darkmode=false, range=[])
+
+Plots the slew rate currents of a Sequence `seq`.
+
+# Arguments
+- `seq`: (`::Sequence`) Sequence
+
+# Keywords
+- `height`: (`::Int64`, `=nothing`) height of the plot
+- `width`: (`::Int64`, `=nothing`) width of the plot
+- `slider`: (`::Bool`, `=true`) boolean to display a slider
+- `darkmode`: (`::Bool`, `=false`) boolean to define colors for darkmode
+- `range`: (`::Vector{Float64}`, `=[]`) time range to be displayed initially
+
+# Returns
+- `p`: (`::PlotlyJS.SyncPlot`) plot of the slew rate currents of the sequence struct `seq`
+
+# Examples
+```julia-repl
+julia> seq_file = joinpath(dirname(pathof(KomaMRI)), "../examples/1.sequences/spiral.seq")
+
+julia> seq = read_seq(seq_file)
+
+julia> plot_slew_rate(seq)
+```
+"""
+function plot_slew_rate(seq; width=nothing, height=nothing, slider=true, show_seq_blocks=false, darkmode=false, range=[], title="")
 	#Times
 	dt = 1
 	t, Δt = KomaMRICore.get_uniform_times(seq, dt)
@@ -381,16 +432,12 @@ function plot_eddy_currents(seq, λ; α=ones(size(λ)), width=nothing, height=no
 	T0 = cumsum([0; ΔT],dims=1)
 	ts = t .+ Δt
 	#Eddy currents per lambda
-	k = zeros(length(t), 3)
-	for (i, l) in enumerate(λ)
-		aux, _ =  KomaMRICore.get_eddy_currents(seq; Δt=dt, λ=l)
-		k .+= α[i] .* aux
-	end
+	k, _ =  KomaMRICore.get_slew_rate(seq; Δt=dt)
 	#Plot eddy currents
 	p = [scatter() for j=1:4]
-	p[1] = scatter(x=ts*1e3, y=k[:,1], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2x")
-	p[2] = scatter(x=ts*1e3, y=k[:,2], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2y")
-	p[3] = scatter(x=ts*1e3, y=k[:,3], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2z")
+	p[1] = scatter(x=ts*1e3, y=k[:,1], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", name="SRx")
+	p[2] = scatter(x=ts*1e3, y=k[:,2], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", name="SRy")
+	p[3] = scatter(x=ts*1e3, y=k[:,3], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", name="SRz")
 	#Layout and config
 	l, config = generate_seq_time_layout_config(title, width, height, range, slider, show_seq_blocks, darkmode; T0)
 	plot(p, l; config)
diff --git a/KomaMRIPlots/test/runtests.jl b/KomaMRIPlots/test/runtests.jl
index 654794fd7..3eb865e17 100644
--- a/KomaMRIPlots/test/runtests.jl
+++ b/KomaMRIPlots/test/runtests.jl
@@ -95,6 +95,11 @@ using TestItems, TestItemRunner
             @test true              #If the previous line fails the test will fail
         end
 
+        @testset "plot_slew_rate" begin
+            plot_slew_rate(seq)
+            @test true
+        end
+
         @testset "plot_seqd" begin
             plot_seqd(seq)
             @test true
diff --git a/Project.toml b/Project.toml
index 79342e006..49d585684 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ version = "0.7.4"
 AssetRegistry = "bf4720bc-e11a-5d0c-854e-bdca1663c893"
 Blink = "ad839575-38b3-5650-b840-f874b8c74a25"
 Interact = "c601a237-2ae4-5e1e-952c-7a85b0c7eef1"
-InteractBase = "d3863d7c-f0c8-5437-a7b4-3ae773c01009"
 KomaMRICore = "4baa4f4d-2ae9-40db-8331-a7d1080e3f4e"
 KomaMRIPlots = "76db0263-63f3-4d26-bb9a-5dba378db904"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
@@ -19,9 +18,9 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 AssetRegistry = "0.1"
 Blink = "0.12"
 Interact = "0.10"
-InteractBase = "=0.10.5"
 KomaMRICore = "0.7"
 KomaMRIPlots = "0.7"
+MAT = "0.10"
 MRIReco = "0.6, 0.7"
 PlotlyJS = "0.18"
 Reexport = "1"
diff --git a/README.md b/README.md
index f1677d450..06a97fef9 100644
--- a/README.md
+++ b/README.md
@@ -11,11 +11,16 @@
 | [![][docr-img]][docr-url] | [![][paper-img1]][paper-url1]  | [![][gh-actions-img]][gh-actions-url] |
 | [![][docd-img]][docd-url] | [![][arXiv-img1]][arXiv-url1]  | [![][codecov-img]][codecov-url]       |
 ##### Submodules
+
 | **KomaMRI.jl**            | **KomaMRICore.jl**             | **KomaMRIPlots.jl**                   |
 |:-------------------------:|:------------------------------:|:-------------------------------------:|
 | [![][gh-actions-img1]][gh-actions-url] | [![][gh-actions-img2]][gh-actions-url] | [![][gh-actions-img3]][gh-actions-url]       |
 | [![][codecov-img1]][codecov-url] | [![][codecov-img2]][codecov-url] | [![][codecov-img3]][codecov-url]       |
-
+ ```
+📦 KomaMRI.jl (UI)                      
+├─ 📦 KomaMRICore.jl (Simulation and IO)
+└─ 📦 KomaMRIPlots.jl (Plots)           
+```
 </div>
 
 [docr-img]: https://img.shields.io/badge/docs-stable-blue.svg
@@ -46,6 +51,8 @@ KomaMRI.jl (formerly MRIsim.jl), whose name comes from the Japanese word for spi
 
 This package is meant to simulate general Magnetic Resonance Imaging (MRI) scenarios that could arise in pulse sequence development. 
 
+> 🟢 **[29 Jun 2023] [KomaMRI.jl's paper](https://onlinelibrary.wiley.com/doi/10.1002/mrm.29635) was part of July's editor's picks in MRM 🥳!**
+
 > 🟢 **[6 Mar 2023] Paper published in MRM 😃!** The open access article is available [here](https://onlinelibrary.wiley.com/doi/10.1002/mrm.29635).
 
 > 🟢 **[8 Dec 2022] KomaMRI v0.7**: Huge code rewrite, this implies improved: performance (now **5x faster**), type stability, extensibility, and more!
diff --git a/docs/src/docstrings.md b/docs/src/docstrings.md
index 524ac5530..13da70e59 100644
--- a/docs/src/docstrings.md
+++ b/docs/src/docstrings.md
@@ -25,10 +25,6 @@ Refer to [API Documentation](api.md#Phantom):
 * [`KomaMRI.brain_phantom2D`](@ref) — Function
 * [`KomaMRI.brain_phantom3D`](@ref) — Function
 
-```@docs
-KomaMRI.heart_phantom
-```
-
 ### `Scanner`
 
 Refer to [API Documentation](api.md#Scanner):
diff --git a/examples/5.fat_sat_low_field/CardiacLowFieldbSSFPSignal.jl b/examples/5.fat_sat_low_field/CardiacLowFieldbSSFPSignal.jl
index 6ef9144a5..c9686245b 100644
--- a/examples/5.fat_sat_low_field/CardiacLowFieldbSSFPSignal.jl
+++ b/examples/5.fat_sat_low_field/CardiacLowFieldbSSFPSignal.jl
@@ -11,16 +11,16 @@ fat_freq = γ*B0*fat_ppm
 off = Array(range(-1, 1, Niso))
 dx = Array(range(-Δx_voxel/2, Δx_voxel/2, Niso))
 #BASED ON T1 and T2 maps
-myocard = Phantom{Float64}(x=dx, T1=700e-3*ones(Niso), T2=60e-3*ones(Niso))
-blood =   Phantom{Float64}(x=dx, T1=1250e-3*ones(Niso), T2=300e-3*ones(Niso)) #1.5 T2 240ms
-fat =     Phantom{Float64}(x=dx, T1=172.93e-3*ones(Niso), T2=135.32e-3*ones(Niso), Δw=2π*(fat_freq .+ 10*off))
+myocard = Phantom{Float64}(x=dx, T1=700e-3*ones(Niso), T2=60e-3*ones(Niso), Δw=2π*20*off)
+blood =   Phantom{Float64}(x=dx, T1=1250e-3*ones(Niso), T2=300e-3*ones(Niso), Δw=2π*20*off) #1.5 T2 240ms
+fat =     Phantom{Float64}(x=dx, T1=172.93e-3*ones(Niso), T2=135.32e-3*ones(Niso), Δw=2π*(fat_freq .+ 20*off))
 obj = myocard+blood+fat
 #Sequence parameters
-# TR = 5.29e-3
-TR = 4.89e-3
+TR = 5.29e-3
+# TR = 4.89e-3
 im_flip_angle = 110
 iNAV_lines = 14
-FatSat_flip_angle = 180
+FatSat_flip_angle = 110
 #Pre-pulses
 T2prep_duration = 50e-3
 Tfatsat = 32e-3
@@ -32,7 +32,7 @@ Tadc = 1e-6
 RR = 1 #1 [s]
 Trf = 1.4e-3  #1 [ms]
 B1 = 1 / (360*γ*Trf)
-im_segments = 22
+im_segments = 20
 iNAV_flip_angle = 3.2
 number_dummy_heart_beats = 3
 #Scanner
@@ -100,21 +100,19 @@ function CMRA_iNAV_bSSFP_cardiac(number_dummy_heart_beats, iNAV_lines, im_segmen
         bssfp = bSSFP(iNAV_lines, im_segments, iNAV_flip_angle, im_flip_angle; sample=sample_recovery[hb+1])
         RRdelay = RR  - dur(bssfp) - dur(t2p) - dur(fatsat)
         seq += t2p
-        seq += Delay(4e-3)
         seq += fatsat
-        seq += Delay(4e-3)         
         seq += bssfp
-        seq += sample_recovery[hb+1] ? ADC(40, RRdelay) : Delay(RRdelay) #Sampling recovery curve
+        seq += sample_recovery[hb+1] ? ADC(80, RRdelay) : Delay(RRdelay) #Sampling recovery curve
     end
     return seq
 end
 
-sim_method = BlochDict(save_Mz=true)
 #Sequence
 seq = CMRA_iNAV_bSSFP_cardiac(number_dummy_heart_beats, iNAV_lines, 
                         im_segments, iNAV_flip_angle, im_flip_angle; sample_recovery=[false, false, false, true])
 #Simulation
-simParams = Dict{String,Any}("return_type"=>"mat", "sim_method"=>sim_method, "Nthreads"=>1)
+sim_method = BlochDict(save_Mz=true)
+simParams = Dict{String,Any}("return_type"=>"mat", "sim_method"=>sim_method, "gpu"=>false, "Nthreads"=>1)
 magnetization = simulate(obj, seq, sys; simParams)
 
 labels = ["Myocardium", "Blood", "Fat"]
@@ -134,9 +132,16 @@ end
 seqd = KomaMRICore.discretize(seq; simParams)
 p3 = scatter(x=seqd.t, y=abs.(seqd.B1)*1e6, name="B1",marker_color="purple",yaxis_range=[0,5])
 add_trace!(p0, p3, row=3, col=1)
-relayout!(p0, xaxis_range=[3, 4],title_text="TR=$(TR*1e3), α=$(im_flip_angle), iNAV_lines=$(iNAV_lines), FatSat α=$(FatSat_flip_angle)")
+add_layout_image!(p0, attr(
+                    source="https://raw.githubusercontent.com/cncastillo/KomaMRI.jl/master/src/ui/assets/Logo.svg",
+                    xref="x domain",
+                    yref="y domain",
+                    x=0.99,
+                    y=1.2,
+                    opacity=0.7,
+                    xanchor="right",
+                    yanchor="top",
+                    sizex=0.15,
+                    sizey=0.15,))
+relayout!(p0, yaxis_range=[0, 0.5], xaxis_range=[3, 4],title_text="TR=$(round(TR*1e3;digits=3)) ms, α=$(im_flip_angle), iNAV_lines=$(iNAV_lines), FatSat α=$(FatSat_flip_angle)")
 p0
-
-
-# plot_seq(seq)
-# p = KomaMRIPlots.plot_seqd(seq; simParams); p
\ No newline at end of file
diff --git a/examples/DiffPrepMCMRF.jl b/examples/DiffPrepMCMRF.jl
index b891c81cc..0c3916f9f 100644
--- a/examples/DiffPrepMCMRF.jl
+++ b/examples/DiffPrepMCMRF.jl
@@ -158,9 +158,9 @@ end
 function write_diffprep_fwf(G1, G2, G3, bmax, Gmax, Smax; filename="./qte_vectors_input.txt", name="Maxwell2", 
     precision::Int=6, dwell_time=6.4e-6, verbose=false)
     open(filename, "w") do io
-        t1 = range(0, maximum(G1.GR.T), step=dwell_time) #length=δ2N(maximum(G1.GR.T))) #step=dwell_time) #
-		t2 = range(0, maximum(G2.GR.T), step=dwell_time) #length=δ2N(maximum(G2.GR.T)))
-        t3 = range(0, maximum(G3.GR.T), step=dwell_time) #length=δ2N(maximum(G3.GR.T)))
+        t1 = range(0, G1.GR.dur[1] - maximum(G1.GR.delay), step=dwell_time) #length=δ2N(maximum(G1.GR.T))) #step=dwell_time) #
+		t2 = range(0, G2.GR.dur[1] - maximum(G2.GR.delay), step=dwell_time) #length=δ2N(maximum(G2.GR.T)))
+        t3 = range(0, G3.GR.dur[1] - maximum(G3.GR.delay), step=dwell_time) #length=δ2N(maximum(G3.GR.T)))
         maxN = max(length(t1), length(t2), length(t3))
         Gx1, Gy1, Gz1 = KomaMRICore.get_grads(G1, Array(t1).+maximum(G1.GR.delay))
 		Gx2, Gy2, Gz2 = KomaMRICore.get_grads(G2, Array(t2).+maximum(G2.GR.delay))
@@ -242,7 +242,7 @@ function write_diffprep_fwf(G1, G2, G3, bmax, Gmax, Smax; filename="./qte_vector
         date = "#Generated on $(now())\r\n"
         vars =  @sprintf "%s %s %s %s %s %s %s\r\n" "#Name"*" "^(length(name)-5) "N1"*" "^(length(string(N1))-1) "N2"*" "^(length(string(N2))-1) "N3"*" "^(length(string(N3))-1) "bval"*" "^(length(string(round(bmax,digits=1)))-3) "Gmax"*" "^(length(string(round(Gmax,digits=1)))-3) "Smax"
         unit =  @sprintf "%s %s %s %s\r\n" "#"*" "^(length(name)+length(string(N1))+length(string(N2))+length(string(N3))+2)  "s/mm2"*" "^(length(string(bval))-3) "mT/m"*" "^(length(string(round(Gmax,digits=1)))-3) "T/m/s"  
-        line =  @sprintf "%s %i %i %i %.1f %.1f %.1f\r\n" name N1 N2 N3 bmax Gmax*1e3 Smax
+        line =  @sprintf "%s %i %i %i %.1f %.1f %.1f\r\n" name N1 N2 N3 bval Gmax*1e3 Smax
         write(io, date)
         write(io, vars)
         write(io, unit)
@@ -260,12 +260,13 @@ end
 #Params.
 dwell_time = 6.4e-6
 Gmax = 62e-3 # mT/m
-Smaxs = [70, 100] # mT/m/ms
-axis_to_calc = ["x", "y", "z", "xyz"] # ["x", "y", "z", "xyz", "xz", "yz"] 
+Smaxs = [70] # mT/m/ms
+axis_to_calc = ["xyz"] # ["x", "y", "z", "xyz", "xz", "yz"] 
 moment_to_calc = [0, 1] #[0, 1, 2]
+# 4 : HS 50ms
 # 5 : HS 55ms
 # 9 : BIR4 50ms
-pulses_to_calc = [5]
+pulses_to_calc = [4] #3,4,5]
 n_dwells = 4
 maxwell = true #Maxwell/concomitant gradient compensation
 gap_left_ms = 0; gap_left = floor(Int64, gap_left_ms * 1e-3 / (n_dwells * dwell_time))
@@ -278,8 +279,9 @@ ec_spectraGradM0 = zeros(length(λ_spectra))
 ecc_λ = [85e-3]#[0.04, 0.1, 0.24, 0.57, 1.4, 3.4, 8.3, 20, 49, 118, 288, 700] * 1e-3
 ecc_α = [1.] #[51.2*ones(3); 10.24*ones(9)] / 100.
 DIF = Sequence()
+DIF_ref = Sequence()
 
-for eddy = [false], gap_left_ms = [0]
+for eddy = [false], gap_left_ms = [0, 1], maxwell = [true, false]
 # for eddy = [false, true], gap_left_ms = [0, 1]
 gap_left = floor(Int64, gap_left_ms * 1e-3 / (n_dwells * dwell_time))
 for Smax = Smaxs
@@ -512,7 +514,9 @@ for k = moment_to_calc #Number of moments to null
         println( "NOT Solved 😢" )
     end
     ## Solution to Sequence object (for plotting)
-    B1 = 13.5e-6
+    t = range(-1.5*rf1/2, 1.5*rf1/2, 80)
+    β = 4e2 #frequency modulation param (rad/s)
+    B1 = 2*13.5e-6 * sech.(β * t)
     R1 = [RF(B1, rf1, 0, δ1);;]
     R2 = [RF(B1, rf2, 0, δ2);;]
     for axis = axis_to_calc
@@ -554,7 +558,7 @@ for k = moment_to_calc #Number of moments to null
             bmax = 3 * objective_value(model)
         end
         ## TO SCANNER
-        path_res = "C:/Users/carlo/Desktop/DPW/G$(floor(Int,Gmax*1e3))_SR$(ceil(Int,Smax))_$axis/"
+        path_res = "/home/ccp/DPW/G$(floor(Int,Gmax*1e3))_SR$(ceil(Int,Smax))_$axis/"
         inv = sum(DIF[1].GR[ax].A) <= 0 #if first grdient's x area is negative, invert 
         global DIF = inv ? -DIF : DIF
         # Plots
@@ -568,28 +572,51 @@ for k = moment_to_calc #Number of moments to null
         #     global ec_spectraB0[i] = get_ECmatrix(DIF; λ, τ_sample=τ*1e-3) * [gx1; gx2; gx3]
         #     global ec_spectraGradM0[i] = get_ECmatrixM0(DIF; λ, τ_sample=τ*1e-3) * [gx1; gx2; gx3]
         # end
-        # p5 = plot_eddy_currents(DIF, ecc_λ; α=ecc_α, slider=false, range=[0,τ])
+        p5 = plot_eddy_currents(DIF, ecc_λ; α=ecc_α, slider=false, range=[0,τ])
         # p6 = plot([
         #     scatter(x=λ_spectra*1e3, y=log.(abs.(ec_spectraB0)), name="EC_B0"), 
         #     scatter(x=λ_spectra*1e3, y=log.(abs.(ec_spectraGradM0)), name="EC_GradM0")
         #     ])
         # p7 = KomaMRIPlots.plot_slew_rate(DIF; slider=false, range=[0,τ])
-        p = [p1; p2 p3 p4]
-        display(p)
-        #savefig(p, path_res*"$seq_name.svg")
+        p = [p1; p2 p3; p5]
+        # display(p)
+        savefig(p, path_res*"$seq_name.svg")
         # Write
         write_diffprep_fwf(DIF[1], DIF[2], DIF[3], bmax, Gmax, Smax; 
                 filename=path_res*"$seq_name.txt", name=seq_name, verbose=false)
-        
         println( "λ0 = $(abs(round(M[1,:]'*gx/Gmax,digits=3))), λ1 = $(abs(round(M[2,:]'*gx/Gmax,digits=3))), λ2 = $(abs(round(M[3,:]'*gx/Gmax,digits=3)))" )
         println( "MX ∫g1²-∫g2²+∫g3²=$(gx1'*MX[1]*gx1 - gx2'*MX[2]*gx2 + gx3'*MX[3]*gx3)")
         println( "b-value: $(round(bmax, digits=2)) s/mm2" )
         println( "Eddy currents B0: $(ecc_α'*EC*[gx1; gx2; gx3])")
         println( "Eddy currents GradM0: $(ecc_α'*ECM0*[gx1; -gx2; gx3])")
         println( seq_name )
+
+        #REFERENCES
+        #REF M0 TRSE
+        ζ = Gmax/Smax_discrete
+        G_waveform = [-Gmax; -Gmax; 0; Gmax; Gmax]
+        t_waveform = [(δ2-4ζ)/2; ζ; ζ; (δ2-4ζ)/2]
+        println("TRSE: $((δ1-2ζ)*1e3) $(t_waveform*1e3) $((δ3-2ζ)*1e3)")
+        global DIF_ref = Sequence([Grad(Gmax, δ1-2ζ, ζ); Grad(Gmax, δ1-2ζ, ζ); Grad(Gmax, δ1-2ζ, ζ);;], R1)
+        global DIF_ref += Sequence([Grad(G_waveform, t_waveform, ζ); Grad(G_waveform, t_waveform, ζ); Grad(G_waveform, t_waveform, ζ);;], R2)
+        global DIF_ref += Sequence([Grad(-Gmax, δ3-2ζ, ζ); Grad(-Gmax, δ3-2ζ, ζ); Grad(-Gmax, δ3-2ζ, ζ);;])
+        #Plot
+        τ = dur(DIF_ref) * 1e3
+        p1 = plot_seq(DIF_ref; slider=false, range=[0,τ], title="TRSE_$durT")
+        p2 = plot_M0(DIF_ref; slider=false, range=[0,τ])
+        p3 = plot_M1(DIF_ref; slider=false, range=[0,τ])
+        p4 = plot_M2(DIF_ref; slider=false, range=[0,τ])
+        p5 = plot_eddy_currents(DIF_ref, ecc_λ; α=ecc_α, slider=false, range=[0,τ])
+        p = [p1; p2 p3; p5]
+        # display(p)
+        savefig(p, path_res*"TRSE_$durT.svg")
+        #Write
+        write_diffprep_fwf(DIF_ref[1], DIF_ref[2], DIF_ref[3], bmax, Gmax, Smax; 
+            filename=path_res*"TRSE_$durT.txt", name="TRSE_$durT", verbose=false)
     end
 end
 end
 end
 end
-println("Finished! 💃")
+
+println("Finished! 💃")
\ No newline at end of file
diff --git a/src/KomaUI.jl b/src/KomaUI.jl
index eef46c9f5..48f378e3e 100644
--- a/src/KomaUI.jl
+++ b/src/KomaUI.jl
@@ -1,6 +1,6 @@
 include("ui/ExportMATFunctions.jl")
 
-function KomaUI(;dark=true,frame=true, phantom_mode="2D", sim=Dict{String,Any}(), rec=Dict{Symbol,Any}(), dev_tools=false)
+function KomaUI(;dark=true,frame=true, phantom_mode="2D", sim=Dict{String,Any}(), rec=Dict{Symbol,Any}(), dev_tools=false, blink_show=true)
 ## ASSETS
 path = @__DIR__
 assets = AssetRegistry.register(dirname(path*"/ui/assets/"))
@@ -39,7 +39,8 @@ global w = Blink.Window(Dict(
     :icon=>path*"/ui/assets/Logo_icon.png",
     "width"=>1200,
     "height"=>800,
-    "webPreferences" => Dict("devTools" => dev_tools)
+    "webPreferences" => Dict("devTools" => dev_tools),
+    :show=>blink_show
     ),async=false);
 ## LOADING BAR
 buffericon = """<div class="spinner-border spinner-border-sm text-light" role="status"></div>"""
@@ -115,7 +116,7 @@ global raw_ismrmrd = RawAcquisitionData(Dict(
 global rawfile = ""
 global image =  [0.0im 0.; 0. 0.]
 global kspace = [0.0im 0.; 0. 0.]
-global matfolder = pwd()
+global matfolder = tempdir()
 #Reco
 default = Dict{Symbol,Any}(:reco=>"direct") #, :iterations=>10, :λ=>1e-5,:solver=>"admm",:regularization=>"TV")
 global recParams = merge(default, rec)
@@ -132,9 +133,42 @@ global sig_obs = Observable{RawAcquisitionData}(raw_ismrmrd)
 global img_obs = Observable{Any}(image)
 global mat_obs = Observable{Any}(matfolder)
 
+#
+handle(w, "matfolder") do args...
+    str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="all", matfilename="")
+    @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);)
+    @js_ w document.getElementById("content").dataset.content = "matfolder"
+end
+handle(w, "matfolderseq") do args...
+    str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="sequence")
+    @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);)
+    @js_ w document.getElementById("content").dataset.content = "matfolderseq"
+end
+handle(w, "matfolderpha") do args...
+    str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="phantom")
+    @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);)
+    @js_ w document.getElementById("content").dataset.content = "matfolderpha"
+end
+handle(w, "matfoldersca") do args...
+    str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="scanner")
+    @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);)
+    @js_ w document.getElementById("content").dataset.content = "matfoldersca"
+end
+handle(w, "matfolderraw") do args...
+    str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="raw")
+    @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);)
+    @js_ w document.getElementById("content").dataset.content = "matfolderraw"
+end
+handle(w, "matfolderima") do args...
+    str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="image")
+    @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);)
+    @js_ w document.getElementById("content").dataset.content = "matfolderima"
+end
+
 ## MENU FUNCTIONS
 handle(w, "index") do args...
     content!(w, "div#content", index)
+    @js_ w document.getElementById("content").dataset.content = "index"
 end
 handle(w, "pulses_seq") do args...
     loading = replace(open(f->read(f, String), path*"/ui/html/loading.html"), "LOADDES"=>"Plotting sequence ...")
@@ -171,11 +205,13 @@ handle(w, "sig") do args...
     loading = replace(open(f->read(f, String), path*"/ui/html/loading.html"), "LOADDES"=>"Plotting raw signal ...")
     content!(w, "div#content", loading)
     include(path*"/ui/SignalGUI.jl")
+    @js_ w document.getElementById("content").dataset.content = "sig"
 end
 handle(w, "reconstruction_absI") do args...
     loading = replace(open(f->read(f, String), path*"/ui/html/loading.html"), "LOADDES"=>"Plotting image magnitude ...")
     content!(w, "div#content", loading)
     include(path*"/ui/ReconGUI_absI.jl")
+    @js_ w document.getElementById("content").dataset.content = "absi"
 end
 handle(w, "reconstruction_angI") do args...
     loading = replace(open(f->read(f, String), path*"/ui/html/loading.html"), "LOADDES"=>"Plotting image phase ...")
@@ -244,6 +280,7 @@ handle(w, "simulate") do args...
     loading = replace(open(f->read(f, String), path*"/ui/html/loading.html"), "LOADDES"=>"Plotting raw signal ...")
     content!(w, "div#content", loading)
     include(path*"/ui/SignalGUI.jl")
+    @js_ w document.getElementById("content").dataset.content = "simulation"
     # @js_ w document.getElementById("simulate!").prop("disabled", false); #Re-enable button
     # @js_ w (@var button = document.getElementById("recon!"); @var bsButton = @new bootstrap.Button(button); vsButton.toggle())
 end
@@ -289,6 +326,7 @@ handle(w, "recon") do args...
     loading = replace(open(f->read(f, String), path*"/ui/html/loading.html"), "LOADDES"=>"Plotting image magnitude ...")
     content!(w, "div#content", loading)
     include(path*"/ui/ReconGUI_absI.jl")
+    @js_ w document.getElementById("content").dataset.content = "reconstruction"
 end
 handle(w, "close") do args...
     global darkmode = nothing
@@ -406,78 +444,6 @@ map!(f->if f!="" #Assigning function of data when load button (filepicker) is ch
         end
     , sig_obs, load_sig)
 w = content!(w, "#sigfilepicker", load_sig, async=false)
-# Folder observable
-load_folder = opendialog(; label = "Save All", properties = ["openDirectory"], icon = "far fa-save")
-map!(f->if f!="" #Assigning function of data when load button (opendialog) is changed
-            global matfolder = f[1]
-            str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="all", matfilename="")
-            @js_ w Toasty("1", "Saved .mat files" , str_toast);
-            matfolder
-        else
-            matfolder #default sequence
-        end
-    , mat_obs, load_folder)
-w = content!(w, "#matfolder", load_folder, async=false)
-
-load_folder_seq = savedialog(; label = "Sequence", defaultPath = "seq.mat", filters = [(; name = "Matlab Data", extensions = ["mat"])])
-map!(f->if f!="" #Assigning function of data when load button (opendialog) is changed
-            global matfolder = dirname(f)
-            str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="sequence", matfilename=basename(f))
-            @js_ w Toasty("1", "Saved .mat files" , str_toast);
-            matfolder
-        else
-            matfolder #default sequence
-        end
-    , mat_obs, load_folder_seq)
-w = content!(w, "#matfolderseq", load_folder_seq, async=false)
-
-load_folder_pha = savedialog(; label = "Phantom", defaultPath = "phantom.mat", filters = [(; name = "Matlab Data", extensions = ["mat"])])
-map!(f->if f!="" #Assigning function of data when load button (opendialog) is changed
-            global matfolder = dirname(f)
-            str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="phantom", matfilename=basename(f))
-            @js_ w Toasty("1", "Saved .mat files" , str_toast);
-            matfolder
-        else
-            matfolder #default sequence
-        end
-    , mat_obs, load_folder_pha)
-w = content!(w, "#matfolderpha", load_folder_pha, async=false)
-
-load_folder_sca = savedialog(; label = "Scanner", defaultPath = "scanner.mat", filters = [(; name = "Matlab Data", extensions = ["mat"])])
-map!(f->if f!="" #Assigning function of data when load button (opendialog) is changed
-            global matfolder = dirname(f)
-            str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="scanner", matfilename=basename(f))
-            @js_ w Toasty("1", "Saved .mat files" , str_toast);
-            matfolder
-        else
-            matfolder #default sequence
-        end
-    , mat_obs, load_folder_sca)
-w = content!(w, "#matfoldersca", load_folder_sca, async=false)
-
-load_folder_raw = savedialog(; label = "Raw", defaultPath = "raw.mat", filters = [(; name = "Matlab Data", extensions = ["mat"])])
-map!(f->if f!="" #Assigning function of data when load button (opendialog) is changed
-            global matfolder = dirname(f)
-            str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="raw", matfilename=basename(f))
-            @js_ w Toasty("1", "Saved .mat files" , str_toast);
-            matfolder
-        else
-            matfolder #default sequence
-        end
-    , mat_obs, load_folder_raw)
-w = content!(w, "#matfolderraw", load_folder_raw, async=false)
-
-load_folder_ima = savedialog(; label = "Image", defaultPath = "image.mat", filters = [(; name = "Matlab Data", extensions = ["mat"])])
-map!(f->if f!="" #Assigning function of data when load button (opendialog) is changed
-            global matfolder = dirname(f)
-            str_toast = export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="image", matfilename=basename(f))
-            @js_ w Toasty("1", "Saved .mat files" , str_toast);
-            matfolder
-        else
-            matfolder #default sequence
-        end
-    , mat_obs, load_folder_ima)
-w = content!(w, "#matfolderima", load_folder_ima, async=false)
 
 #Update Koma version
 version = string(KomaMRICore.__VERSION__)
diff --git a/src/ui/ExportMATFunctions.jl b/src/ui/ExportMATFunctions.jl
index ac539ac3a..66265f0d7 100644
--- a/src/ui/ExportMATFunctions.jl
+++ b/src/ui/ExportMATFunctions.jl
@@ -70,18 +70,18 @@ end
 
 function export_2_mat_raw(raw_ismrmrd, matfolder; matfilename="raw.mat");
     if haskey(raw_ismrmrd.params, "userParameters")
-        dictForMat = Dict()
-        dictUserParams = raw_ismrmrd.params["userParameters"]
-        for (key, value) in dictUserParams
+        dict_for_mat = Dict()
+        dict_user_params = raw_ismrmrd.params["userParameters"]
+        for (key, value) in dict_user_params
             if key == "Δt_rf"
-                dictForMat["dt_rf"] = value
+                dict_for_mat["dt_rf"] = value
             elseif key == "Δt"
-                dictForMat["dt_rf"] = value
+                dict_for_mat["dt_rf"] = value
             else
-                dictForMat[key] = value
+                dict_for_mat[key] = value
             end
         end
-        matwrite(joinpath(matfolder, "sim_params.mat"), dictForMat)
+        matwrite(joinpath(matfolder, "sim_params.mat"), dict_for_mat)
 
         not_Koma = raw_ismrmrd.params["systemVendor"] != "KomaMRI.jl"
         t = Float64[]
@@ -111,18 +111,19 @@ function export_2_mat_raw(raw_ismrmrd, matfolder; matfilename="raw.mat");
 
 end
 
-function export_2_mat_image(image, recParams, matfolder; matfilename="image.mat")
-    if haskey(recParams, :reconSize)
-        recParams_dict = Dict("reco" => recParams[:reco],
-                            "Nx" => recParams[:reconSize][1],
-                            "Ny" => recParams[:reconSize][2])
-        matwrite(joinpath(matfolder, "rec_params.mat"), Dict("rec_params" => recParams_dict))
+function export_2_mat_image(image, rec_params, matfolder; matfilename="image.mat")
+    if haskey(rec_params, :reconSize)
+        dict_rec_params = Dict("reco" => rec_params[:reco],
+                            "Nx" => rec_params[:reconSize][1],
+                            "Ny" => rec_params[:reconSize][2])
+        matwrite(joinpath(matfolder, "rec_params.mat"), Dict("rec_params" => dict_rec_params))
     end
 
     matwrite(joinpath(matfolder, matfilename), Dict("image" => image))
 end
 
-function export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolder; type="all", matfilename="data.mat")
+function export_2_mat(seq, phantom, sys, raw_ismrmrd, rec_params, image, matfolder; type="all", matfilename="data.mat")
+    head = splitext(matfilename)[1]
     if type=="all"
         export_2_mat_sequence(seq, matfolder)
         export_2_mat_kspace(seq, matfolder)
@@ -130,20 +131,19 @@ function export_2_mat(seq, phantom, sys, raw_ismrmrd, recParams, image, matfolde
         export_2_mat_phantom(phantom, matfolder)
         export_2_mat_scanner(sys, matfolder)
         export_2_mat_raw(raw_ismrmrd, matfolder)
-        export_2_mat_image(image, recParams, matfolder)
+        export_2_mat_image(image, rec_params, matfolder)
     elseif type=="sequence"
-        head = splitext(matfilename)[1]
 		export_2_mat_sequence(seq, matfolder; matfilename=(head*"_sequence.mat"))
         export_2_mat_kspace(seq, matfolder; matfilename=(head*"_kspace.mat"))
         export_2_mat_moments(seq, matfolder; matfilename=(head*"_moments.mat"))
 	elseif type=="phantom"
-		export_2_mat_phantom(phantom, matfolder; matfilename)
+		export_2_mat_phantom(phantom, matfolder; matfilename=(head*"_phantom.mat"))
     elseif type=="scanner"
-		export_2_mat_scanner(sys, matfolder; matfilename)
+		export_2_mat_scanner(sys, matfolder; matfilename=(head*"_scanner.mat"))
     elseif type=="raw"
-		export_2_mat_raw(raw_ismrmrd, matfolder; matfilename)
+		export_2_mat_raw(raw_ismrmrd, matfolder; matfilename=(head*"_raw.mat"))
     elseif type=="image"
-		export_2_mat_image(image, recParams, matfolder; matfilename)
+		export_2_mat_image(image, rec_params, matfolder; matfilename=(head*"_image.mat"))
 	end
 
     strToast = "<ul><li><b>Name:</b> " * matfilename * "</li><li><b>Path:</b> " * matfolder * "</li></ul>"
diff --git a/src/ui/PhantomGUI.jl b/src/ui/PhantomGUI.jl
index 500df4efd..fb47e1115 100644
--- a/src/ui/PhantomGUI.jl
+++ b/src/ui/PhantomGUI.jl
@@ -25,3 +25,4 @@ columnbuttons = Observable{Any}(makebuttons(phantom))
 map!(makebuttons, columnbuttons, pha_obs)
 ui = dom"div"(vbox(columnbuttons, plt))
 content!(w, "div#content", ui)
+@js_ w document.getElementById("content").dataset.content = "phantom"
diff --git a/src/ui/PulsesGUI_M0.jl b/src/ui/PulsesGUI_M0.jl
index e6cb2149e..6d9b0477e 100644
--- a/src/ui/PulsesGUI_M0.jl
+++ b/src/ui/PulsesGUI_M0.jl
@@ -2,3 +2,4 @@ plt = Observable{Any}(plot_M0(seq;darkmode))
 ui = dom"div"(plt)
 map!(p->plot_M0(p;darkmode), plt, seq_obs)
 content!(w, "div#content", ui)
+@js_ w document.getElementById("content").dataset.content = "m0"
diff --git a/src/ui/PulsesGUI_M1.jl b/src/ui/PulsesGUI_M1.jl
index c654271f3..93e08d2af 100644
--- a/src/ui/PulsesGUI_M1.jl
+++ b/src/ui/PulsesGUI_M1.jl
@@ -2,3 +2,4 @@ plt = Observable{Any}(plot_M1(seq;darkmode))
 ui = dom"div"(plt)
 map!(p->plot_M1(p;darkmode), plt, seq_obs)
 content!(w, "div#content", ui)
+@js_ w document.getElementById("content").dataset.content = "m1"
diff --git a/src/ui/PulsesGUI_M2.jl b/src/ui/PulsesGUI_M2.jl
index 451987c8b..d10690c8c 100644
--- a/src/ui/PulsesGUI_M2.jl
+++ b/src/ui/PulsesGUI_M2.jl
@@ -2,3 +2,4 @@ plt = Observable{Any}(plot_M2(seq;darkmode))
 ui = dom"div"(plt)
 map!(p->plot_M2(p;darkmode), plt, seq_obs)
 content!(w, "div#content", ui)
+@js_ w document.getElementById("content").dataset.content = "m2"
diff --git a/src/ui/PulsesGUI_kspace.jl b/src/ui/PulsesGUI_kspace.jl
index 49859fd2a..f539fde40 100644
--- a/src/ui/PulsesGUI_kspace.jl
+++ b/src/ui/PulsesGUI_kspace.jl
@@ -2,3 +2,4 @@ plt = Observable{Any}(plot_kspace(seq;darkmode))
 ui = dom"div"(plt)
 map!(p->plot_kspace(p;darkmode), plt, seq_obs)
 content!(w, "div#content", ui)
+@js_ w document.getElementById("content").dataset.content = "kspace"
diff --git a/src/ui/PulsesGUI_seq.jl b/src/ui/PulsesGUI_seq.jl
index 634ce22fa..1b0d283e1 100644
--- a/src/ui/PulsesGUI_seq.jl
+++ b/src/ui/PulsesGUI_seq.jl
@@ -2,3 +2,4 @@ plt = Observable{Any}(plot_seq(seq; darkmode, range=[0 30]))
 ui = dom"div"(plt)
 map!(p->plot_seq(p; darkmode, range=[0 30]), plt, seq_obs)
 content!(w, "div#content", ui)
+@js_ w document.getElementById("content").dataset.content = "sequence"
diff --git a/src/ui/RecParams_view.jl b/src/ui/RecParams_view.jl
index 8312738e8..9408a260e 100644
--- a/src/ui/RecParams_view.jl
+++ b/src/ui/RecParams_view.jl
@@ -1,3 +1,4 @@
 ui = plot_dict(recParams)
 title = """<h1 style="padding: 8px 16px; color: #868888;">Reconstruction parameters</h1>"""
 content!(w, "div#content", title*ui)
+@js_ w document.getElementById("content").dataset.content = "recparams"
diff --git a/src/ui/ReconGUI_absK.jl b/src/ui/ReconGUI_absK.jl
index 2f719e9cf..e74993f7c 100644
--- a/src/ui/ReconGUI_absK.jl
+++ b/src/ui/ReconGUI_absK.jl
@@ -9,3 +9,4 @@ plt = Observable{Any}(plotInteract())
 map!(t-> plotInteract(), plt, img_obs)
 ui = dom"div"(plt)
 content!(w, "div#content", ui)
+@js_ w document.getElementById("content").dataset.content = "absk"
diff --git a/src/ui/ReconGUI_angI.jl b/src/ui/ReconGUI_angI.jl
index f44c27347..ce1162858 100644
--- a/src/ui/ReconGUI_angI.jl
+++ b/src/ui/ReconGUI_angI.jl
@@ -9,3 +9,4 @@ plt = Observable{Any}(plotInteract())
 map!(t-> plotInteract(), plt, img_obs)
 ui = dom"div"(plt)
 content!(w, "div#content", ui)
+@js_ w document.getElementById("content").dataset.content ="angi"
diff --git a/src/ui/ScannerParams_view.jl b/src/ui/ScannerParams_view.jl
index 7edefe493..11b14a468 100644
--- a/src/ui/ScannerParams_view.jl
+++ b/src/ui/ScannerParams_view.jl
@@ -12,3 +12,4 @@ sys_dict = Dict("B0" => sys.B0,
 plt = plot_dict(sys_dict)
 title = """<h1 style="padding: 8px 16px; color: #868888;">Scanner parameters</h1>"""
 content!(w, "div#content", title*plt)
+@js_ w document.getElementById("content").dataset.content = "scanneparams"
diff --git a/src/ui/SimParams_view.jl b/src/ui/SimParams_view.jl
index 7b05362d7..2ca5f94b1 100644
--- a/src/ui/SimParams_view.jl
+++ b/src/ui/SimParams_view.jl
@@ -1,3 +1,4 @@
 plt = plot_dict(simParams)
 title = """<h1 style="padding: 8px 16px; color: #868888;">Simulation parameters</h1>"""
 content!(w, "div#content", title*plt)
+@js_ w document.getElementById("content").dataset.content = "simparams"
diff --git a/src/ui/css/custom.css b/src/ui/css/custom.css
index 84a16feb4..61c5ac32d 100644
--- a/src/ui/css/custom.css
+++ b/src/ui/css/custom.css
@@ -164,4 +164,7 @@ td, th, tr{
 
 .pname {
   margin-left: 24px;
+  overflow: hidden;
+  text-overflow: ellipsis;
+  white-space: nowrap;
 }
diff --git a/src/ui/html/sidebar.html b/src/ui/html/sidebar.html
index 0385e7113..bf5e51582 100644
--- a/src/ui/html/sidebar.html
+++ b/src/ui/html/sidebar.html
@@ -11,18 +11,18 @@
           <span class="icon-wave-square me-2"></span><span class="w-100">Sequence</span>
       </a>
         <ul class="dropdown-menu dropdown-menu-dark shadow" aria-labelledby="seq-dropdown" style="position: absolute;">
-          <li><a class="dropdown-item" href="#" onclick='Blink.msg("pulses_seq", 1)'>View Sequence</a></li>
+          <li><a id="button_pulses_seq" class="dropdown-item" href="#" onclick='Blink.msg("pulses_seq", 1)'>View Sequence</a></li>
           <!-- Sequence (MPS) -->
           <!--<li><a class="dropdown-item disabled" href="#" onclick='Blink.msg("pulses_seqxyz", 1)'>Sequence (XYZ)</a></li>-->
-          <li><a class="dropdown-item" href="#" onclick='Blink.msg("pulses_kspace", 1)'>View $k$-space</a></li>
+          <li><a id="button_pulses_kspace" class="dropdown-item" href="#" onclick='Blink.msg("pulses_kspace", 1)'>View $k$-space</a></li>
           <li>
             <a class="dropdown-item" href="#">
               View Moments &raquo;
             </a>
             <ul class="dropdown-menu dropdown-menu-dark dropdown-submenu">
-              <li><a class="dropdown-item" href="#" onclick='Blink.msg("pulses_M0", 1)'>$M_0(t)$</a></li>
-              <li><a class="dropdown-item" href="#" onclick='Blink.msg("pulses_M1", 1)'>$M_1(t)$</a></li>
-              <li><a class="dropdown-item" href="#" onclick='Blink.msg("pulses_M2", 1)'>$M_2(t)$</a></li>
+              <li><a id="button_pulses_M0" class="dropdown-item" href="#" onclick='Blink.msg("pulses_M0", 1)'>$M_0(t)$</a></li>
+              <li><a id="button_pulses_M1" class="dropdown-item" href="#" onclick='Blink.msg("pulses_M1", 1)'>$M_1(t)$</a></li>
+              <li><a id="button_pulses_M2" class="dropdown-item" href="#" onclick='Blink.msg("pulses_M2", 1)'>$M_2(t)$</a></li>
             </ul>
           </li>
           <li><hr class="dropdown-divider"></li>
@@ -36,7 +36,7 @@
           <span class="icon-brain me-2"></span><span class="w-100">Phantom</span>
       </a>
       <ul class="dropdown-menu dropdown-menu-dark shadow" aria-labelledby="pha-dropdown">
-          <li><a class="dropdown-item" href="#" onclick='Blink.msg("phantom", 1)'>View Phantom</a></li>
+          <li><a id="button_phantom" class="dropdown-item" href="#" onclick='Blink.msg("phantom", 1)'>View Phantom</a></li>
           <li><hr class="dropdown-divider"></li>
           <li class="dropdown-interact" id="phafilepicker"></li>
       </ul>
@@ -48,7 +48,7 @@
         <span class="icon-magnet me-2"></span><span class="w-100">Scanner</span>
       </a>
       <ul class="dropdown-menu dropdown-menu-dark shadow" aria-labelledby="pha-dropdown">
-        <li><a class="dropdown-item" href="#" onclick='Blink.msg("scanner", 1)'>View Scanner</a></li>
+        <li><a id="button_scanner" class="dropdown-item" href="#" onclick='Blink.msg("scanner", 1)'>View Scanner</a></li>
       </ul>
     </li>
     <p id="scaname" class="pname">default.scanner</p>
@@ -61,7 +61,7 @@
         </button>
         <div class="dropdown-menu dropdown-menu-dark shadow" aria-labelledby="btnDropdownDemo">
           <a class="dropdown-item disabled" href="#"><i class="bi bi-upload me-2"></i>Load Options</a>
-          <a class="dropdown-item" href="#" onclick='Blink.msg("sim_params", 1)'><i class="bi bi-gear-fill me-2"></i>View Options</a>
+          <a class="dropdown-item" href="#" id="button_sim_params" onclick='Blink.msg("sim_params", 1)'><i class="bi bi-gear-fill me-2"></i>View Options</a>
         </div>
       </div>
     </li>
@@ -71,7 +71,7 @@
           <span class="icon-oscilloscope me-2"></span><span class="w-100">RawData</span>
       </a>
       <ul class="dropdown-menu dropdown-menu-dark shadow" aria-labelledby="sig-dropdown">
-          <li><a class="dropdown-item" href="#" onclick='Blink.msg("sig", 1)'>View RawData</a></li>
+          <li><a id="button_sig" class="dropdown-item" href="#" onclick='Blink.msg("sig", 1)'>View RawData</a></li>
           <li><hr class="dropdown-divider"></li>
           <li class="dropdown-interact" id="sigfilepicker"></li>
       </ul>
@@ -86,7 +86,7 @@
         </button>
         <div class="dropdown-menu dropdown-menu-dark shadow" aria-labelledby="btnDropdownDemo">
           <a class="dropdown-item disabled" href="#"><i class="bi bi-upload me-2"></i>Load Options</a>
-          <a class="dropdown-item" href="#"  onclick='Blink.msg("rec_params", 1)'><i class="bi bi-gear-fill me-2"></i>View Options</a>
+          <a class="dropdown-item" href="#" id="button_rec_params" onclick='Blink.msg("rec_params", 1)'><i class="bi bi-gear-fill me-2"></i>View Options</a>
         </div>
       </div>
     </li>
@@ -96,9 +96,9 @@
           <span class="icon-images me-2"></span><span class="w-100">Image</span>
       </a>
       <ul class="dropdown-menu dropdown-menu-dark shadow" aria-labelledby="rec-dropdown">
-          <li><a class="dropdown-item" href="#" onclick='Blink.msg("reconstruction_absI", 1)'>View $|\mathrm{Image}|$</a></li>
-          <li><a class="dropdown-item" href="#" onclick='Blink.msg("reconstruction_angI", 1)'>View $\angle\mathrm{Image}$</a></li>
-          <li><a class="dropdown-item" href="#" onclick='Blink.msg("reconstruction_absK", 1)'>View $\log(|k$-data$|+1)$</a></li>
+          <li><a id="button_reconstruction_absI" class="dropdown-item" href="#" onclick='Blink.msg("reconstruction_absI", 1)'>View $|\mathrm{Image}|$</a></li>
+          <li><a id="button_reconstruction_angI" class="dropdown-item" href="#" onclick='Blink.msg("reconstruction_angI", 1)'>View $\angle\mathrm{Image}$</a></li>
+          <li><a id="button_reconstruction_absK" class="dropdown-item" href="#" onclick='Blink.msg("reconstruction_absK", 1)'>View $\log(|k$-data$|+1)$</a></li>
       </ul>
     </li>
     <hr>
@@ -108,13 +108,21 @@
           <span class="bi-download me-2"></span><span class="w-100">Export to .mat</span>
       </a>
       <ul class="dropdown-menu dropdown-menu-dark shadow" aria-labelledby="exp-dropdown">
-          <li><a class="dropdown-interact" id="matfolder" href="#" >All</a></li>
+          <!--<li><a class="dropdown-interact" id="matfolder" href="#" >All</a></li>-->
+          <li><a id="button_matfolder" class="dropdown-item" href="#" onclick='Blink.msg("matfolder", 1)'>All</a></li>
           <li><hr class="dropdown-divider"></li>
+          <li><a id="button_matfolderseq" class="dropdown-item" href="#" onclick='Blink.msg("matfolderseq", 1)'>Sequence</a></li>
+          <li><a id="button_matfolderpha" class="dropdown-item" href="#" onclick='Blink.msg("matfolderpha", 1)'>Phantom</a></li>
+          <li><a id="button_matfoldersca" class="dropdown-item" href="#" onclick='Blink.msg("matfoldersca", 1)'>Scanner</a></li>
+          <li><a id="button_matfolderraw" class="dropdown-item" href="#" onclick='Blink.msg("matfolderraw", 1)'>RawData</a></li>
+          <li><a id="button_matfolderima" class="dropdown-item" href="#" onclick='Blink.msg("matfolderima", 1)'>Image</a></li>
+          <!--
           <li><a class="dropdown-interact" id="matfolderseq" href="#">Sequence</a></li>
           <li><a class="dropdown-interact" id="matfolderpha" href="#">Phantom</a></li>
           <li><a class="dropdown-interact" id="matfoldersca" href="#">Scanner</a></li>
           <li><a class="dropdown-interact" id="matfolderraw" href="#">RawData</a></li>
           <li><a class="dropdown-interact" id="matfolderima" href="#">Image</a></li>
+          -->
       </ul>
     </li>
   </ul>
@@ -197,4 +205,4 @@
   }
 </script>
 <!-- For index.html -->
-<div class="container" style="padding: 0px !important;" id="content">
\ No newline at end of file
+<div class="container" style="padding: 0px !important;" id="content" data-content="index">
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index dcfe93079..26149747f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -36,114 +36,120 @@ using TestItems, TestItemRunner
 end
 
 @testitem "KomaUI" tags=[:koma] begin
+
     using Blink, Interact
-    include(joinpath(@__DIR__, "../src/reconstruction/Recon.jl"))
-
-    w = Window(Blink.Dict(:show => false), async=false);
-    body!(w, """<div class="container" style="padding: 0px !important;" id="content">""");
-
-    phantom = brain_phantom2D()
-    sys = Scanner()
-    B1 = sys.B1; durRF = π/2/(2π*γ*B1) #90-degree hard excitation pulse
-    EX = PulseDesigner.RF_hard(B1, durRF, sys; G=[0,0,0])
-    N = 101
-    FOV = 23e-2
-    EPI = PulseDesigner.EPI(FOV, N, sys)
-    TE = 30e-3
-    d1 = TE-dur(EPI)/2-dur(EX)
-    d1 = d1 > 0 ? d1 : 0
-    if d1 > 0 DELAY = Delay(d1) end
-    seq = d1 > 0 ? EX + DELAY + EPI : EX + EPI #Only add delay if d1 is positive (enough space)
-    seq.DEF["TE"] = round(d1 > 0 ? TE : TE - d1, digits=4)*1e3
-    simParams = KomaMRICore.default_sim_params()
 
-    path = @__DIR__
-    fraw = ISMRMRDFile(path*"/test_files/Koma_signal.mrd")
-    darkmode = true
-    raw_ismrmrd = RawAcquisitionData(fraw)
+    function with_timeout(f::Function, timeout)
+        c = Channel{Any}(1)
+        @async begin
+            put!(c, f())
+        end
+        @async begin
+            sleep(timeout)
+            put!(c, nothing)
+        end
+        take!(c)
+    end
 
-    acq = AcquisitionData(raw_ismrmrd)
-    Nx, Ny = raw_ismrmrd.params["reconSize"][1:2]
-    recParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(Nx,Ny), :densityWeighting=>true)
-    aux = reconstruction(acq, recParams)
-    image  = reshape(aux.data, Nx, Ny, :)
-    kspace = KomaMRI.fftc(image)
+    function attempt_to_open_koma_ui(n_attempts, timeout_sec)
+        for cnt = 1:n_attempts
+            @info "Trying to open the KomaUI-Window ..."
+            w = with_timeout(()->KomaUI(dev_tools=true), timeout_sec)
+            @info "Number of KomaUI-Window attempts: $cnt"
+            if !isnothing(w)
+                @info "KomaUI-Window successfully opened"
+                return w
+            end
+        end
+    end
+
+    n_attempts = 5
+    timeout_sec = 120
+    w = attempt_to_open_koma_ui(n_attempts, timeout_sec)
 
-    seq_obs = Observable{Sequence}(seq)
-    pha_obs = Observable{Phantom}(phantom)
-    sig_obs = Observable{RawAcquisitionData}(raw_ismrmrd)
-    img_obs = Observable{Any}(image)
+    @testset "Open UI" begin
+        @test "index" == @js w document.getElementById("content").dataset.content
+    end
 
     @testset "PulsesGUI" begin
-        include(joinpath(@__DIR__, "../src/ui/PulsesGUI_seq.jl"))
-        include(joinpath(@__DIR__, "../src/ui/PulsesGUI_kspace.jl"))
-        include(joinpath(@__DIR__, "../src/ui/PulsesGUI_M0.jl"))
-        include(joinpath(@__DIR__, "../src/ui/PulsesGUI_M1.jl"))
-        include(joinpath(@__DIR__, "../src/ui/PulsesGUI_M2.jl"))
-        @test true
+        @js w document.getElementById("button_pulses_seq").click()
+        @test "sequence" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_pulses_kspace").click()
+        @test "kspace" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_pulses_M0").click()
+        @test "m0" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_pulses_M1").click()
+        @test "m1" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_pulses_M2").click()
+        @test "m2" == @js w document.getElementById("content").dataset.content
     end
 
     @testset "PhantomGUI" begin
-        include(joinpath(@__DIR__, "../src/ui/PhantomGUI.jl"))
-        @test true
+        @js w document.getElementById("button_phantom").click()
+        @test "phantom" == @js w document.getElementById("content").dataset.content
     end
 
     @testset "ParamsGUI" begin
-        include(joinpath(@__DIR__, "../src/ui/ScannerParams_view.jl"))
-        include(joinpath(@__DIR__, "../src/ui/SimParams_view.jl"))
-        include(joinpath(@__DIR__, "../src/ui/RecParams_view.jl"))
-        @test true
+        @js w document.getElementById("button_scanner").click()
+        @test "scanneparams" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_sim_params").click()
+        @test "simparams" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_rec_params").click()
+        @test "recparams" == @js w document.getElementById("content").dataset.content
     end
 
-    @testset "ReconGUI" begin
-        include(joinpath(@__DIR__, "../src/ui/ReconGUI_absI.jl"))
-        include(joinpath(@__DIR__, "../src/ui/ReconGUI_angI.jl"))
-        include(joinpath(@__DIR__, "../src/ui/ReconGUI_absK.jl"))
-        @test true
+    @testset "Simulation" begin
+        @js w document.getElementById("simulate!").click()
+        @test "simulation" == @js w document.getElementById("content").dataset.content
     end
 
     @testset "SignalGUI" begin
-        include(joinpath(@__DIR__, "../src/ui/SignalGUI.jl"))
-        @test true
+        @js w document.getElementById("button_sig").click()
+        @test "sig" == @js w document.getElementById("content").dataset.content
     end
 
-    @testset "Open UI" begin
-        KomaUI()
-        @test true
+    @testset "Reconstruction" begin
+        @js w document.getElementById("recon!").click()
+        @test "reconstruction" == @js w document.getElementById("content").dataset.content
     end
-end
 
-@testitem "Auxiliar Functions" tags=[:koma] begin
-    using MAT
-    include(joinpath(@__DIR__, "../src/ui/ExportMATFunctions.jl"))
-    @testset "ExportMATFunctions" begin
-        phantom = brain_phantom2D()
-        sys = Scanner()
-        B1 = sys.B1; durRF = π/2/(2π*γ*B1) #90-degree hard excitation pulse
-        EX = PulseDesigner.RF_hard(B1, durRF, sys; G=[0,0,0])
-        N = 101
-        FOV = 23e-2
-        EPI = PulseDesigner.EPI(FOV, N, sys)
-        TE = 30e-3
-        d1 = TE-dur(EPI)/2-dur(EX)
-        d1 = d1 > 0 ? d1 : 0
-        if d1 > 0 DELAY = Delay(d1) end
-        seq = d1 > 0 ? EX + DELAY + EPI : EX + EPI #Only add delay if d1 is positive (enough space)
-        seq.DEF["TE"] = round(d1 > 0 ? TE : TE - d1, digits=4)*1e3
-        path = @__DIR__
-        fraw = ISMRMRDFile(path*"/test_files/Koma_signal.mrd")
-        raw = RawAcquisitionData(fraw)
-        acq = AcquisitionData(raw)
-        Nx, Ny = raw.params["reconSize"][1:2]
-        recParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(Nx,Ny), :densityWeighting=>true)
-        aux = reconstruction(acq, recParams)
-        image  = reshape(aux.data, Nx, Ny, :)
-        export_2_mat(seq, phantom, sys, raw, recParams, image, pwd(); type="all")
-        export_2_mat(seq, phantom, sys, raw, recParams, image, pwd(); type="sequence")
-        export_2_mat(seq, phantom, sys, raw, recParams, image, pwd(); type="phantom")
-        export_2_mat(seq, phantom, sys, raw, recParams, image, pwd(); type="scanner")
-        export_2_mat(seq, phantom, sys, raw, recParams, image, pwd(); type="raw")
-        export_2_mat(seq, phantom, sys, raw, recParams, image, pwd(); type="image")
-        @test true
+    @testset "ReconGUI" begin
+        @js w document.getElementById("button_reconstruction_absI").click()
+        @test "absi" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_reconstruction_angI").click()
+        @test "angi" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_reconstruction_absK").click()
+        @test "absk" == @js w document.getElementById("content").dataset.content
     end
+
+    @testset "ExportToMAT" begin
+        @js w document.getElementById("button_matfolder").click()
+        @test "matfolder" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_matfolderseq").click()
+        @test "matfolderseq" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_matfolderpha").click()
+        @test "matfolderpha" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_matfoldersca").click()
+        @test "matfoldersca" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_matfolderraw").click()
+        @test "matfolderraw" == @js w document.getElementById("content").dataset.content
+
+        @js w document.getElementById("button_matfolderima").click()
+        @test "matfolderima" == @js w document.getElementById("content").dataset.content
+    end
+
+    close(w)
+
 end