Convert RF pulses into components

Also, SincPulse now can only have fixed number of zero crossings.

src/components/components.jl

 module Components
 include("abstract_types.jl")
 include("gradient_waveforms/gradient_waveforms.jl")
+include("instant_gradients.jl")
+include("pulses/pulses.jl")
+
+import .AbstractTypes: BaseComponent, GradientWaveform, EventComponent, RFPulseComponent, ReadoutComponent
+import .GradientWaveforms: ConstantGradientBlock, ChangingGradientBlock, NoGradientBlock
+import .InstantGradients: InstantGradient1D, InstantGradient3D
+import .Pulses: GenericPulse, InstantRFPulse, SincPulse, ConstantPulse
 
 end
\ No newline at end of file

src/components/gradient_waveforms/gradient_waveforms.jl

@@ -14,7 +14,7 @@ Each part of this gradient waveform can compute:
 - [`qvec`]: area under curve in each dimension
 - [`bmat_gradient`]: diffusion weighting (scalar in 1D or matrix in 3D).
 """
-module Gradients
+module GradientWaveforms
 
 include("changing_gradient_blocks.jl")
 include("constant_gradient_blocks.jl")

src/components/instant_gradients.jl

 module InstantGradients
 import StaticArrays: SVector, SMatrix
-import ...Variables: VariableType, duration, qvec, bmat_gradient, get_free_variable, set_simple_constraints, qval, effective_time
+import ...Variables: VariableType, duration, qvec, bmat_gradient, get_free_variable, set_simple_constraints!, qval, effective_time
 import ..AbstractTypes: EventComponent, GradientWaveform
 
 """
@@ -34,7 +34,7 @@ end
 
 function InstantGradient1D(; orientation=[1, 0, 0], group=nothing, qval=nothing, variables...)
     res = InstantGradient1D(qval, orientation, group)
-    set_simple_constraints(res, variables)
+    set_simple_constraints!(res, variables)
     return res
 end
 
