BlochSim.BlochDynamicsMatrixType
BlochDynamicsMatrix(R1, R2, Δω)
BlochDynamicsMatrix{T}()
BlochDynamicsMatrix()

Create a mutable BlochDynamicsMatrix object.

Properties

• R1::Real: Spin-lattice relaxation rate
• R2::Real: Spin-spin relaxation rate
• Δω::Real: Off-resonance frequency
BlochSim.BlochMatrixType
BlochMatrix(a11, a21, a31, a12, a22, a32, a13, a23, a33)
BlochMatrix{T}()
BlochMatrix()

Create a mutable BlochMatrix object representing a fixed-size 3×3 matrix.

Properties

• aij::Real: Matrix entry (i,j) ∈ {1, 2, 3}²
BlochSim.BlochMcConnellDynamicsMatrixType
BlochMcConnellDynamicsMatrix(A, E)
BlochMcConnellDynamicsMatrix{T}(N)
BlochMcConnellDynamicsMatrix(N)

Create a BlochMcConnellDynamicsMatrix object with N compartments.

Properties

• A::NTuple{N,BlochDynamicsMatrix{<:Real}}: List of BlochDynamicsMatrixes that make up the main block diagonal of the BlochMcConnellDynamicsMatrix
• E::NTuple{M,ExchangeDynamicsMatrix{<:Real}}: List of ExchangeDynamicsMatrixes that describe exchange between the different compartments; these matrices make up the remaining M = N^2 - N blocks of the BlochMcConnellDynamicsMatrix, sorted by column-major ordering
BlochSim.BlochMcConnellMatrixType
BlochMcConnellMatrix(A)
BlochMcConnellMatrix{T}(N)
BlochMcConnellMatrix(N)

Create a BlochMcConnellMatrix object with N compartments and representing a fixed-size 3N×3N matrix.

Properties

• A::NTuple{N,NTuple{N,BlochMatrix{<:Real}}}: List of 3×3 matrices that comprise the blocks of the BlochMcConnellMatrix; A[i][j] is the (i,j)th block

Examples

julia> B = BlochMcConnellMatrix(3)
BlochMcConnellMatrix{Float64,3}:
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> fill!(B, 0.5)

julia> B
BlochMcConnellMatrix{Float64,3}:
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5
BlochSim.ExchangeDynamicsMatrixType
ExchangeDynamicsMatrix(r)
ExchangeDynamicsMatrix{T}()
ExchangeDynamicsMatrix()

Create a mutable ExchangeDynamicsMatrix object.

Properties

• r::Real: Exchange rate from one compartment to another
BlochSim.ExcitationMatrixType
ExcitationMatrix(A)
ExcitationMatrix{T}()
ExcitationMatrix()

Create an ExcitationMatrix object. Multiplying by a MagnetizationMC object has the effect of multiplying each component of the multi-compartment magnetization by the ExcitationMatrix.

Properties

• A::BlochMatrix{<:Real}: Matrix used to describe the excitation

Examples

julia> E = ExcitationMatrix(BlochMatrix(0, 1, 0, 1, 0, 0, 0, 0, 1))
ExcitationMatrix{Int64}:
0  1  0
1  0  0
0  0  1

julia> M = MagnetizationMC((1, 2, 3), (4, 5, 6))
2-compartment Magnetization vector with eltype Int64:
Compartment 1:
Mx = 1
My = 2
Mz = 3
Compartment 2:
Mx = 4
My = 5
Mz = 6

julia> E * M
2-compartment Magnetization vector with eltype Int64:
Compartment 1:
Mx = 2
My = 1
Mz = 3
Compartment 2:
Mx = 5
My = 4
Mz = 6
BlochSim.FreePrecessionMatrixType
FreePrecessionMatrix(E1, E2, θ)
FreePrecessionMatrix{T}()
FreePrecessionMatrix()

Create a mutable FreePrecessionMatrix object encoding the effects of relaxation and off-resonance precession.

Properties

• E1::Real: T1 relaxation
• E2::Real: T2 relaxation
• θ::Real: Angle of off-resonance precession (rad)

