sinc_pulses.jl

module SincPulses
import JuMP: @constraint
import QuadGK: quadgk
import ....BuildSequences: global_model
import ....Variables: duration, amplitude, effective_time, flip_angle, phase, inverse_bandwidth, VariableType, set_simple_constraints!, frequency, make_generic, get_free_variable
import ...AbstractTypes: RFPulseComponent, split_timestep
import ..GenericPulses: GenericPulse

"""
    SincPulse(; Nzeros=3, apodise=true, variables...)

Represents a radio-frequency pulse with a sinc-like amplitude and constant frequency.

## Parameters
- `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
- [`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 <: RFPulseComponent
    apodise :: Bool
    Nzeros :: Tuple{Integer, Integer}
    norm_flip_angle :: Tuple{Float64, Float64}
    amplitude :: VariableType
    phase :: VariableType
    frequency :: VariableType
    lobe_duration :: VariableType
    group :: Union{Nothing, Symbol}
end

function SincPulse(; 
    Nzeros=3, apodise=true,
    amplitude=nothing, phase=nothing, frequency=nothing, lobe_duration=nothing, group=nothing, kwargs...
) 
    if Nzeros isa Number
        Nzeros = (Nzeros, Nzeros)
    end
    res = SincPulse(
        apodise,
        Nzeros,
        integral_nzero.(Nzeros, apodise),
        [get_free_variable(value) for value in (amplitude, phase, frequency, lobe_duration)]...,
        group
    )
    if !(res.amplitude isa Number)
        @constraint global_model() res.amplitude >= 0
    end
    if !(res.lobe_duration isa Number)
        @constraint global_model() res.lobe_duration >= 0
    end
    set_simple_constraints!(res, kwargs)
    return res
end

function normalised_function(x, Nleft, Nright; apodise=false)
    if iszero(x)
        return 1.
    end
    if apodise
        if x < 0
            return (0.54 + 0.46 * cos(π * x / Nleft)) * sin(π * x) / (π * x)
        else
            return (0.54 + 0.46 * cos(π * x / Nright)) * sin(π * x) / (π * x)