function sys_lim = ensure_limits(sys, bandlimit, rolloff, method)
    %ENSURE_LIMITS ensures that the system is strictly proper and
    %bandlimited by moving higher-frequency poles to the bandlimit point or
    %adding new poles. It also ensures that integrators in the system are
    %shifted into the RHP at roll-off frequency.

    arguments
        sys lti
        bandlimit (1,1) double = inf
        rolloff (1,1) double = 0
        method (1,1) string = "move"
    end

    bandlimit = 2*pi * bandlimit;
    rolloff = 2*pi * rolloff;

    sys_lim = zpk(sys);

    p = sys_lim.P{1};
    p = p - rolloff * (real(p) == 0); % Move uncontrollable modes on imaginary axis back to rolloff frequency.
    
    if isfinite(bandlimit)
        idx = (abs(p) > bandlimit);

        if method == "move"
            sys_lim.K = sys_lim.K / prod(abs(p(idx)) / bandlimit);  % Compensate gain
            p(idx) = sign(p(idx)) * bandlimit;  % Ensure bandlimit by moving poles beyond bandlimit to bandlimit frequency.
        elseif method == "cut"
            sys_lim.K = sys_lim.K / prod(abs(p(idx)));
            p = p(~idx);
        else
            error("Invalid method '%s'!", method);
        end

        rel_deg = length(sys_lim.Z{1}) - length(p);
        p = [p; -bandlimit * ones(max(0, rel_deg+1), 1)];  % Ensure no feedthrough w -> z by adding poles at bandlimit frequency till strictly proper.
        sys_lim.K = sys_lim.K * bandlimit^max(0, rel_deg+1);  % Compensate gain
    end
    sys_lim.P{1} = p;
end