Examples

julia> T1 = 1000; T2 = 100; Δω = π/600; t = 100;

julia> F = FreePrecessionMatrix(exp(-t / T1), exp(-t / T2) * cos(Δω * t), exp(-t / T2) * sin(Δω * t))
FreePrecessionMatrix{Float64}:
E1 = 0.9048374180359595
E2 = 0.31859294158449203

julia> Matrix(F)
3×3 Matrix{Float64}:
0.313219  0.058272  0.0
-0.058272  0.313219  0.0
0.0       0.0       0.904837
BlochSim.GradientType
Gradient(x, y, z)

Create a Gradient object representing x, y, and z B0 gradients. Units are G/cm.

Properties

• x::Real: x gradient
• y::Real: y gradient
• z::Real: z gradient
BlochSim.GradientSpoilingType
GradientSpoiling(grad, Tg) <: AbstractSpoiling
GradientSpoiling(gx, gy, gz, Tg)

Represents gradient spoiling, e.g., applying a gradient grad = Gradient(gx, gy, gz) for time Tg (ms). grad can be a Gradient (or gx, gy, and gz can be scalars), representing a constant gradient, or grad can be a collection of Gradients (or gx, gy, and gz can be collections of values), representing a gradient waveform with a constant time step.

BlochSim.IdealSpoilingType
IdealSpoiling() <: AbstractSpoiling

Represents ideal spoiling, i.e., setting the transverse (x and y) components of a spin's magnetization to 0.

BlochSim.IdealSpoilingMatrixType
idealspoiling = IdealSpoilingMatrix()

Matrix representing ideal spoiling. Multiplying by a Magnetization or MagnetizationMC has the effect of setting the x and y components to 0.

Examples

julia> idealspoiling * MagnetizationMC((1, 1, 1), (2, 2, 2))
2-compartment Magnetization vector with eltype Int64:
Compartment 1:
Mx = 0
My = 0
Mz = 1
Compartment 2:
Mx = 0
My = 0
Mz = 2
BlochSim.InstantaneousRFType
InstantaneousRF(α, θ = 0) <: AbstractRF

Represents an idealized instantaneous RF pulse with flip angle α and phase θ.

BlochSim.MESEBlochSimType
mese! = MESEBlochSim(TR, TE, nechoes, [rfex, rfref, rephaser, crusher, spoiling])
mese!(spin, [workspace])

Simulate a multi-echo spin echo (MESE) scan on spin, overwriting the spin's magnetization vector. Returns a Vector with the magnetization vectors at each echo.

Arguments

• TR::Real: Repetition time (ms)
• TE::Real: First echo time, and echo spacing (ms); the first echo time is measured from the middle of the excitation pulse
• nechoes::Integer: Number of echoes to readout
• rfex::AbstractRF = InstantaneousRF(π/2): Excitation RF pulse
• rfref::AbstractRF = InstantaneousRF(π, -π/2): Refocussing RF pulse
• rephaser::Union{<:GradientSpoiling,Nothing} = nothing: Slice-select excitation rephasing gradient
• crusher::Union{<:GradientSpoiling,Nothing} = nothing: Crusher gradient (placed on either side of each refocussing pulse)
• spoiling::Union{IdealSpoiling,<:GradientSpoiling,Nothing} = IdealSpoiling(): Type of spoiling to apply

workspace isa MESEBlochSimWorkspace.

BlochSim.MagnetizationType
Magnetization(x, y, z)
Magnetization{T}()
Magnetization()

Create a mutable Magnetization object representing a 3D magnetization vector.

Properties

• x::Real: x component of magnetization vector
• y::Real: y component of magnetization vector
• z::Real: z component of magnetization vector

Examples

julia> M = Magnetization()
Magnetization vector with eltype Float64:
Mx = 0.0
My = 0.0
Mz = 0.0

julia> M.x = 1; M.y = 2; M.z = 3; M
Magnetization vector with eltype Float64:
Mx = 1.0
My = 2.0
Mz = 3.0
BlochSim.MagnetizationMCType
MagnetizationMC((x1, y1, z1), ...)
MagnetizationMC{T}(N)
MagnetizationMC(N)

