Unitful.jl to create Sequences, Phantoms, and Scanner types.

[cncastillo]

This would greatly improve the readability of the "Pulse Sequence" without needing to guess the units. Although we haven't tested it yet, it should not be difficult.

[cncastillo]

I think this is related to the goal of adding automatic differentiation.

[cncastillo]

I already started doing this for the Grads in #321

using Parameters
abstract type MRISequenceEvent end
NumOrVect = Union{T, Vector{T}} where {T<:Number}

@with_kw_noshow mutable struct Grad{GradAmp<:NumOrVect, TimeDur<:NumOrVect} <: MRISequenceEvent
    # Grad properties
    A    ::GradAmp = zero(eltype(GradAmp))
    T    ::TimeDur = zero(eltype(TimeDur))
    rise ::eltype(TimeDur) = zero(eltype(T))
    fall ::eltype(TimeDur) = rise
    delay::eltype(TimeDur) = zero(eltype(T))
    first::eltype(GradAmp) = zero(eltype(A))
    last ::eltype(GradAmp) = first
    # Checking validity of Grad
    @assert length(A) >= length(T)
    @assert all(T .>= zero(eltype(T)))
    @assert rise   >= zero(eltype(T))
    @assert fall   >= zero(eltype(T))
    @assert delay  >= zero(eltype(T))
end

# Gradient type aliases for dispatching
TrapezoidalGrad      = Grad{G, T} where {G, T}
UniformlySampledGrad = Grad{Vector{G}, T} where {G, T}
TimeShapedGrad       = Grad{Vector{G}, Vector{T}} where {G, T}
ValidGrad = Union{TrapezoidalGrad{G, T}, UniformlySampledGrad{G, T}, TimeShapedGrad{G, T}} where {G, T}

# MRISequenceEvent functions
# Event duration
dur(g::Grad) = g.delay + g.rise + sum(g.T) + g.fall
dur(G::Vector{Grad{G, T}}) where {G, T} = max(dur.(G)...)
# Event amplitude
ampl(g::TrapezoidalGrad{G,T})      where {G, T} = vcat(zero(G), g.first, g.A, g.A, g.last)
ampl(g::TimeShapedGrad{G,T})       where {G, T} = vcat(zero(G), g.first, g.A, g.last)
ampl(g::UniformlySampledGrad{G,T}) where {G, T} = vcat(zero(G), g.first, g.A, g.last)
# Event time
time(g::TrapezoidalGrad{G,T})      where {G, T} = cumsum(vcat(zero(T), g.delay, g.rise, g.T, g.fall))
time(g::TimeShapedGrad{G,T})       where {G, T} = cumsum(vcat(zero(T), g.delay, g.rise, g.T, g.fall))
time(g::UniformlySampledGrad{G,T}) where {G, T} = begin
    NA = length(g.A) - 1
    ΔT = repeat([g.T / NA], NA) 
    return cumsum(vcat(zero(T), g.delay, g.rise, ΔT, g.fall))
end

[cncastillo]

For Unitful support

using Unitful
using Unitful: 𝐈, 𝐌, 𝐓, 𝐋
using Unitful: T, m, mT, s, ms, Hz, MHz
Base.show(io::IO, g::ValidGrad{QA, QT}) where {QA <: Quantity, QT <: Quantity} = begin
    printfield(g, f) = begin
        val = getfield(g, f)
        if eltype(val) <: QA && !all(abs.(val) .== zero(QA)) 
            return "$f=$(ustrip(val)) $(unit(QA))"
        elseif eltype(val) <: QT && !all(abs.(val) .== zero(QT)) 
            return "$f=$(ustrip(val)) $(unit(QT))"
        else
            return ""
        end
    end
    text = [printfield(g, f) for f in fieldnames(Grad)]
    text = join(text[text.!=""], ", ")
    print(io, "Grad($text)")
end

γ = 42.5774688MHz/T

# Defining gradient units
MagneticFieldGradient = Union{Quantity{T, 𝐈^-1*𝐌*𝐓^-2*𝐋^-1, U}, Level{L, S, Quantity{T, 𝐈^-1*𝐌*𝐓^-2*𝐋^-1, U}} where {L, S}} where {T, U}

# Promotion rules
Unitful.promote_to_derived()
Unitful.@derived_dimension BFieldGradient 𝐈^-1*𝐌*𝐓^-2*𝐋^-1 true
Unitful.@derived_dimension Frequency 𝐓^-1 true
Unitful.promote_unit(::S, ::T) where {S<:BFieldGradientFreeUnits, T<:BFieldGradientFreeUnits} = Unitful.T / Unitful.m
Unitful.promote_unit(::S, ::T) where {S<:FrequencyFreeUnits,      T<:FrequencyFreeUnits}      = Unitful.Hz
end