diff --git a/Project.toml b/Project.toml index 9c5ae506a12baa91242ca3d4347ed57965f92a29..321acf7e5eaa13a7026df6b30951330300985e0c 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.0.0" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" diff --git a/src/gradients/fixed_gradients.jl b/src/gradients/fixed_gradients.jl index 27d904a4d646767239a68323499d3eeb28e3d622..5ee85f8b79f9379d79586d4ffb5cce600e973d50 100644 --- a/src/gradients/fixed_gradients.jl +++ b/src/gradients/fixed_gradients.jl @@ -1,7 +1,9 @@ module FixedGradients -import ...BuildingBlocks: GradientBlock, fixed, BuildingBlock -import ...Variables: variables, duration, qval +import Printf: @sprintf +import LinearAlgebra: norm +import ...BuildingBlocks: GradientBlock, fixed, BuildingBlock, BuildingBlockPrinter +import ...Variables: variables, duration, qval, gradient_strength, slew_rate """ @@ -49,7 +51,36 @@ variables(::Type{<:FixedGradient}) = [] duration(fg::FixedGradient) = maximum(fg.time) -Base.show(io::IO, fb::FixedGradient) = print(io, "FixedGradient for $(duration(fb)) ms") +gradient_strength(fg::FixedGradient) = max(maximum(abs.(fg.Gx)), maximum(abs.(fg.Gy)), maximum(abs.(fg.Gz))) + +slew_rate(fg::FixedGradient) = maximum(map((fg.Gx, fg.Gy, fg.Gz)) do gradient + return maximum(abs.(gradient[2:end] - gradient[1:end-1]) ./ (fg.time[2:end] - fg.time[1:end-1])) +end) + +qvec(fg::FixedGradient) = map((fg.Gx, fg.Gy, fg.Gz)) do gradient + weights_double = fg.time[2:end] - fg.time[1:end-1] + weights = [weights_double[1] / 2, ((weights_double[1:end-1] + weights_double[2:end]) / 2)..., weights_double[end]/2] + return sum(weights .* gradient) +end + +qval(fg::FixedGradient) = norm(qvec(fg)) + + +function Base.show(io::IO, printer::BuildingBlockPrinter{<:FixedGradient}) + fp = printer.bb + pp(value::Number) = @sprintf("%.3g", value) + + print(io, "FixedGradient(") + if isnothing(printer.start_time) + print(io, "duration=", pp(duration(fp))) + else + print(io, "t=", pp(printer.start_time), "-", pp(printer.start_time + duration(fp))) + end + for fn in (qval, gradient_strength, slew_rate) + print(io, ", $(nameof(fn))=", pp(fn(fp))) + end + print(io, ")") +end """ FixedInstantGradient(orientation, qval) diff --git a/src/pulses/fixed_pulses.jl b/src/pulses/fixed_pulses.jl index 49de3fe0bfbba7b86cd5bb05727e95796bc23d39..bbc6332fe34a2d2a5506124404b8ac4eb8d7c4d4 100644 --- a/src/pulses/fixed_pulses.jl +++ b/src/pulses/fixed_pulses.jl @@ -1,7 +1,8 @@ module FixedPulses +import Printf: @sprintf import ...BuildingBlocks: RFPulseBlock, fixed, BuildingBlockPrinter -import ...Variables: variables, duration, amplitude, effective_time +import ...Variables: variables, duration, amplitude, effective_time, flip_angle """ @@ -22,6 +23,7 @@ struct FixedPulse <: RFPulseBlock function FixedPulse(time::AbstractVector{<:Number}, amplitude::AbstractVector{<:Number}, phase::AbstractVector{<:Number}) @assert length(time) == length(amplitude) @assert length(time) == length(phase) + @assert all(time .>= 0) new(Float64.(time), Float64.(amplitude), Float64.(phase)) end end @@ -30,14 +32,29 @@ function FixedPulse(time::AbstractVector{<:Number}, amplitude::AbstractVector{<: return FixedPulse(time, amplitude, (time .- time[1]) .* (frequency * 360) .+ phase) end -variables(::Type{<:FixedPulse}) = [] +variables(::Type{<:FixedPulse}) = [duration, flip_angle, effective_time] duration(fp::FixedPulse) = maximum(fp.time) amplitude(fp::FixedPulse) = maximum(abs.(fp.amplitude)) -effective_time(pulse::FixedPulse) = pulse.time[argmax(abs(pulse.amplitude))] +effective_time(pulse::FixedPulse) = pulse.time[argmax(abs.(pulse.amplitude))] +function flip_angle(pulse::FixedPulse) + weights_double = pulse.time[2:end] - pulse.time[1:end-1] + weights = [weights_double[1] / 2, ((weights_double[1:end-1] + weights_double[2:end]) / 2)..., weights_double[end]/2] + return sum(weights .* pulse.amplitude) * 360 +end -Base.show(io::IO, fp::FixedPulse) = print(io, "FixedPulse for $(duration(fp)) ms (peak at $(effective_time(fp)) ms)") +function Base.show(io::IO, printer::BuildingBlockPrinter{<:FixedPulse}) + fp = printer.bb + pp(value::Number) = @sprintf("%.3g", value) + print(io, "FixedPulse(") + if isnothing(printer.start_time) + print(io, "duration=", pp(duration(fp))) + else + print(io, "t=", pp(printer.start_time), "-", pp(printer.start_time + duration(fp))) + end + print(io, ", flip_angle=", pp(flip_angle(fp)), ", effective_time=", pp(effective_time(fp)), ")") +end """