Skip to content

Commit

Permalink
Merge pull request #23 from StevenWhitaker/mese
Browse files Browse the repository at this point in the history 
Add rephasing, crushing, and spoiling to `MESEBlochSim`.
  • Loading branch information
@StevenWhitaker
StevenWhitaker authored Aug 13, 2022
2 parents e25ee10 + 492553a commit 68eba01
Show file tree
Hide file tree
Showing 11 changed files with 485 additions and 147 deletions.
160 changes: 114 additions & 46 deletions src/blochmatrix.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,21 @@ function Base.show(io::IO, ::MIME"text/plain", A::BlochMatrix{T}) where {T}

end

function make_identity!(A::BlochMatrix{T}) where {T}

A.a11 = one(T)
A.a21 = zero(T)
A.a31 = zero(T)
A.a12 = zero(T)
A.a22 = one(T)
A.a32 = zero(T)
A.a13 = zero(T)
A.a23 = zero(T)
A.a33 = one(T)
return nothing

end

function Base.fill!(A::BlochMatrix, v)

A.a11 = v
Expand Down Expand Up @@ -125,7 +140,7 @@ function Base.Matrix(A::BlochDynamicsMatrix{T}) where {T}
end

"""
FreePrecessionMatrix(E1, E2cosθ, E2sinθ)
FreePrecessionMatrix(E1, E2, θ)
FreePrecessionMatrix{T}()
FreePrecessionMatrix()
Expand All @@ -134,8 +149,8 @@ relaxation and off-resonance precession.
# Properties
- `E1::Real`: T1 relaxation
- `E2cosθ::Real`: T2 relaxation and off-resonance precession
- `E2sinθ::Real`: T2 relaxation and off-resonance precession
- `E2::Real`: T2 relaxation
- `θ::Real`: Angle of off-resonance precession (rad)
# Examples
```jldoctest
Expand All @@ -144,8 +159,8 @@ 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
E2cosθ = 0.31859294158449203
E2sinθ = 0.18393972058572114
E2 = 0.36787944117144233
θ = 0.5235987755982988
julia> Matrix(F)
3×3 Matrix{Float64}:
Expand All @@ -156,31 +171,44 @@ julia> Matrix(F)
"""
mutable struct FreePrecessionMatrix{T<:Real} <: AbstractBlochMatrix{T}
E1::T
E2cosθ::T
E2sinθ::T
E2::T
θ::T
end

FreePrecessionMatrix{T}() where {T} = FreePrecessionMatrix(zero(T), zero(T), zero(T))
FreePrecessionMatrix{T}() where {T} = FreePrecessionMatrix(one(T), one(T), zero(T))
FreePrecessionMatrix() = FreePrecessionMatrix{Float64}()
FreePrecessionMatrix(E1, E2cosθ, E2sinθ) = FreePrecessionMatrix(promote(E1, E2cosθ, E2sinθ)...)
FreePrecessionMatrix(E1, E2, θ) = FreePrecessionMatrix(promote(E1, E2, θ)...)

function Base.show(io::IO, ::MIME"text/plain", A::FreePrecessionMatrix{T}) where {T}

print(io, "FreePrecessionMatrix{$T}:")
print(io, "\n E1 = ", A.E1)
print(io, "\n E2cosθ = ", A.E2cosθ)
print(io, "\n E2sinθ = ", A.E2sinθ)
print(io, "\n E2 = ", A.E2)
print(io, "\n θ = ", A.θ, " rad")

end

function Base.Matrix(A::FreePrecessionMatrix{T}) where {T}
function make_identity!(A::FreePrecessionMatrix{T}) where {T}

A.E1 = one(T)
A.E2 = one(T)
A.θ = zero(T)
return nothing

end

function Base.Matrix(A::FreePrecessionMatrix)

(sinθ, cosθ) = sincos(A.θ)
E2cosθ = A.E2 * cosθ
E2sinθ = A.E2 * sinθ
T = promote_type(typeof(E2cosθ), typeof(E2sinθ))
mat = Matrix{T}(undef, 3, 3)
mat[1,1] = A.E2cosθ
mat[2,1] = -A.E2sinθ
mat[1,1] = E2cosθ
mat[2,1] = -E2sinθ
mat[3,1] = zero(T)
mat[1,2] = A.E2sinθ
mat[2,2] = A.E2cosθ
mat[1,2] = E2sinθ
mat[2,2] = E2cosθ
mat[3,2] = zero(T)
mat[1,3] = zero(T)
mat[2,3] = zero(T)
Expand All @@ -192,11 +220,12 @@ end