Create a MagnetizationMC object representing N 3D magnetization vectors in the same physical location.

One can access the ith component magnetization vector by indexing the MagnetizationMC object. Furthermore, iterating the MagnetizationMC object iterates through each of the component magnetization vectors.

Properties

• M::NTuple{N,Magnetization{<:Real}}: List of component magnetization vectors

Examples

julia> M = MagnetizationMC((1, 2, 3), (4, 5, 6))
2-compartment Magnetization vector with eltype Int64:
Compartment 1:
Mx = 1
My = 2
Mz = 3
Compartment 2:
Mx = 4
My = 5
Mz = 6

julia> M[2]
Magnetization vector with eltype Int64:
Mx = 4
My = 5
Mz = 6

julia> foreach(m -> (m.x = 0; m.y = 1; m.z = 2), M)

julia> M
2-compartment Magnetization vector with eltype Int64:
Compartment 1:
Mx = 0
My = 1
Mz = 2
Compartment 2:
Mx = 0
My = 1
Mz = 2
BlochSim.PositionType
Position(x, y, z)

Create a mutable Position object representing a 3D location. Units are cm.

Properties

• x::Real: x position
• y::Real: y position
• z::Real: z position
BlochSim.RFType
RF(waveform, Δt, [Δθ], [grad]) <: AbstractRF

Represents an RF pulse with the given (possibly complex-valued) waveform (G) and time step Δt (ms). Δθ is additional phase added to the waveform (defaults to 0), and grad is the B0 gradient that is turned on during the RF pulse (defaults to Gradient(0, 0, 0), i.e., turned off).

Properties

• α::Vector{<:Real}: Instantaneous flip angles (rad) at each time point; computed from the magnitude of waveform
• θ::Vector{<:Real}: Instantaneous phase (rad) at each time point; computed from the phase of waveform
• Δt::Real: Time step (ms)
• Δθ_initial::Real: Phase added to θ before any phase-cycling increment has been applied
• Δθ::Ref{<:Real}: Phase to be added to θ; can be updated to simulate phase-cycling/RF spoiling
• grad: Gradient applied during the RF pulse
• ::Gradient: Constant gradient
• ::Vector{<:Gradient}: Time-varying gradient
BlochSim.RFSpoilingType
RFSpoiling(Δθ = 117°) <: AbstractSpoiling

Represents RF spoiling, i.e., quadratically incrementing the phase of the RF pulses from TR to TR.

BlochSim.SPGRBlochSimType
spgr! = SPGRBlochSim(TR, TE, rf, [spoiling], [nTR], [save_transients])
spgr!(spin, [workspace])

Simulate a spoiled gradient-recalled echo (SPGR) scan on spin, overwriting the spin's magnetization vector. The resultant magnetization is stored in spin.M. If nTR > 0 and save_transients === true, then spgr!(...) returns a Vector with the magnetization vectors at the echo time for each of the nTR simulated TRs.

Arguments

• TR::Real: Repetition time (ms)
• TE::Real: Echo time (ms)
• rf:
• ::Real: Excitation flip angle (rad) (assumes an instantaneous RF pulse)
• ::AbstractRF: Excitation RF pulse
• spoiling::AbstractSpoiling = IdealSpoiling(): Type of spoiling to apply
• nTR::Val = Val(0): Number of TRs to simulate; Val(0) indicates to simulate a steady-state scan
• save_transients::Val = Val(false): Whether or not to return the magnetization vectors at the TE for each of the nTR simulated TRs; does nothing if nTR == 0

workspace isa SPGRBlochSimWorkspace.

BlochSim.SpinType
Spin([M], M0, T1, T2, Δf, [pos]) <: AbstractSpin

Create an object that represents a single spin.

Properties

• M::Magnetization = Magnetization(0, 0, M0): Magnetization vector
• M0::Real: Equilibrium magnetization
• T1::Real: Spin-lattice recovery time constant (ms)
• T2::Real: Spin-spin recovery time constant (ms)
• Δf::Real: Off-resonance frequency (Hz)
• pos::Position = Position(0, 0, 0): Spatial location (cm)
• N::Int = 1: Number of compartments