@@ -67,7 +67,7 @@ function InstantGradient3D(; qvec=[nothing, nothing, nothing], group=nothing, va
         qvec = [nothing, nothing, nothing]
     end
     res = InstantGradient3D(get_free_variable.(qvec), group)
-    set_simple_constraints(res, variables)
+    set_simple_constraints!(res, variables)
     return res
 end
 

src/pulses/constant_pulses.jl → src/components/pulses/constant_pulses.jl

 module ConstantPulses
-import JuMP: VariableRef, @constraint, @variable, value
-import ...BuildingBlocks: RFPulseBlock, set_simple_constraints!
-import ...Variables: variables, get_free_variable, flip_angle, phase, amplitude, frequency, start_time, end_time, VariableType, duration, effective_time, inverse_bandwidth
-import ...BuildSequences: global_model
-import ..GenericPulses: GenericPulse, make_generic
+import JuMP: @constraint
+import ...AbstractTypes: RFPulseComponent
+import ....BuildSequences: global_model
+import ....Variables: duration, amplitude, effective_time, flip_angle, phase, inverse_bandwidth, VariableType, set_simple_constraints!, frequency
+import ..GenericPulses: GenericPulse
 
 """
     ConstantPulse(; variables...)
 
 Represents an radio-frequency pulse with a constant amplitude and frequency (i.e., a rectangular function).
 
-## Properties
+## Parameters
+- `group`: name of the group to which this pulse belongs. This is used for scaling or adding phases/off-resonance frequencies.
+
+## Variables
 - [`flip_angle`](@ref): rotation expected for on-resonance spins in degrees.
 - [`duration`](@ref): duration of the RF pulse in ms.
 - [`amplitude`](@ref): amplitude of the RF pulse in kHz.
 - [`phase`](@ref): phase at the start of the RF pulse in degrees.
 - [`frequency`](@ref): frequency of the RF pulse relative to the Larmor frequency (in kHz).
 """
-struct ConstantPulse <: RFPulseBlock
+struct ConstantPulse <: RFPulseComponent
     amplitude :: VariableType
     duration :: VariableType
     phase :: VariableType
     frequency :: VariableType
-    scale :: Union{Nothing, Symbol}
+    group :: Union{Nothing, Symbol}
 end
 
-function ConstantPulse(amplitude=nothing, duration=nothing, phase=nothing, frequency=nothing, scale=nothing, kwargs...) 
+function ConstantPulse(; amplitude=nothing, duration=nothing, phase=nothing, frequency=nothing, group=nothing, kwargs...) 
     res = ConstantPulse(
         [get_free_variable(value) for value in (amplitude, duration, phase, frequency)]...,
-        scale
-
+        group
     )
     @constraint global_model() res.amplitude >= 0
     set_simple_constraints!(res, kwargs)

src/pulses/generic_pulses.jl → src/components/pulses/generic_pulses.jl

 module GenericPulses
 
-import ...BuildingBlocks: RFPulseBlock, fixed, make_generic
-import ...Variables: variables, duration, amplitude, effective_time, flip_angle
+import ...AbstractTypes: RFPulseComponent
+import ....Variables: duration, amplitude, effective_time, flip_angle
 
 
 """
@@ -9,14 +9,17 @@ import ...Variables: variables, duration, amplitude, effective_time, flip_angle
     GenericPulse(time, amplitude; phase=0., frequency=0., effective_time=<halfway>)
 
 Create a Pulse profile that has been fully defined by N control point.
+
 All arguments should be arrays of the same length N defining these control points.
 
+This pulse has no free variables.
+
 - `time`: time since the start of this [`BuildingBlock`](@ref) in ms.
 - `amplitude`: amplitude of the RF pulse at every timepoint in kHz.
 - `phase`: phase of the RF pulse at every timpoint in degrees. If not set explicitly it will be determined by the provided starting `phase` (degrees) and the `frequency` (kHz).
 - `effective_time`: the time that the RF pulse should be considered to have taken place when computing a `Pathway` (defaults: whenever half of the final flip angle has been achieved for on-resonance spins).
 """
-struct GenericPulse <: RFPulseBlock
+struct GenericPulse <: RFPulseComponent
     time :: Vector{Float64}
     amplitude :: Vector{Float64}
     phase :: Vector{Float64}

src/components/pulses/instant_pulses.jl

+module InstantPulses
+import JuMP: @constraint
+import ...AbstractTypes: RFPulseComponent
+import ....BuildSequences: global_model
+import ....Variables: duration, amplitude, effective_time, flip_angle, phase, inverse_bandwidth, VariableType
+
+"""
+    InstantPulse(; flip_angle=nothing, phase=nothing, group=nothing)
+
+Return an instant RF pulse that rotates all spins by `flip_angle` around an axis that has an angle of `phase` with the X-Y plane.
+
+## Parameters
+- `group`: name of the group to which this pulse belongs. This is used for scaling or adding phases/off-resonance frequencies.
+
+## Variables
+- [`flip_angle`](@ref): angle by which spins are rotated in degrees.
+- [`phase`](@ref): angle of axis around which spins are rotated in degrees.
+"""
+struct InstantPulse <: RFPulseComponent
+    flip_angle :: VariableType
+    phase :: VariableType
+    group :: Union{Nothing, Symbol}
+end
+
+function InstantPulse(; flip_angle=nothing, phase=nothing, group=nothing) 
+    res = InstantPulse(
+        get_free_variable(flip_angle),
+        get_free_variable(phase),
+        group
+    )
+    @constraint global_model() res.flip_angle >= 0
+    return res
+end
+
+flip_angle(instant::InstantPulse) = instant.flip_angle
+phase(instant::InstantPulse) = instant.phase
+duration(::InstantPulse) = 0.
+effective_time(::InstantPulse) = 0.
+inverse_bandwidth(::InstantPulse) = 0.
+
+
+make_generic(block::InstantPulse) = block
+
+end
\ No newline at end of file

src/pulses/pulses.jl → src/components/pulses/pulses.jl

@@ -4,10 +4,9 @@ include("instant_pulses.jl")
 include("constant_pulses.jl")
 include("sinc_pulses.jl")
 
-import ..Variables: effective_time
-import ..BuildingBlocks: RFPulseBlock
+import ..AbstractTypes: RFPulseComponent
 import .GenericPulses: GenericPulse
-import .InstantPulses: InstantRFPulseBlock
+import .InstantPulses: InstantPulse
 import .ConstantPulses: ConstantPulse
 import .SincPulses: SincPulse
 

src/pulses/sinc_pulses.jl → src/components/pulses/sinc_pulses.jl

 module SincPulses
-
-import JuMP: VariableRef, @constraint, @variable
+import JuMP: @constraint
 import QuadGK: quadgk
-import Polynomials: fit, Polynomial
-import ...BuildingBlocks: RFPulseBlock, set_simple_constraints!
-import ...Variables: flip_angle, phase, amplitude, frequency, VariableType, variables, get_free_variable, duration, effective_time, inverse_bandwidth
-import ...BuildSequences: global_model
-import ..GenericPulses: GenericPulse, make_generic
+import ....BuildSequences: global_model
+import ....Variables: duration, amplitude, effective_time, flip_angle, phase, inverse_bandwidth, VariableType, set_simple_constraints!, frequency
+import ...AbstractTypes: RFPulseComponent
+import ..GenericPulses: GenericPulse
 
 """
-    SincPulse(; symmetric=true, max_Nlobes=nothing, apodise=true, variables...)
+    SincPulse(; Nzeros=3, apodise=true, variables...)
 
-Represents an radio-frequency pulse with a constant amplitude and frequency.
+Represents a radio-frequency pulse with a sinc-like amplitude and constant frequency.
 
 ## Parameters
-- `symmetric`: by default the sinc pulse will be symmetric (i.e., `N_left` == `N_right`). Set `symmetric=false` to independently control `N_left` and `N_right`.
-- `max_Nlobes`: by default the computed [`flip_angle`](@ref) is only approximated as `amplitude` * `lobe_duration`. By setting `max_Nlobes` the flip_angle will be given by the actual integral of the sinc function, which will be more accurate for small number of lobes. However, the number of lobes will not be able to exceed this `max_Nlobes`.
+- `Nzeros`: Number of zero-crossings on each side of the sinc pulse. Can be set to a tuple with two values to have a different number of zero crossings on the left and the right of the sinc pulse.
 - `apodise`: if true (default) applies a Hanning apodising window to the sinc pulse.
+- `group`: name of the group to which this pulse belongs. This is used for scaling or adding phases/off-resonance frequencies.
 
 ## Variables
-- [`N_left`](@ref): number of zero-crossings of the sinc function before the main peak (minimum of 1).
-- [`N_right`](@ref): number of zero-crossings of the sinc function after the main peak (minimum of 1).
-- [`flip_angle`](@ref): rotation expected for on-resonance spins in degrees (only approximate if `max_Nlobes` is not set).
+- [`flip_angle`](@ref): rotation expected for on-resonance spins in degrees.
 - [`duration`](@ref): duration of the RF pulse in ms.
 - [`amplitude`](@ref): amplitude of the RF pulse in kHz.
 - [`phase`](@ref): phase at the start of the RF pulse in degrees.
 - [`frequency`](@ref): frequency of the RF pulse relative to the Larmor frequency (in kHz).
 - [`bandwidth`](@ref): width of the rectangular function in frequency space (in kHz). If the `duration` is short (compared with 1/`bandwidth`), this bandwidth will only be approximate.
 """
-struct SincPulse <: RFPulseBlock
+struct SincPulse <: RFPulseComponent
     symmetric :: Bool
     apodise :: Bool
-    nlobe_integral :: Polynomial
-    N_left :: VariableType
-    N_right :: VariableType
+    Nzeros :: Tuple{Integer, Integer}
+    norm_flip_angle :: Tuple{Float64, Float64}
     amplitude :: VariableType
     phase :: VariableType
     frequency :: VariableType
     lobe_duration :: VariableType
-    scale :: Union{Nothing, Symbol}
+    group :: Union{Nothing, Symbol}
 end
 
 function SincPulse(; 
-    symmetric=true, max_Nlobes=nothing, apodise=true, N_lobes=nothing, N_left=nothing, N_right=nothing, 
-    amplitude=nothing, phase=nothing, frequency=nothing, lobe_duration=nothing, scale=nothing, kwargs...
+    Nzeros=3, apodise=true,
+    amplitude=nothing, phase=nothing, frequency=nothing, lobe_duration=nothing, group=nothing, kwargs...
 ) 
-    if symmetric
-        N_lobes = get_free_variable(N_lobes)
-        @assert isnothing(N_left) && isnothing(N_right) "N_left and N_right cannot be set if symmetric=true (default)"
-        N_left_var = N_right_var = N_lobes
-    else
-        @assert isnothing(N_lobes) "N_lobes cannot be set if symmetric=true (default)"
-        N_left_var = get_free_variable(N_left)
-        N_right_var = get_free_variable(N_right)
+    if Nzeros isa Number
+        Nzeros = (Nzeros, Nzeros)
     end
     res = SincPulse(
         symmetric,
         apodise,
-        nlobe_integral_params(max_Nlobes, apodise),
-        N_left_var,
-        N_right_var,
+        Nzeros,
+        integral_nzero.(Nzeros, apodise),
         [get_free_variable(value) for value in (amplitude, phase, frequency, lobe_duration)]...,
-        scale
+        group
     )
-    model = global_model()
-    @constraint model res.amplitude >= 0
-    @constraint model res.N_left >= 1
-    if !symmetric
-        @constraint model res.N_right >= 1
-    end
+    @constraint global_model() res.amplitude >= 0
     set_simple_constraints!(res, kwargs)
     return res
 end
@@ -84,22 +67,18 @@ function normalised_function(x; apodise=false)
     end
 end
 
-function nlobe_integral_params(Nlobe_max, apodise=false)
-    if isnothing(Nlobe_max)
-        return fit([1], [0.5], 0)
-    end
+function integral_nzero(Nzeros, apodise)
     f = x -> normalised_function(x; apodise=apodise)
-    integral_values = [quadgk(f, 0, i)[1] for i in 1:Nlobe_max]
-    return fit(1:Nlobe_max, integral_values, Nlobe_max - 1)
+    return quadgk(f, 0, Nzeros)[1]
 end
 
 amplitude(pulse::SincPulse) = pulse.amplitude
-N_left(pulse::SincPulse) = pulse.N_left
-N_right(pulse::SincPulse) = pulse.N_right
+N_left(pulse::SincPulse) = pulse.Nzeros[1]
+N_right(pulse::SincPulse) = pulse.Nzeros[2]
 duration(pulse::SincPulse) = (N_left(pulse) + N_right(pulse)) * lobe_duration(pulse)
 phase(pulse::SincPulse) = pulse.phase
 frequency(pulse::SincPulse) = pulse.frequency
-flip_angle(pulse::SincPulse) = (pulse.nlobe_integral(N_left(pulse)) + pulse.nlobe_integral(N_right(pulse))) * amplitude(pulse) * lobe_duration(pulse) * 360
+flip_angle(pulse::SincPulse) = (pulse.norm_flip_angle[1] + pulse.norm_flip_angle[2]) * amplitude(pulse) * lobe_duration(pulse) * 360
 lobe_duration(pulse::SincPulse) = pulse.lobe_duration
 inverse_bandwidth(pulse::SincPulse) = lobe_duration(pulse)
 effective_time(pulse::SincPulse) = N_left(pulse) * lobe_duration(pulse)

src/pulses/instant_pulses.jl

-module InstantPulses
-import JuMP: @constraint, @variable, VariableRef, value
-import ...BuildingBlocks: RFPulseBlock, make_generic
-import ...Variables: flip_angle, phase, start_time, variables, duration, get_free_variable, VariableType, effective_time, inverse_bandwidth
-import ...BuildSequences: global_model
-
-struct InstantRFPulseBlock <: RFPulseBlock
-    flip_angle :: VariableType
-    phase :: VariableType
-    scale :: Union{Nothing, Symbol}
-end
-
-function InstantRFPulseBlock(; flip_angle=nothing, phase=nothing, scale=nothing) 
-    res = InstantRFPulseBlock(
-        get_free_variable(flip_angle),
-        get_free_variable(phase),
-        scale
-    )
-    @constraint global_model() res.flip_angle >= 0
-    return res
-end
-
-flip_angle(instant::InstantRFPulseBlock) = instant.flip_angle
-phase(instant::InstantRFPulseBlock) = instant.phase
-duration(::InstantRFPulseBlock) = 0.
-effective_time(::InstantRFPulseBlock) = 0.
-inverse_bandwidth(::InstantRFPulseBlock) = 0.
-
-
-make_generic(block::InstantRFPulseBlock) = block
-
-end
\ No newline at end of file