function Base.copyto!(dst::BlochMatrix{T}, src::FreePrecessionMatrix) where {T}

dst.a11 = src.E2cosθ
dst.a21 = -src.E2sinθ
(sinθ, cosθ) = sincos(src.θ)
dst.a11 = src.E2 * cosθ
dst.a21 = -src.E2 * sinθ
dst.a31 = zero(T)
dst.a12 = src.E2sinθ
dst.a22 = src.E2cosθ
dst.a12 = src.E2 * sinθ
dst.a22 = src.E2 * cosθ
dst.a32 = zero(T)
dst.a13 = zero(T)
dst.a23 = zero(T)
Expand All @@ -205,6 +234,15 @@ function Base.copyto!(dst::BlochMatrix{T}, src::FreePrecessionMatrix) where {T}

end

function Base.copyto!(dst::FreePrecessionMatrix, src::FreePrecessionMatrix)

dst.E1 = src.E1
dst.E2 = src.E2
dst.θ = src.θ
return nothing

end

"""
ExchangeDynamicsMatrix(r)
ExchangeDynamicsMatrix{T}()
Expand Down Expand Up @@ -404,6 +442,19 @@ end

getblock(A::BlochMcConnellMatrix, i, j) = A.A[i][j]

function make_identity!(A::BlochMcConnellMatrix{T,N}) where {T,N}

for j = 1:N, i = 1:N
if i == j
make_identity!(getblock(A, i, j))
else
fill!(getblock(A, i, j), 0)
end
end
return nothing

end

function Base.fill!(A::BlochMcConnellMatrix{T,N}, v) where {T,N}

for i = 1:N, j = 1:N
Expand Down Expand Up @@ -622,9 +673,10 @@ end

function Base.:*(A::FreePrecessionMatrix, M::Magnetization)

(sinθ, cosθ) = sincos(A.θ)
return Magnetization(
A.E2cosθ * M.x + A.E2sinθ * M.y,
A.E2cosθ * M.y - A.E2sinθ * M.x,
A.E2 * cosθ * M.x + A.E2 * sinθ * M.y,
A.E2 * cosθ * M.y - A.E2 * sinθ * M.x,
A.E1 * M.z
)

Expand Down Expand Up @@ -666,13 +718,14 @@ end

function Base.:*(A::BlochMatrix, B::FreePrecessionMatrix)

(sinθ, cosθ) = sincos(B.θ)
return BlochMatrix(
A.a11 * B.E2cosθ - A.a12 * B.E2sinθ,
A.a21 * B.E2cosθ - A.a22 * B.E2sinθ,
A.a31 * B.E2cosθ - A.a32 * B.E2sinθ,
A.a11 * B.E2sinθ + A.a12 * B.E2cosθ,
A.a21 * B.E2sinθ + A.a22 * B.E2cosθ,
A.a31 * B.E2sinθ + A.a32 * B.E2cosθ,
A.a11 * B.E2 * cosθ - A.a12 * B.E2 * sinθ,
A.a21 * B.E2 * cosθ - A.a22 * B.E2 * sinθ,
A.a31 * B.E2 * cosθ - A.a32 * B.E2 * sinθ,
A.a11 * B.E2 * sinθ + A.a12 * B.E2 * cosθ,
A.a21 * B.E2 * sinθ + A.a22 * B.E2 * cosθ,
A.a31 * B.E2 * sinθ + A.a32 * B.E2 * cosθ,
A.a13 * B.E1,
A.a23 * B.E1,
A.a33 * B.E1
Expand Down Expand Up @@ -934,8 +987,9 @@ end

function LinearAlgebra.mul!(M2::Magnetization, A::FreePrecessionMatrix, M1::Magnetization)

M2.x = A.E2cosθ * M1.x + A.E2sinθ * M1.y
M2.y = A.E2cosθ * M1.y - A.E2sinθ * M1.x
(sinθ, cosθ) = sincos(A.θ)
M2.x = A.E2 * cosθ * M1.x + A.E2 * sinθ * M1.y
M2.y = A.E2 * cosθ * M1.y - A.E2 * sinθ * M1.x
M2.z = A.E1 * M1.z
return nothing

Expand Down Expand Up @@ -1031,12 +1085,13 @@ end

function LinearAlgebra.mul!(C::BlochMatrix, A::BlochMatrix, B::FreePrecessionMatrix)

C.a11 = A.a11 * B.E2cosθ - A.a12 * B.E2sinθ
C.a21 = A.a21 * B.E2cosθ - A.a22 * B.E2sinθ
C.a31 = A.a31 * B.E2cosθ - A.a32 * B.E2sinθ
C.a12 = A.a11 * B.E2sinθ + A.a12 * B.E2cosθ
C.a22 = A.a21 * B.E2sinθ + A.a22 * B.E2cosθ
C.a32 = A.a31 * B.E2sinθ + A.a32 * B.E2cosθ
(sinθ, cosθ) = sincos(B.θ)
C.a11 = A.a11 * B.E2 * cosθ - A.a12 * B.E2 * sinθ
C.a21 = A.a21 * B.E2 * cosθ - A.a22 * B.E2 * sinθ
C.a31 = A.a31 * B.E2 * cosθ - A.a32 * B.E2 * sinθ
C.a12 = A.a11 * B.E2 * sinθ + A.a12 * B.E2 * cosθ
C.a22 = A.a21 * B.E2 * sinθ + A.a22 * B.E2 * cosθ
C.a32 = A.a31 * B.E2 * sinθ + A.a32 * B.E2 * cosθ
C.a13 = A.a13 * B.E1
C.a23 = A.a23 * B.E1
C.a33 = A.a33 * B.E1
Expand All @@ -1046,14 +1101,15 @@ end

function LinearAlgebra.mul!(C::BlochMatrix, A::FreePrecessionMatrix, B::BlochMatrix)

C.a11 = A.E2cosθ * B.a11 + A.E2sinθ * B.a21
C.a21 = A.E2cosθ * B.a21 - A.E2sinθ * B.a11
(sinθ, cosθ) = sincos(A.θ)
C.a11 = A.E2 * cosθ * B.a11 + A.E2 * sinθ * B.a21
C.a21 = A.E2 * cosθ * B.a21 - A.E2 * sinθ * B.a11
C.a31 = A.E1 * B.a31
C.a12 = A.E2cosθ * B.a12 + A.E2sinθ * B.a22
C.a22 = A.E2cosθ * B.a22 - A.E2sinθ * B.a12
C.a12 = A.E2 * cosθ * B.a12 + A.E2 * sinθ * B.a22
C.a22 = A.E2 * cosθ * B.a22 - A.E2 * sinθ * B.a12
C.a32 = A.E1 * B.a32
C.a13 = A.E2cosθ * B.a13 + A.E2sinθ * B.a23
C.a23 = A.E2cosθ * B.a23 - A.E2sinθ * B.a13
C.a13 = A.E2 * cosθ * B.a13 + A.E2 * sinθ * B.a23
C.a23 = A.E2 * cosθ * B.a23 - A.E2 * sinθ * B.a13
C.a33 = A.E1 * B.a33
return nothing

Expand Down Expand Up @@ -1112,8 +1168,10 @@ end

function LinearAlgebra.mul!(C::BlochMatrix{T}, A::FreePrecessionMatrix, B::FreePrecessionMatrix) where {T}

C.a11 = A.E2cosθ * B.E2cosθ - A.E2sinθ * B.E2sinθ
C.a12 = A.E2cosθ * B.E2sinθ + A.E2sinθ * B.E2cosθ
(Asinθ, Acosθ) = sincos(A.θ)
(Bsinθ, Bcosθ) = sincos(B.θ)
C.a11 = A.E2 * Acosθ * B.E2 * Bcosθ - A.E2 * Asinθ * B.E2 * Bsinθ
C.a12 = A.E2 * Acosθ * B.E2 * Bsinθ + A.E2 * Asinθ * B.E2 * Bcosθ
C.a13 = zero(T)
C.a21 = -C.a12
C.a22 = C.a11
Expand All @@ -1125,6 +1183,15 @@ function LinearAlgebra.mul!(C::BlochMatrix{T}, A::FreePrecessionMatrix, B::FreeP

end

function LinearAlgebra.mul!(C::FreePrecessionMatrix, A::FreePrecessionMatrix, B::FreePrecessionMatrix)

C.E1 = A.E1 * B.E1
C.E2 = A.E2 * B.E2
C.θ = A.θ + B.θ
return nothing

end

function LinearAlgebra.mul!(C::BlochMatrix, A::FreePrecessionMatrix, B::ExcitationMatrix)

mul!(C, A, B.A)
Expand Down Expand Up @@ -1274,8 +1341,9 @@ end

function muladd!(M2::Magnetization, A::FreePrecessionMatrix, M1::Magnetization)

M2.x += A.E2cosθ * M1.x + A.E2sinθ * M1.y
M2.y += A.E2cosθ * M1.y - A.E2sinθ * M1.x
(sinθ, cosθ) = sincos(A.θ)
M2.x += A.E2 * cosθ * M1.x + A.E2 * sinθ * M1.y
M2.y += A.E2 * cosθ * M1.y - A.E2 * sinθ * M1.x
M2.z += A.E1 * M1.z
return nothing

Expand Down
87 changes: 80 additions & 7 deletions src/freeprecess.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,13 +55,9 @@ Simulate free-precession, overwriting `A` and `B` (in-place version of
"""
function freeprecess!(A, B, t, M0, T1, T2, Δf)