Examples

julia> Spin(1, 1000, 100, 0, Position(1, 2, 3))
Spin{Float64}:
M = Magnetization(0.0, 0.0, 1.0)
M0 = 1.0
T1 = 1000.0 ms
T2 = 100.0 ms
Δf = 0.0 Hz
pos = Position(1.0, 2.0, 3.0) cm

julia> Spin(Magnetization(1, 2, 3), 1, 1000, 100, 0)
Spin{Float64}:
M = Magnetization(1.0, 2.0, 3.0)
M0 = 1.0
T1 = 1000.0 ms
T2 = 100.0 ms
Δf = 0.0 Hz
pos = Position(0.0, 0.0, 0.0) cm
BlochSim.SpinMCType
SpinMC([M], M0, frac, T1, T2, Δf, τ, [pos]) <: AbstractSpin

Create an object that represents a single spin with multiple compartments.

Properties

• M::MagnetizationMC = Meq: Magnetization vector
• Meq::MagnetizationMC = MagnetizationMC((0, 0, frac[1] * M0), ...): Equilibrium magnetization vector
• M0::Real: Equilibrium magnetization
• frac::Tuple{<:Real}: Water fraction of each compartment
• T1::Tuple{<:Real}: Spin-lattice recovery time constants (ms)
• T2::Tuple{<:Real}: Spin-spin recovery time constants (ms)
• Δf::Tuple{<:Real}: Off-resonance frequencies (Hz)
• r::Tuple{Tuple{<:Real}}: Exchange rates (1/ms); r[i][j] is the exchange rate from compartment i to compartment j
• pos::Position = Position(0, 0, 0): Spatial location (cm)
• N::Int = length(frac): Number of compartments

Note

The SpinMC constructor takes τ (inverse exchange rate, or residence time) as input, not r. Furthermore, τ is given as a Tuple with N^2 - N elements, arranged like (τ12, τ13, ..., τ1N, τ21, τ23, ..., τ2N, ...).

Examples

julia> SpinMC(1, (0.2, 0.8), (400, 1000), (20, 100), (15, 0), (100, 25))
SpinMC{Float64,2}:
M = MagnetizationMC((0.0, 0.0, 0.2), (0.0, 0.0, 0.8))
M0 = 1.0
frac = (0.2, 0.8)
T1 = (400.0, 1000.0) ms
T2 = (20.0, 100.0) ms
Δf = (15.0, 0.0) Hz
r = ((0.0, 0.01), (0.04, 0.0)) 1/ms
pos = Position(0.0, 0.0, 0.0) cm
BlochSim.add!Method
add!(C, A, B)

Compute A + B, storing the result in C.

BlochSim.add!Method
add!(A, B)

Compute A + B, storing the result in A.

BlochSim.applydynamics!Method
applydynamics!(spin, BtoM, A, [B])

Apply dynamics to the given spin, overwriting the spin's magnetization vector. BtoM is used to store intermediate results (and is thus overwritten).

Examples

julia> s = Spin(1, 1000, 100, 3.75); s.M
Magnetization vector with eltype Float64:
Mx = 0.0
My = 0.0
Mz = 1.0

julia> BtoM = Magnetization();

julia> (A,) = excite(s, InstantaneousRF(π/2)); applydynamics!(s, BtoM, A); s.M
Magnetization vector with eltype Float64:
Mx = 1.0
My = 0.0
Mz = 6.123233995736766e-17

julia> (A, B) = freeprecess(s, 100); applydynamics!(s, BtoM, A, B); s.M
Magnetization vector with eltype Float64:
Mx = -0.2601300475114444
My = -0.2601300475114445
Mz = 0.09516258196404054
BlochSim.combine!Method
combine!(A, B, A1, B1, A2, B2)
combine!(A, A1, A2)

Combine the matrices and vectors that describe the dynamics of a spin into one matrix and one vector, overwriting A and B. The dynamics described by A1 and B1 apply first, then those described by A2 and B2. In other words, A = A2 * A1 and B = A2 * B1 + B2.

