Add running cost functional for smoothness

[goerz]

Especially for GRAPE, running costs that penalize the derivative of the control can be useful to force smooth pulses. A prototype implementation would be

function J_a_smoothness(pulsevals, tlist)
    N = length(tlist) - 1  # number of intervals
    L = length(pulsevals) ÷ N
    @assert length(pulsevals) == N * L
    J_a = 0.0
    for l = 1:L
        for n = 2:N
            J_a += (pulsevals[(l-1) * N + n] - pulsevals[(l-1) * N + n - 1])^2
        end
    end
    return 0.5 * J_a
end

function grad_J_a_smoothness(pulsevals, tlist)
    ∇J_a = zeros(length(pulsevals))
    N = length(tlist) - 1  # number of intervals
    L = length(pulsevals) ÷ N
    for l = 1:L
        for n = 1:N
            ∇J_a[(l-1) * N + n] = 0.0
            uₙ = pulsevals[(l-1) * N + n]
            if n > 1
                uₙ₋₁ = pulsevals[(l-1) * N + n - 1]
                ∇J_a[(l-1) * N + n] += (uₙ - uₙ₋₁)
            end
            if n < N
                uₙ₊₁ = pulsevals[(l-1) * N + n + 1]
                ∇J_a[(l-1) * N + n] += (uₙ - uₙ₊₁)
            end
        end
    end
    return ∇J_a
end

However, this may need some tweaking for "normalization", so that the value of the functional is independent of time sampling. See the attached Pluto notebook.

[goerz]

@jabir-ch In case you need this code before I get around making a PR for this