diff --git a/src/gradients/changing_gradient_blocks.jl b/src/gradients/changing_gradient_blocks.jl new file mode 100644 index 0000000000000000000000000000000000000000..b42487167186a16e86eb90da167bef52f370b671 --- /dev/null +++ b/src/gradients/changing_gradient_blocks.jl @@ -0,0 +1,72 @@ +module ChangingGradientBlocks +import StaticArrays: SVector +import ...Variables: VariableType, variables +import ...BuildingBlocks: GradientBlock +import JuMP: @constraint, @variable, Model, owner_model +import ...Variables: qvec, bmat, gradient_strength, slew_rate, duration, variables, VariableType +import ...BuildingBlocks: GradientBlock, fixed, RFPulseBlock + +struct AbstractChangingGradientBlock{T<:VariableType} <: GradientBlock + gradient_strength_start :: SVector{3, T} + gradient_strength_end :: SVector{3, T} + duration :: T + rotate :: Union{Nothing, Symbol} + scale :: Union{Nothing, Symbol} +end + +""" + ChangingGradientBlock(duration, gradient_strength, rotate, scale) + +Underlying type for any linearly changing part in a gradient waveform. + +The corresponding type in a fixed gradient is [`FixedChangingGradientBlock`](@ref) + +Do not create this directly, it will be created by the gradient waveforms. +""" +const ChangingGradientBlock = AbstractChangingGradientBlock{VariableType} + +""" + FixedChangingGradientBlock(duration, gradient_strength, rotate, scale) + +Fixed equivalent of [`ChangingGradientBlock`](@ref). +""" +const FixedChangingGradientBlock = AbstractChangingGradientBlock{Float64} + +duration(cgb::AbstractChangingGradientBlock) = cgb.duration + +slew_rate(cgb::AbstractChangingGradientBlock) = (cgb.gradient_strength_end .- cgb.gradient_strength_start) ./ duration(cgb) +qvec(cgb::AbstractChangingGradientBlock) = (cgb.gradient_strength_start .+ cgb.gradient_strength_end) .* (duration(cgb)/2) + +function bmat(cgb::AbstractChangingGradientBlock, qstart) + # grad = (g1 * (duration - t) + g2 * t) / duration + # = g1 + (g2 - g1) * t / duration + # q = qstart + g1 * t + (g2 - g1) * t^2 / (2 * duration) + # \int dt (qstart + t * grad)^2 = + # qstart^2 * duration + + # qstart * g1 * duration^2 + + # qstart * (g2 - g1) * duration^2 / 3 + + # g1^2 * duration^3 / 3 + + # g1 * (g2 - g1) * duration^3 / 4 + + # (g2 - g1)^2 * duration^3 / 10 + return ( + qstart .* qstart' .* duration(cgb) .+ + duration(cgb)^2 .* qstart .* ( + 2 .* cgb.gradient_strength_start' .+ + cgb.gradient_strength_end') ./ 3 .+ + bmat(cgb) + ) +end + +function bmat(cgb::AbstractChangingGradientBlock) + diff = cgb.gradient_strength_end .- cgb.gradient_strength_start + return ( + cgb.gradient_strength_start .* cgb.gradient_strength_start' ./ 3 .+ + cgb.gradient_strength_start .* diff' ./ 4 .+ + diff .* diff' ./ 10 + ) .* duration(cgb)^3 +end + + +variables(::Type{<:AbstractChangingGradientBlock}) = [duration, qvec, bmat, slew_rate] + +end diff --git a/src/gradients/constant_gradient_blocks.jl b/src/gradients/constant_gradient_blocks.jl new file mode 100644 index 0000000000000000000000000000000000000000..2f3eea804f3ccdf5a93581a9f904f2a25493b737 --- /dev/null +++ b/src/gradients/constant_gradient_blocks.jl @@ -0,0 +1,58 @@ +module ConstantGradientBlocks +import StaticArrays: SVector +import ...Variables: VariableType, variables +import ...BuildingBlocks: GradientBlock +import JuMP: @constraint, @variable, Model, owner_model +import ...Variables: qvec, bmat, gradient_strength, slew_rate, duration, variables, VariableType +import ...BuildingBlocks: GradientBlock, fixed, RFPulseBlock + +struct AbstractConstantGradientBlock{T<:VariableType} <: GradientBlock + gradient_strength :: SVector{3, T} + duration :: T + rotate :: Union{Nothing, Symbol} + scale :: Union{Nothing, Symbol} +end + +""" + ConstantGradientBlock(duration, gradient_strength, rotate, scale) + +Underlying type for any flat part in a gradient waveform. +It can overlap with an [`RFPulseBlock`], in which case it will be a `SliceSelectPulse` + +The corresponding type in a fixed gradient is [`FixedConstantGradientBlock`](@ref) + +Do not create this directly, it will be created by the gradient waveforms. +""" +const ConstantGradientBlock = AbstractConstantGradientBlock{VariableType} + +""" + FixedConstantGradientBlock(duration, gradient_strength, rotate, scale) + +Fixed equivalent of [`ConstantGradientBlock`](@ref). +""" +const FixedConstantGradientBlock = AbstractConstantGradientBlock{Float64} + +duration(cgb::AbstractConstantGradientBlock) = cgb.duration +gradient_strength(cgb::AbstractConstantGradientBlock) = cgb.gradient_strength +slew_rate(::AbstractConstantGradientBlock) = zero(SVector{3, Float64}) +qvec(cgb::AbstractConstantGradientBlock) = duration(cgb) .* gradient_strength(cgb) + +function bmat(cgb::AbstractConstantGradientBlock) + return (grad .* grad') .* duration(cgb)^3 ./3 +end + +function bmat(cgb::AbstractConstantGradientBlock, qstart) + # \int dt (qstart + t * grad)^2 = + # qstart^2 * duration + + # qstart * grad * duration^2 + + # grad * grad * duration^3 / 3 + + return ( + qstart .* qstart' .* duration(cgb) .+ + qstart .* grad' .* duration(cgb)^2 .+ + bmat(cgb) + ) +end + +variables(::Type{<:AbstractConstantGradientBlock}) = [duration, gradient_strength, qvec, bmat] + +end