Examples

julia> s = Spin(1, 1000, 100, 3.75);

julia> A = BlochMatrix(); B = Magnetization();

julia> (A1, B1) = excite(s, InstantaneousRF(π/2));

julia> (A2, B2) = freeprecess(s, 100);

julia> combine!(A, B, A1, B1, A2, B2); A * s.M + B
Magnetization vector with eltype Float64:
Mx = -0.2601300475114444
My = -0.2601300475114445
Mz = 0.09516258196404054
BlochSim.div!Method
div!(A, a)

Compute A / a, storing the result in A.

BlochSim.exciteFunction
excite(spin, rf::InstantaneousRF, [nothing])
excite(spin, rf::RF, [workspace])

Simulate excitation for the given spin. Returns (A, B) such that A * M + B applies excitation to the magnetization M. If isnothing(B) (as is the case for InstantaneousRFs), then A * M applies excitation to M.

For RF objects, workspace isa ExcitationWorkspace. For SpinMC objects, use workspace = ExcitationWorkspace(spin, nothing) to use an approximate matrix exponential to solve the Bloch-McConnell equation.

For an in-place version, see excite!.

Examples

julia> s = Spin(1, 1000, 100, 3.75); s.M
Magnetization vector with eltype Float64:
Mx = 0.0
My = 0.0
Mz = 1.0

julia> (A,) = excite(s, InstantaneousRF(π/2, π/4)); A * s.M
Magnetization vector with eltype Float64:
Mx = 0.7071067811865476
My = -0.7071067811865475
Mz = 6.123233995736766e-17
BlochSim.excite!Method
excite!(spin, ...)

Apply excitation to the given spin, overwriting the spin's magnetization vector.

BlochSim.excite!Method
excite!(A, [nothing], spin, rf::InstantaneousRF, [nothing])
excite!(A, B, spin, rf::RF, [workspace])

Simulate excitation, overwriting A and B (in-place version of excite).

BlochSim.expm!Method
expm!(expA, A, [workspace])

Compute the matrix exponential of A, storing it in expA.

workspace isa MatrixExponentialWorkspace.

BlochSim.freeprecessFunction
freeprecess(spin, t, [nothing])
freeprecess(spinmc, t, [workspace])
freeprecess(spinmc, t, grad, [workspace])

Simulate free-precession for the given spin for time t ms, optionally in the presence of a B0 gradient. Returns (A, B) such that A * M + B applies free-precession to the magnetization M.

For SpinMC objects, workspace isa BlochMcConnellWorkspace. Pass in nothing instead to use an approximate matrix exponential to solve the Bloch-McConnell equation.

For an in-place version, see freeprecess!.

Examples

julia> s = Spin(Magnetization(1, 0, 0), 1, 1000, 100, 3.75)
Spin{Float64}:
M = Magnetization(1.0, 0.0, 0.0)
M0 = 1.0
T1 = 1000.0 ms
T2 = 100.0 ms
Δf = 3.75 Hz
pos = Position(0.0, 0.0, 0.0) cm

julia> (A, B) = freeprecess(s, 100); A * s.M + B
Magnetization vector with eltype Float64:
Mx = -0.2601300475114444
My = -0.2601300475114445
Mz = 0.09516258196404048

julia> s = Spin(Magnetization(1, 0, 0), 1, 1000, 100, 0, Position(0, 0, 3.75))
Spin{Float64}:
M = Magnetization(1.0, 0.0, 0.0)
M0 = 1.0
T1 = 1000.0 ms
T2 = 100.0 ms
Δf = 0.0 Hz
pos = Position(0.0, 0.0, 3.75) cm

julia> (A, B) = freeprecess(s, 100, Gradient(0, 0, 1/GAMBAR)); A * s.M + B
Magnetization vector with eltype Float64:
Mx = -0.2601300475114444
My = -0.2601300475114445
Mz = 0.09516258196404048
BlochSim.freeprecess!Method
freeprecess!(A, B, t, M0, T1, T2, Δf)
freeprecess!(A, B, spin, t, [nothing])
freeprecess!(A, B, spinmc, t, [workspace])
freeprecess!(A, B, spin, t, grad, [nothing])
freeprecess!(A, B, spinmc, t, grad, [workspace])

