function sys = build_uncertain_laser_model(sys_params, id, scaling, time_unit)
%build_uncertain_laser_model Constructs a parameterized dynamic model for a laser-lock setup from uncertain system parameters.
% 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_phi should be scaled according to sys_params.scaling). If specified, uncertain system 
% parameters are prefixed with "`id`_" to avoid naming conflicts. Returns new system structure with rescaled models.
    
    arguments
        sys_params
        id (1,:) char = '' %mustBeTextScalar
        scaling (1,1) double {mustBeFinite, mustBeInteger} = sys_params.scaling
        time_unit (1,:) char {mustBeMember(time_unit, {'seconds', 'milliseconds', 'microseconds'})} = 'seconds'
    end
    
    sys = sys_params;  % Make copy of system

    % Set new scaling
    sys.scaling = scaling;
    sys.scaling_str = num2eng(10^-sys.scaling, true, false, false, 5, true);
    sys.scaling_str = extractAfter(sys.scaling_str, find(sys.scaling_str == ' ', 1, 'last')) + "seconds";
    sys.k_phi = sys.k_phi * 10^(sys.scaling - sys_params.scaling);
    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);
    
    if strlength(id) > 0
        sys.id = id;
        sys.k_omega.name = [id '_kw'];
        sys.pz_wn.name = [id '_wn'];
        sys.pz_rwn.name = [id '_rwn'];
        sys.pz_z.name = [id '_z'];
        sys.pz_rz.name = [id '_rz'];
    end
    
    sys.G_pzt = zpk([], -2*pi * sys.pzt_fixed_pole_cfreq, 1);
    if isfinite(sys.pzt_lp_cfreq)
     	sys.G_pzt.P{1} = [sys.G_pzt.P{1}; -2*pi * sys.pzt_lp_cfreq];
    end
    sys.G_pzt.K = sys.G_pzt.K * sys.pzt_gain / dcgain(sys.G_pzt);  % Set DC-Gain to sys.pzt_gain
    sys.G_pzt = chgTimeUnit(sys.G_pzt, time_unit);
    sys.G_pzt.InputName = 'u'; sys.G_pzt.OutputName = 'U';
    sys.G_pzt.InputUnit = 'normalized'; sys.G_pzt.OutputUnit = 'volts';
    
    sys.G_pz = tf([1, 2 * sys.pz_z * sys.pz_rz * sys.pz_wn * sys.pz_rwn, (sys.pz_wn * sys.pz_rwn)^2], [1, 2 * sys.pz_z * sys.pz_wn, sys.pz_wn^2]);
    sys.G_pz = sys.G_pz / dcgain(sys.G_pz); % DC gain of G_pz is absorbed in k_omega
    sys.G_pz = chgTimeUnit(sys.G_pz, time_unit);
    
    % If we use sys.G_pz, we need to divide by sys.G_pz to obtain sys.G_l, which creates an explicit model, which we don't want.
    G_pz = tf([1, 2 * sys.pz_z * sys.pz_wn, sys.pz_wn^2], [1, 2 * sys.pz_z * sys.pz_rz * sys.pz_wn * sys.pz_rwn, (sys.pz_wn * sys.pz_rwn)^2]);
    G_pz = G_pz / dcgain(G_pz); % DC gain of G_pz is absorbed in k_omega
    G_pz = chgTimeUnit(G_pz, time_unit);
    
    % !!! Order of operations has influence on scaling of uncertain params in LFT form
    sys.G_l = 10^sys.scaling / sys.omega_nom * sys.k_omega;
    sys.G_l = sys.G_l * tf(1, [1 0], TimeUnit=time_unit);
    %sys.G_l = sys.G_l * zpk([], -2*pi*1e-1, 1);
    sys.G_l = sys.G_l * G_pz;
    sys.G_l.InputName = 'U'; sys.G_l.OutputName = 'phi';
    
    sys.G = connect(sys.G_pzt, sys.G_l, sys.Fd, sumblk('y = phi + d'), {'w', 'u'}, {'y'}); % Connect to create full model

    sys.C = chgTimeUnit(sys.C, time_unit);
    sys.C.InputName = 'e'; sys.C.OutputName = 'u';
    sys.C.InputUnit = 'normalized'; sys.C.OutputUnit = 'normalized';
    
    sys.CL = connect(sys.G, sys.C, sumblk('e = r - y'), {'r', 'w'}, {'y'});
end
