% BUILD_OSCILLATOR_MODEL
% Updates system parameters for a generic oscillator.
%
% Model time unit is seconds unless changed by `time_unit`, model output will be scaled in
% 10^sys_params.scaling unless changed by `scaling`. It is assumed that all parameters in sys_params are consistent
% (e.g. sys_params.k_det should be scaled according to sys_params.scaling). Returns new system structure with
% rescaled models.
%
% Arguments
% ---------
% sys_params    : Standardized structure of laser-lock setup parameters.
% Scaling       : How the units of the resulting model should be scaled.
% TimeUnit      : Time unit in which the LTI models are scaled. Default: "seconds"
% FilterMinDecay: Ensures stability of colouring filters by assigning real part to -FilterMinDecay.
% Bandlimit     : Ensures spectral radius of system matrix by moving all poles beyond Bandlimit back to Bandlimit.
%
% Outputs
% -------
% sys: New standardized structure of oscillator parameters.
%
% SPDX-FileCopyrightText: 2023 Deutsches Elektronen-Synchrotron DESY, MSK, LbSync <lbsync-expert@desy.de>
% SPDX-License-Identifier: GPL-3.0-or-later
%
% Author: Maximilian Schütte <maximilian.schuette@desy.de>

function sys = build_oscillator_model(sys_params, args)
    arguments
        sys_params
        args.Scaling        (1,1) double {mustBeFinite, mustBeInteger} = sys_params.scaling
        args.TimeUnit       (1,1) string                               = "seconds"
        args.FilterMinDecay (1,1) double {mustBeNonnegative}           = 0
        args.Bandlimit      (1,1) double {mustBePositive}              = inf
    end

    if mod(args.Scaling, 3)
        error("Scaling must be divisible by 3!");
    end

    time_unit = args.TimeUnit;
    time_unit = ltipack.matchKey(time_unit, ltipack.getValidTimeUnits());
    if isempty(time_unit), error(message('Control:ltiobject:chgunit1')); end
    tconv = tunitconv("seconds", time_unit);

    sys = sys_params;  % Make copy of system

    % Set default parameters, if not provided
    if ~isfield(sys, "pade_order"), sys.pade_order = 1; end
    if ~isfield(sys, "input_scaling"), sys.input_scaling = 1; end

    % Set new scaling
    sys.scaling = args.Scaling;
    sys.scaling_str = num2sip(10^-sys.scaling, [], true);
    sys.scaling_str = extractAfter(sys.scaling_str, find(sys.scaling_str == ' ', 1, 'last')) + "seconds";

    sys.Fd = sys.Fd * 10^(sys.scaling - sys_params.scaling);
    sys.Fd.InputName = "w"; sys.Fd.OutputName = "d";
    sys.Fd.OutputUnit = sys.scaling_str;
    sys.Fd = chgTimeUnit(sys.Fd, time_unit);
    sys.Fd = ensure_limits(sys.Fd, args.Bandlimit / tconv, 1 / tconv, "cut");
    sys.Fd.UserData = sys_params.Fd.UserData;

    sys.Pcol = sys.Fd;
    sys.Pcol.InputName = "w"; sys.Pcol.OutputName = "y";
end