Simulate free-precession, overwriting A and B (in-place version of freeprecess).

BlochSim.freeprecess!Method
freeprecess!(spin, ...)

Apply free-precession to the given spin, overwriting the spin's magnetization vector.

BlochSim.freeprecessMethod
freeprecess(t, M0, T1, T2, Δf)

Simulate free-precession, i.e., relaxation and off-resonance precession. Returns (A, B) such that A * M + B applies free-precession to the magnetization M.

For an in-place version, see freeprecess!.

Arguments

• t::Real: Duration of free-precession (ms)
• M0::Real: Equilibrium magnetization
• T1::Real: Spin-lattice recovery time constant (ms)
• T2::Real: Spin-spin recovery time constant (ms)
• Δf::Real: Off-resonance frequency (Hz)

Examples

julia> (A, B) = freeprecess(100, 1, 1000, 100, 3.75); A * Magnetization(1, 0, 0) + B
Magnetization vector with eltype Float64:
Mx = -0.2601300475114444
My = -0.2601300475114445
Mz = 0.09516258196404048
BlochSim.gradient_frequencyMethod
gradient_frequency(grad, pos)

Compute the off-resonance frequency in Hz induced by the given B0 gradient grad at position pos.

BlochSim.muladd!Method
muladd!(C, A, B)

Compute A * B + C, storing the result in C.

BlochSim.rotatethetaFunction
rotatetheta(α::Real = π/2, θ::Real = 0)

Return a 3 × 3 BlochMatrix for (left-handed) flip angle α about an axis in the x-y plane that makes (left-handed) angle θ with the negative y-axis.

For an in-place version, see rotatetheta!.

Examples

julia> BlochSim.rotatetheta() * Magnetization(0, 0, 1) # Rotate towards positive x-axis
Magnetization vector with eltype Float64:
Mx = 1.0
My = 0.0
Mz = 6.123233995736766e-17

julia> BlochSim.rotatetheta(π/2, π/2) * Magnetization(0, 0, 1) # Rotate towards negative y-axis
Magnetization vector with eltype Float64:
Mx = 6.123233995736766e-17
My = -1.0
Mz = 6.123233995736766e-17
BlochSim.rotatetheta!Method
rotatetheta!(A, α, θ)

Simulate left-handed rotation by angle α about an axis in the x-y plane that makes left-handed angle θ with the negative y-axis, overwriting A.

This function is an in-place version of rotatetheta.

BlochSim.signalMethod
signal(spin)
signal(M)

Return the signal detected from the given spin or magnetization vector.

Examples

julia> signal(Spin(Magnetization(1, 2, 3), 1, 1000, 100, 0))
1.0 + 2.0im

julia> signal(MagnetizationMC((1, 2, 3), (1, 1, 1)))
2 + 3im
BlochSim.spoilFunction
spoil(spin, spoiling, [nothing])
spoil(spinmc, spoiling, [workspace])

Simulate gradient or ideal spoiling for the given spin. Returns (A, B), such that A * M + B applies spoiling to the magnetization M. If B is nothing (as is the case for IdealSpoiling), then A * M applies spoiling, and if both A and B are nothing (as is the case for RFSpoiling) then there is no spoiling.

For SpinMC objects and for GradientSpoiling and RFandGradientSpoiling, workspace isa BlochMcConnellWorkspace. Pass in nothing instead to use an approximate matrix exponential to solve the Bloch-McConnell equation.

Note

This function only simulates gradient or ideal spoiling, not RF spoiling. RF spoiling must be implemented by updating the phase of the RF pulse(s) in your sequence from TR to TR.

For an in-place version, see spoil!. 

BlochSim.spoil!Function
spoil!(A, B, spin, spoiling, [nothing])
spoil!(A, B, spinmc, spoiling, [workspace])

Simulate gradient or ideal spoiling, overwriting A and B (in-place version of spoil`).