E2 = exp(-t / T2)
θ = 2π * Δf * t / 1000
(s, c) = sincos(θ)

A.E1 = exp(-t / T1)
A.E2cosθ = E2 * c
A.E2sinθ = E2 * s
A.E2 = exp(-t / T2)
A.θ = 2π * Δf * t / 1000

B.x = 0
B.y = 0
Expand Down Expand Up @@ -129,7 +125,7 @@ function freeprecess(spin::Spin, t, ::Nothing = nothing)

end

function freeprecess(spin::SpinMC{T,N}, t, workspace = BlochMcConnellWorkspace(spin)) where {T,N}
function freeprecess(spin::SpinMC{T,N}, t, workspace::Union{Nothing,<:BlochMcConnellWorkspace} = BlochMcConnellWorkspace(spin)) where {T,N}

A = BlochMcConnellMatrix{T}(N)
B = MagnetizationMC{T}(N)
Expand Down Expand Up @@ -212,6 +208,83 @@ function freeprecess!(

end

struct FreePrecessionWorkspace{T1,T2,T3}
Af::T1
Bf::T2
tmpA::T1
tmpB::T2
bm_workspace::T3
end

function FreePrecessionWorkspace(
spin::AbstractSpin,
bm_workspace = spin isa Spin ? nothing : BlochMcConnellWorkspace(spin)
)

FreePrecessionWorkspace(typeof(spin), bm_workspace)

end

function FreePrecessionWorkspace(
spin::Union{Type{Spin{T}},Type{SpinMC{T,N}}},
bm_workspace = spin <: Spin ? nothing : BlochMcConnellWorkspace(spin)
) where{T,N}

if spin <: Spin
Af = FreePrecessionMatrix{T}()
Bf = Magnetization{T}()
tmpA = FreePrecessionMatrix{T}()
tmpB = Magnetization{T}()
else
Af = BlochMcConnellMatrix{T}(N)
Bf = MagnetizationMC{T}(N)
tmpA = BlochMcConnellMatrix{T}(N)
tmpB = MagnetizationMC{T}(N)
end
FreePrecessionWorkspace(Af, Bf, tmpA, tmpB, bm_workspace)

end

function freeprecess(spin::Spin, t, grads, workspace = FreePrecessionWorkspace(spin))

A = FreePrecessionMatrix{eltype(spin)}()
B = Magnetization{eltype(spin)}()
freeprecess!(A, B, spin, t, grads, workspace)
return (A, B)

end

function freeprecess(spin::SpinMC{T,N}, t, grads, workspace = FreePrecessionWorkspace(spin)) where {T,N}

A = BlochMcConnellMatrix{T}(N)
B = MagnetizationMC{T}(N)
freeprecess!(A, B, spin, t, grads, workspace)
return (A, B)

end

function freeprecess!(
A,
B,
spin,
t,
grads, # Collection of Gradients
workspace::FreePrecessionWorkspace = FreePrecessionWorkspace(spin)
)

Δt = t / length(grads)
make_identity!(workspace.tmpA)
fill!(workspace.tmpB, 0)
for grad in grads
freeprecess!(workspace.Af, workspace.Bf, spin, Δt, grad, workspace.bm_workspace)
combine!(A, B, workspace.tmpA, workspace.tmpB, workspace.Af, workspace.Bf)
copyto!(workspace.tmpA, A)
copyto!(workspace.tmpB, B)
end
return nothing

end

# Exact matrix exponential
function expm!(expAt, workspace, spin, t, gradfreq = 0)

Expand Down
Loading

0 comments on commit 68eba01

Please sign in to comment.