% Initialization of the Simulink model. The script loads all important
% signals and parameter to the workspace.
%
% Fabian Bredemeier - IAV GmbH
% License: GPL v3
%% Define log index and start + stop time
if ~exist('log_idx', 'var')
    log_idx = 2;
end
tStart = 0;
tEnd = inf;

%% Setting simMode
% 0 = Model validation - Identified models are used and full control loop
% is active
% 1 = Rate Controller validation - No models used, measured input signals
% of rate controller are used
% 2 = Attitude Controller validation - No models used, measured input
% signals of attitude controller are used
% 3 = Altitude Controller validation - No models used, measured input
% signals of attitude controller are used
% 4 = Optimization of controller - Not implemented yet
if ~exist('simMode', 'var')
    simMode = 0;
end

%% Check for existance of data

% Check if sid data is existing
if exist('sidLogs', 'var') % in Workspace
    disp('.mat file already loaded. Reading parameter and signals...');
elseif exist('sid.mat', 'file') % in a .mat file with all other axis tests
    load('sid.mat');
    disp('Loaded data from sid.mat');
elseif ~exist(loadedDataConfig)
    error('Could not find configuration struct. Aborting script.')
else
    error('Could not find data. Aborting script.')
end

%% Config and variable declaration

% Sample rate of simInulation
Fs = 400;
dt = 0.0025;
simIn.param.dt = 0.0025;

% Load data
sid = sidLogs(log_idx).data;
% Check if relevant messages are defined
if (numel(sid) == 0)
    error('Ardupilog object sid is empty. Aborting.');
end

% Get mode of current test run
mode = loadedDataConfig.modes(log_idx);
% Get manual throttle flag for current test run
thr_man = loadedDataConfig.modesConfig(log_idx).thr_man;
% Get filter messages of current test run
msgs = loadedDataConfig.modesConfig(log_idx).filter_msgs;
% Set axis
axis = loadedDataConfig.axisList(log_idx);
% Check if config parameters are complete
if (isempty(mode) || isempty(thr_man) || isempty(msgs))
   error('Config variables (mode, thr_man, msgs) not complete. Aborting.'); 
end

% Define validation and optimization flag of simInulation
simIn.param.mdlVal = simMode == 0;
simIn.param.rateCtrlVal = simMode == 1;
simIn.param.attCtrlVal = simMode == 2 || simMode == 3;
simIn.param.altCtrlVal = simMode == 3 || mode == 2;
simIn.param.optCtrl = simMode == 4;

%% General Copter settings
simIn.param.LOOP_RATE = getParamVal(sid, 'SCHED_LOOP_RATE');
simIn.param.dt = round(double(1 / simIn.param.LOOP_RATE), 4); % Sample time resulting from loop rate. Round to four decimal digits
simIn.param.FRAME_CLASS = getParamVal(sid, 'FRAME_CLASS');
simIn.param.FRAME_TYPE = getParamVal(sid, 'FRAME_TYPE');
[simIn.param.NUM_MOTORS, simIn.param.AXIS_FAC_MOTORS] = getMotorMixerFactors(simIn.param.FRAME_CLASS, simIn.param.FRAME_TYPE);
simIn.param.EKF_TYPE = getParamVal(sid, 'AHRS_EKF_TYPE');
simIn.param.GRAVITY_MSS = single(9.80665);
simIn.param.dt_THR_LOOP = 1 / 50; % Sample time of throttle loop in Copter.cpp. Defined in Copter.cpp, line 94
simIn.param.ANGLE_MAX = getParamVal(sid, 'ANGLE_MAX'); % Maximum lean angle in all flight modes
simIn.param.PILOT_Y_EXPO = getParamVal(sid, 'PILOT_Y_EXPO'); % Acro yaw expo to allow faster rotation when stick at edges
simIn.param.PILOT_Y_RATE = getParamVal(sid, 'PILOT_Y_RATE');
simIn.param.PILOT_Y_RATE_TC = getParamVal(sid, 'PILOT_Y_RATE_TC'); % Pilot yaw rate control input time constant
simIn.param.MODE.ModeNr = mode; % Flight mode

% Abort script if frame configuration is not defined yet
if (isempty(simIn.param.AXIS_FAC_MOTORS) || simIn.param.NUM_MOTORS == 0)
    error("Body configuration of copter is not defined yet!");
end

%% Simulation duration
% Define end time of simInulation
if simIn.param.mdlVal || simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.altCtrlVal
    if isinf(tEnd) % Run full data
        simIn.param.timeEnd = sid.RATE.TimeS(end);
    else % Only run part of data
        duration = tEnd-tStart;
        iEnd = int16(duration/dt);
        simIn.param.timeEnd = sid.RATE.TimeS(iEnd);
        clear duration iEnd
    end
elseif simIn.param.optCtrl
    simIn.param.timeEnd = OPT.tEndOpt;
end

%% Inertial Sensor
simIn.param.INS.GYRO_FILTER = getParamVal(sid, 'INS_GYRO_FILTER'); % Filter cutoff frequency for gyro lowpass filter
simIn.param.INS.GYRO_RATE = 400; % Gyro sampling rate. Actually define by INS_GYRO_RATE, but simInulation is run with 400 Hz

%% AHRS
trim_x = getParamVal(sid, 'AHRS_TRIM_X');
trim_y = getParamVal(sid, 'AHRS_TRIM_Y');
simIn.param.AHRS.FC2BF = single(dcmFromEuler(trim_x, trim_y, 0)); % Rotation matrix from FC to Body Frame based on AP_AHRS constructor
simIn.param.AHRS.BF2FC = simIn.param.AHRS.FC2BF';

clear trim_x trim_y
%% RCIN - Radio Input and parameters for Attitude Control
rollCh = getParamVal(sid, 'RCMAP_ROLL');
pitchCh = getParamVal(sid, 'RCMAP_PITCH');
yawCh = getParamVal(sid, 'RCMAP_YAW');
thrCh = getParamVal(sid, 'RCMAP_THROTTLE');

simIn.param.RCIN.DZ_RLL = getParamVal(sid, ['RC' num2str(rollCh) '_DZ']);
simIn.param.RCIN.MIN_RLL = getParamVal(sid, ['RC' num2str(rollCh) '_MIN']);
simIn.param.RCIN.MAX_RLL = getParamVal(sid, ['RC' num2str(rollCh) '_MAX']);
simIn.param.RCIN.TRIM_RLL = getParamVal(sid, ['RC' num2str(rollCh) '_TRIM']);
simIn.param.RCIN.REVERSED_RLL = getParamVal(sid, ['RC' num2str(rollCh) '_REVERSED']);
simIn.param.RCIN.DZ_PIT = getParamVal(sid, ['RC' num2str(pitchCh) '_DZ']);
simIn.param.RCIN.MIN_PIT = getParamVal(sid, ['RC' num2str(pitchCh) '_MIN']);
simIn.param.RCIN.MAX_PIT = getParamVal(sid, ['RC' num2str(pitchCh) '_MAX']);
simIn.param.RCIN.TRIM_PIT = getParamVal(sid, ['RC' num2str(pitchCh) '_TRIM']);
simIn.param.RCIN.REVERSED_PIT = getParamVal(sid, ['RC' num2str(pitchCh) '_REVERSED']);
simIn.param.RCIN.DZ_YAW = getParamVal(sid, ['RC' num2str(yawCh) '_DZ']);
simIn.param.RCIN.MIN_YAW = getParamVal(sid, ['RC' num2str(yawCh) '_MIN']);
simIn.param.RCIN.MAX_YAW = getParamVal(sid, ['RC' num2str(yawCh) '_MAX']);
simIn.param.RCIN.TRIM_YAW = getParamVal(sid, ['RC' num2str(yawCh) '_TRIM']);
simIn.param.RCIN.REVERSED_YAW = getParamVal(sid, ['RC' num2str(yawCh) '_REVERSED']);
simIn.param.RCIN.DZ_THR = getParamVal(sid, ['RC' num2str(thrCh) '_DZ']);
simIn.param.RCIN.MIN_THR = getParamVal(sid, ['RC' num2str(thrCh) '_MIN']);
simIn.param.RCIN.MAX_THR = getParamVal(sid, ['RC' num2str(thrCh) '_MAX']);
simIn.param.RCIN.TRIM_THR = getParamVal(sid, ['RC' num2str(thrCh) '_TRIM']);
simIn.param.RCIN.REVERSED_THR = getParamVal(sid, ['RC' num2str(thrCh) '_REVERSED']);

simIn.param.RCIN.ROLL_PITCH_YAW_INPUT_MAX = 4500; % Yaw, roll, pitch input range, defined in config.h

% Check if RCIN log is available
if simIn.param.mdlVal || simIn.param.attCtrlVal || simIn.param.altCtrlVal
    if isfield(sid, 'RCIN')
        iVec = (sid.RCIN.TimeS >= tStart & sid.RCIN.TimeS <= tEnd);
        simIn.signals.RCIN.Time = single(sid.RCIN.TimeS(iVec)-tStart);
        simIn.signals.RCIN.RollInput = eval(['single(sid.RCIN.C' num2str(rollCh) '(iVec))']);
        simIn.signals.RCIN.PitchInput = eval(['single(sid.RCIN.C' num2str(pitchCh) '(iVec))']);
        simIn.signals.RCIN.YawInput = eval(['single(sid.RCIN.C' num2str(yawCh) '(iVec))']);
        simIn.signals.RCIN.ThrottleInput = eval(['single(sid.RCIN.C' num2str(thrCh) '(iVec))']);
    else
        if ~isfield(sid, 'RCIN') && axis == 13 
            warning(['RCIN message cannot be found in log ' num2str(log_idx) '! Validating vertical motion model.']);
        else
            warning(['RCIN message cannot be found in log ' num2str(log_idx) '! Ommiting attitude controller and models in simInulation. Only rate controller is validated']);
            simIn.param.rateCtrlVal = 1;
            simIn.param.mdlVal = 0;
        end
        % Set RCIN signals to zero
        simIn.signals.RCIN.Time = single(0:0.1:simIn.param.timeEnd);
        simIn.signals.RCIN.RollInput = single(zeros(length(simIn.signals.RCIN.Time),1));
        simIn.signals.RCIN.PitchInput = single(zeros(length(simIn.signals.RCIN.Time),1));
        simIn.signals.RCIN.YawInput = single(zeros(length(simIn.signals.RCIN.Time),1));
        simIn.signals.RCIN.ThrottleInput = single(zeros(length(simIn.signals.RCIN.Time),1));
    end
end

%% Attitude
% Yaw Angle signals are a modulo of 360°, so this operation has to be
% inverted
iVec = (sid.ATT.TimeS >= tStart & sid.ATT.TimeS <= tEnd);
simIn.signals.ATT.Time = single(sid.ATT.TimeS(iVec)-tStart);
simIn.signals.ATT.Roll = single(sid.ATT.Roll(iVec));
simIn.signals.ATT.Pitch = single(sid.ATT.Pitch(iVec));
simIn.signals.ATT.Yaw = single(demodYawAngle(sid.ATT.Yaw(iVec)));
simIn.signals.ATT.DesRoll = single(sid.ATT.DesRoll(iVec));
simIn.signals.ATT.DesPitch = single(sid.ATT.DesPitch(iVec));
simIn.signals.ATT.DesYaw = single(demodYawAngle(sid.ATT.DesYaw(iVec)));

% The actual yaw angle is calculated from the DCM with
% Matrix3<T>::to_euler() with the four-quadrant arcus tangens atan2(), for 
% example in AP_AHRS_View::update(), so the controller can only handle 
% values in the range [-pi, pi]. 
simIn.signals.ATT.Yaw = simIn.signals.ATT.Yaw - 360 * (simIn.signals.ATT.Yaw > 180);
simIn.signals.ATT.DesYaw = simIn.signals.ATT.DesYaw - 360 * (simIn.signals.ATT.DesYaw > 180);

% Initial DCM matrix based on initial Euler angles
% Calculation based on Matrix3<T>::from_euler()
for i=1:length(iVec)
    if iVec(i) == 1
        iInit = i;
        break;
    end 
end
rollAngInit = single(sid.ATT.Roll(iInit))*pi/180;
pitchAngInit = single(sid.ATT.Pitch(iInit))*pi/180;
yawAngInit = single(simIn.signals.ATT.Yaw(1))*pi/180;
simIn.init.ATT.Roll = rollAngInit;
simIn.init.ATT.Pitch = pitchAngInit;
simIn.init.ATT.Yaw = yawAngInit;
simIn.init.ATT.DCM = dcmFromEuler(rollAngInit, pitchAngInit, yawAngInit);

% Initialize attitudeTarget according to
% AC_AttitudeControl::reset_yaw_target_and_rate(), which is called in
% ModeStabilize::run() when copter is landed
attitudeTargetUpdate = fromAxisAngle([0;0;yawAngInit]);
simIn.init.ATT.attitudeTarget = quatMult(attitudeTargetUpdate, [1;0;0;0]);
simIn.init.ATT.eulerRateTar = [0;0;0];

% Test signals in SITL
% Attitude Controller Targets
% simIn.signals.ATT.attTarMeasW = single(sid.ATAR.attW);
% simIn.signals.ATT.attTarMeasX = single(sid.ATAR.attX);
% simIn.signals.ATT.attTarMeasY = single(sid.ATAR.attY);
% simIn.signals.ATT.attTarMeasZ = single(sid.ATAR.attZ);
% simIn.signals.ATT.eulerAngTarMeasX = single(sid.ATAR.angX);
% simIn.signals.ATT.eulerAngTarMeasY = single(sid.ATAR.angY);
% simIn.signals.ATT.eulerAngTarMeasZ = single(sid.ATAR.angZ);
% simIn.signals.ATT.eulerRatTarMeasX = single(sid.ATAR.ratX);
% simIn.signals.ATT.eulerRatTarMeasY = single(sid.ATAR.ratY);
% simIn.signals.ATT.eulerRatTarMeasZ = single(sid.ATAR.ratZ);
% 
% simIn.signals.ATT.CtrlRollIn = single(sid.CTIN.rllCtrl);
% simIn.signals.ATT.CtrlPitchIn = single(sid.CTIN.pitCtrl);
% simIn.signals.ATT.CtrlYawIn = single(sid.CTIN.yawCtrl);
% simIn.signals.ATT.yawRateDesMeas = single(sid.ATIN.yawDes);

clear rollAngInit pitchAngInit yawAngInit attitudeTargetUpdate iInit

%% Attitude Controller General
simIn.param.ATC.RATE_FF_ENAB = getParamVal(sid, 'ATC_RATE_FF_ENAB');
simIn.param.ATC.INPUT_TC = getParamVal(sid, 'ATC_INPUT_TC'); % Attitude control input time constant
simIn.param.ATC.ANGLE_LIMIT_MIN = single(10); % Min lean angle so that vehicle can maintain limited control, defined in AC_AttitudeControl.cpp as AC_ATTITUDE_CONTROL_ANGLE_LIMIT_MIN
simIn.param.ATC.ANGLE_LIMIT_THROTTLE_MAX = single(0.8); % Max throttle used to limit lean angle so that vehicle does not lose altitude, defined in AC_AttitudeControl.h as AC_ATTITUDE_CONTROL_ANGLE_LIMIT_THROTTLE_MAX
simIn.param.ATC.ANG_LIM_TC = getParamVal(sid, 'ATC_ANG_LIM_TC'); % Angle Limit (to maintain altitude) Time Constant
% Roll Angle Controller
simIn.param.ATC.ANG_RLL_P = getParamVal(sid, 'ATC_ANG_RLL_P');
simIn.param.ATC.ACCEL_R_MAX = getParamVal(sid, 'ATC_ACCEL_R_MAX');
simIn.param.ATC.RATE_R_MAX = getParamVal(sid, 'ATC_RATE_R_MAX');
simIn.param.ATC.ACCEL_RP_CONTROLLER_MIN_RADSS = 40*pi/180; % Maximum body-frame acceleration limit for the stability controller, defined in AC_AttitudeControl.h
simIn.param.ATC.ACCEL_RP_CONTROLLER_MAX_RADSS = 720*pi/180;
% Pitch Angle Controller
simIn.param.ATC.ANG_PIT_P = getParamVal(sid, 'ATC_ANG_PIT_P');
simIn.param.ATC.ACCEL_P_MAX = getParamVal(sid, 'ATC_ACCEL_P_MAX');
simIn.param.ATC.RATE_P_MAX = getParamVal(sid, 'ATC_RATE_P_MAX');
% Yaw Angle Controller
simIn.param.ATC.ANG_YAW_P = getParamVal(sid, 'ATC_ANG_YAW_P');
simIn.param.ATC.ACCEL_Y_MAX = getParamVal(sid, 'ATC_ACCEL_Y_MAX');
simIn.param.ATC.RATE_Y_MAX = getParamVal(sid, 'ATC_RATE_Y_MAX');
simIn.param.ATC.ACCEL_Y_CONTROLLER_MAX_RADSS = single(120*pi/180); % Maximum body-frame acceleration limit for the stability controller, defined in AC_AttitudeControl.h
simIn.param.ATC.ACCEL_Y_CONTROLLER_MIN_RADSS = single(10*pi/180);
simIn.param.ATC.THRUST_ERROR_ANGLE = single(30*pi/180); % Thrust angle error above which yaw corrections are limited. Defined in AC_AttitudeControl.h

% Thrust and throttle calculation parameters
simIn.param.MODE.THR_MAN_FLG = thr_man;    % Flag for manual throttle, which depends on the current mode. 
if thr_man
    simIn.param.ATC.THR_FILT_CUTOFF = getParamVal(sid, 'PILOT_THR_FILT'); % Cutoff frequency of throttle filter for modes with manual throttle. parameter is given to function call of AC_AttitudeControl_Multi::set_throttle_out().
else
    simIn.param.ATC.THR_FILT_CUTOFF = 2.0; % Cutoff frequency of throttle filter for modes with z position control (Defined in AC_PosControl.h, line 33)
end
simIn.param.ATC.ANGLE_BOOST = getParamVal(sid, 'ATC_ANGLE_BOOST');   % Enabling of angle boost to increase output throttle. Used in AC_AttitudeConrtol_Multi::get_throttle_boosted()
simIn.param.ATC.THR_MIX_MAN = getParamVal(sid, 'ATC_THR_MIX_MAN');   % Throttle vs. attitude priorisation during manual flight
simIn.param.ATC.THR_MIX_MIN = getParamVal(sid, 'ATC_THR_MIX_MIN');   % Throttle vs. attitude priorisation used when landing
simIn.param.ATC.THR_MIX_MAX = getParamVal(sid, 'ATC_THR_MIX_MAX');   % Throttle vs. attitude priorisation during active flights
simIn.param.ATC.THR_MIX_DFLT = single(0.5); % Default value for the priorization of throttle use between attitude control and throttle demand by pilot (or autopilot). Defined in AC_AttitudeControl.h, line 44 
if thr_man
    simIn.init.ATC.THR_MIX = simIn.param.ATC.THR_MIX_MAN;
else
    simIn.init.ATC.THR_MIX = simIn.param.ATC.THR_MIX_DFLT;
end
simIn.param.MOT.HOVER_LEARN = getParamVal(sid, 'MOT_HOVER_LEARN');     % Enable/Disable automatic learning of hover throttle (0=Disabled, 1=Learn, 2=Learn and Save
simIn.param.MOT.THST_HOVER = getParamVal(sid, 'MOT_THST_HOVER');    % Motor thrust needed to hover. Default value of _throttle_hover.
simIn.param.MOT.THST_HOVER = simIn.param.MOT.THST_HOVER(1); % Assign first value of array to default value which is the correct value in the case, that the parameter is defined twice in the param file.
simIn.param.MOT.THST_HOVER_TC = single(10.0);                                         % Time constant used to update estimated hover throttle, found as AP_MOTORS_THST_HOVER_TC in AP_MotorsMulticopter.h
simIn.param.MOT.THST_HOVER_MIN = single(0.125);                                       % Minimum possible hover throttle, found as AP_MOTORS_THST_HOVER_MIN in AP_MotorsMulticopter.h
simIn.param.MOT.THST_HOVER_MAX = single(0.6875);                                      % Maximum possible hover throttle, found as AP_MOTORS_THST_HOVER_MAX in AP_MotorsMulticopter.h
simIn.param.MOT.BAT_VOLT_MAX = getParamVal(sid, 'MOT_BAT_VOLT_MAX');   % Maximum voltage above which no additional scaling on thrust is performed
simIn.param.MOT.BAT_VOLT_MIN = getParamVal(sid, 'MOT_BAT_VOLT_MIN');   % Minimum voltage below which no additional scaling on thrust is performed
simIn.param.MOT.BAT_CURR_MAX = getParamVal(sid, 'MOT_BAT_CURR_MAX');   % Maximum current over which maximum throttle is limited (and no further scaling is performed)
simIn.param.MOT.THST_EXPO = getParamVal(sid, 'MOT_THST_EXPO');             % Motor thrust curve exponent (0.0 for linear to 1.0 for second order curve)
simIn.param.MOT.BAT_CURR_TC = getParamVal(sid, 'MOT_BAT_CURR_TC');               % Time constant used to limit the maximum current
simIn.param.MOT.YAW_HEADROOM = getParamVal(sid, 'MOT_YAW_HEADROOM');             % Yaw control is goven at least this PWM in microseconds range

% Throttle inputs
iVec = (sid.CTUN.TimeS >= tStart & sid.CTUN.TimeS <= tEnd);
simIn.signals.CTUN.Time = single(sid.CTUN.TimeS(iVec)--tStart);
simIn.signals.CTUN.ThrIn = single(sid.CTUN.ThI(iVec)); % Throttle In, pilots input to Attitude Control (attitude_control->get_throttle_in())
simIn.signals.CTUN.ThrOut = single(sid.CTUN.ThO(iVec)); % Throttle Out, throttle given to motor mixer after filtering (motors->get_throttle())
simIn.signals.CTUN.AngBst = single(sid.CTUN.ABst(iVec)); % Extra amount of throttle, added to throttle_in in AC_AttitudeControl_Multi::get_throttle_boosted()
simIn.signals.CTUN.ThrHov = single(sid.CTUN.ThH(iVec)); % Estimated throttle required to hover throttle in range 0-1 (AP_MotorsMulticopter::get_throttle_hover())

% Throttle inital values
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal || simIn.param.altCtrlVal
    simIn.init.CTUN.ThrOut = single(sid.CTUN.ThO(1));
    simIn.init.CTUN.ThrHov = single(sid.CTUN.ThH(1));
else
    simIn.init.CTUN.ThrOut = single(0);
    simIn.init.CTUN.ThrHov = single(simIn.param.MOT.THST_HOVER);
end

% Altitude signals
simIn.signals.CTUN.DSAlt = single(sid.CTUN.DSAlt(iVec)); % Desired rangefinder altitude
simIn.signals.CTUN.DAlt = single(sid.CTUN.DAlt(iVec)); % Desired altitude
simIn.signals.CTUN.Alt = single(sid.CTUN.Alt(iVec)); % Achieved altitude (inav)
simIn.signals.CTUN.SAlt = single(sid.CTUN.SAlt(iVec)); % Achieved rangefinder altitude
simIn.signals.CTUN.TAlt = single(sid.CTUN.TAlt(iVec)); % Terrain altitude
simIn.signals.CTUN.BAlt = single(sid.CTUN.BAlt(iVec)); % Barometric altitude
simIn.signals.CTUN.CRt = single(sid.CTUN.CRt(iVec)); % Climb rate inav
simIn.signals.CTUN.DCRt = single(sid.CTUN.DCRt(iVec)); % Desired climb rate 

% Altitude inital values
simIn.init.CTUN.Alt = single(sid.CTUN.Alt(1));
simIn.init.CTUN.CRt = single(sid.CTUN.CRt(1)*0.01);

% Intermediate throttle calculation quantities - only read if message
% exists
if isfield(sid, 'MOTQ')
    iVec = (sid.MOTQ.TimeS >= tStart & sid.MOTQ.TimeS <= tEnd);
    simIn.signals.MOTQ.Time = single(sid.MOTQ.TimeS(iVec)-tStart);
    simIn.signals.MOTQ.ThrAvgMax = single(sid.MOTQ.ThAvgMax(iVec));
    simIn.signals.MOTQ.ThrThstMaxMeas = single(sid.MOTQ.ThThstMax(iVec));
    simIn.signals.MOTQ.CompGain = single(sid.MOTQ.CompGain(iVec));
    simIn.signals.MOTQ.AirDensityRatio = single(sid.MOTQ.AirDensRat(iVec));
    simIn.signals.MOTQ.ThrMixOut = single(sid.MOTQ.ThrOut(iVec));
else
    iVec = (sid.CTUN.TimeS >= tStart & sid.CTUN.TimeS <= tEnd);
    simIn.signals.MOTQ.Time = single(sid.CTUN.TimeS(iVec)-tStart);
    simIn.signals.MOTQ.ThrAvgMax = single(zeros(length(simIn.signals.MOTQ.Time),1));
    simIn.signals.MOTQ.ThrThstMaxMeas = single(zeros(length(simIn.signals.MOTQ.Time),1));
    simIn.signals.MOTQ.CompGain = single(zeros(length(simIn.signals.MOTQ.Time),1));
    simIn.signals.MOTQ.AirDensityRatio = single(zeros(length(simIn.signals.MOTQ.Time),1));
    simIn.signals.MOTQ.ThrMixOut = single(zeros(length(simIn.signals.MOTQ.Time),1));
end


% Battery inputs
iVec = (sid.BAT.TimeS >= tStart & sid.BAT.TimeS <= tEnd);
simIn.signals.BAT.Time = single(sid.BAT.TimeS(iVec)-tStart);
simIn.signals.BAT.RestVoltEst = single(sid.BAT.VoltR(iVec)); % Estimated resting voltage of the battery
simIn.signals.BAT.Curr = single(sid.BAT.Curr(iVec)); % Measured battery current
simIn.signals.BAT.Res = single(sid.BAT.Res(iVec)); % Estimated battery resistance
simIn.signals.BAT.Volt = single(sid.BAT.Volt(iVec)); % Measured voltage resistance

% Battery inital values
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal || simIn.param.altCtrlVal
    simIn.init.BAT.RestVoltEst = single(sid.BAT.VoltR(1));
else
    simIn.init.BAT.RestVoltEst = single(0);
end

% Baro inputs
iVec = (sid.BARO.TimeS >= tStart & sid.BARO.TimeS <= tEnd);
simIn.signals.BARO.Time = single(sid.BARO.TimeS(iVec)-tStart);
simIn.signals.BARO.Alt = single(sid.BARO.Alt(iVec));     % Calculated altitude by the barometer
simIn.signals.BARO.AirPress = single(sid.BARO.Press(iVec));  % Measured atmospheric pressure by the barometer
simIn.signals.BARO.GndTemp = single(sid.BARO.GndTemp(iVec)); % Temperature on ground, specified by parameter or measured while on ground

% Attitude controller outputs
iVec = (sid.MOTB.TimeS >= tStart & sid.MOTB.TimeS <= tEnd);
simIn.signals.MOTB.Time = single(sid.MOTB.TimeS(iVec)-tStart);
simIn.signals.MOTB.LiftMax = single(sid.MOTB.LiftMax(iVec)); % Maximum motor compensation gain, calculated in AP_MotorsMulticopter::update_lift_max_from_batt_voltage()
simIn.signals.MOTB.ThrLimit = single(sid.MOTB.ThLimit(iVec)); % Throttle limit set due to battery current limitations, calculated in AP_MotorsMulticopter::get_current_limit_max_throttle()
simIn.signals.MOTB.ThrAvMx = single(sid.MOTB.ThrAvMx(iVec)); % Throttle average max

% SID Inputs
% Create array for sid signals, each element containing the signal struct
if (mode == 25)
    iVec = (sid.SIDD.TimeS >= tStart & sid.SIDD.TimeS <= tEnd);
    for i=1:13
        simIn.signals.SIDD.Time{i} = single(sid.SIDD.TimeS(iVec)-tStart);
        if (i == axis)
            simIn.signals.SIDD.Chirp{i} = single(sid.SIDD.Targ(iVec)); % SID test signal
        else
            simIn.signals.SIDD.Chirp{i} = single(zeros(length(simIn.signals.SIDD.Time{i}),1));
        end
    end
    simIn.signals.SIDD.AccX = single(sid.SIDD.Ax(iVec));
    simIn.signals.SIDD.AccY = single(sid.SIDD.Ay(iVec));
    simIn.signals.SIDD.AccZ = single(sid.SIDD.Az(iVec));
    simIn.signals.SIDD.GyrX = single(sid.SIDD.Gx(iVec));
    simIn.signals.SIDD.GyrY = single(sid.SIDD.Gy(iVec));
    simIn.signals.SIDD.GyrZ = single(sid.SIDD.Gz(iVec));
else
    iVec = (sid.RATE.TimeS >= tStart & sid.RATE.TimeS <= tEnd);
    for i=1:13
        simIn.signals.SIDD.Time{i} = single(sid.RATE.TimeS(iVec)-tStart);
        simIn.signals.SIDD.Chirp{i} = single(zeros(length(simIn.signals.SIDD.Time{i}),1));
    end
    simIn.signals.SIDD.AccX = single(zeros(length(simIn.signals.SIDD.Time{1}),1));
    simIn.signals.SIDD.AccY = single(zeros(length(simIn.signals.SIDD.Time{1}),1));
    simIn.signals.SIDD.AccZ = single(zeros(length(simIn.signals.SIDD.Time{1}),1));
    simIn.signals.SIDD.GyrX = single(zeros(length(simIn.signals.SIDD.Time{1}),1));
    simIn.signals.SIDD.GyrY = single(zeros(length(simIn.signals.SIDD.Time{1}),1));
    simIn.signals.SIDD.GyrZ = single(zeros(length(simIn.signals.SIDD.Time{1}),1));
end

%% Roll Rate Controller
iVec = (sid.PIDR.TimeS >= tStart & sid.PIDR.TimeS <= tEnd);
simIn.signals.PIDR.Time = single(sid.PIDR.TimeS(iVec)-tStart);
% Inputs
simIn.signals.PIDR.Tar = single(sid.PIDR.Tar(iVec)); % target values filtered
simIn.signals.PIDR.Act = single(sid.PIDR.Act(iVec)); % actual values
simIn.signals.PIDR.Err = single(sid.PIDR.Err(iVec)); % error values
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    % Clock of Slew Limiter in ms
    % Use tracked logging time of PID message in validation modes
    simIn.signals.PIDR.ClockDMod = single(1000*simIn.signals.PIDR.Time); 
else
    simIn.signals.PIDR.ClockDMod = single(1000 * simIn.param.dt * (0:length(simIn.signals.PIDR.Time)-1)');
end
    
% Outputs
simIn.signals.PIDR.FF = single(sid.PIDR.FF(iVec));
simIn.signals.PIDR.P = single(sid.PIDR.P(iVec));
simIn.signals.PIDR.I = single(sid.PIDR.I(iVec));
simIn.signals.PIDR.D = single(sid.PIDR.D(iVec));
simIn.signals.PIDR.DMod = single(sid.PIDR.Dmod(iVec));
simIn.signals.PIDR.ILimit = single(sid.PIDR.Limit(iVec));
simIn.signals.PIDR.SRate = single(sid.PIDR.SRate(iVec)); % Output slew rate of the slew limiter, stored in _output_slew_rate in AC_PID.cpp, line 161

% parameters - Read from PARM.Value cell array with logical indexing
simIn.param.PIDR.TC = getParamVal(sid, 'ATC_INPUT_TC');
simIn.param.PIDR.FLTT_f = getParamVal(sid, 'ATC_RAT_RLL_FLTT');
simIn.param.PIDR.FLTE_f = getParamVal(sid, 'ATC_RAT_RLL_FLTE');
simIn.param.PIDR.FLTD_f = getParamVal(sid, 'ATC_RAT_RLL_FLTD');
simIn.param.PIDR.P = getParamVal(sid, 'ATC_RAT_RLL_P');
simIn.param.PIDR.I = getParamVal(sid, 'ATC_RAT_RLL_I');
simIn.param.PIDR.D = getParamVal(sid, 'ATC_RAT_RLL_D');
simIn.param.PIDR.IMAX = getParamVal(sid, 'ATC_RAT_RLL_IMAX');
simIn.param.PIDR.FF = getParamVal(sid, 'ATC_RAT_RLL_FF'); %0.05;
simIn.param.PIDR.SR_MAX = getParamVal(sid, 'ATC_RAT_RLL_SMAX'); %5.0;
simIn.param.PIDR.SR_TAU = single(1.0); % Slew Rate Tau - not yet available as a parameter of the copter. Set to 1.0 by default (AC_PID.h, line 24).
simIn.param.PIDR.SR_FLT_f = single(25.0);  % Slew Rate lowpass filter cutoff frequency. Defined in SlewLimiter.h, line 14.

% Inital inputs
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    simIn.init.PIDR.P = single(sid.PIDR.P(1));
    simIn.init.PIDR.I = single(sid.PIDR.I(1));
    simIn.init.PIDR.D = single(sid.PIDR.D(1));
    simIn.init.PIDR.Tar = single(sid.RATE.RDes(1)*pi/180); % Convert to radian due to conversion to degree for logging (AP_AHRS_View::Write_Rate)
    simIn.init.PIDR.TarFilt = single(sid.PIDR.Tar(1));
    simIn.init.PIDR.ErrFilt = single(sid.PIDR.Err(1));
    simIn.init.PIDR.SROut = single(sid.PIDR.SRate(1)); % Initial value of the slew rate determined by the slew limiter. Used for both modifier_slew_rate and output_slew_rate.
    if (sid.PIDR.D(1) ~= 0) % Prevent division by zero
        simIn.init.PIDR.DerFilt = single(sid.PIDR.D(1))/simIn.param.PIDR.D;
    else
        simIn.init.PIDR.DerFilt = single(0);
    end
else % Set to zero if we do not want to validate model with test run
    simIn.init.PIDR.P = single(0);
    simIn.init.PIDR.I = single(0);
    simIn.init.PIDR.D = single(0);
    simIn.init.PIDR.Tar = single(0);
    simIn.init.PIDR.TarFilt = single(0);
    simIn.init.PIDR.ErrFilt = single(0);
    simIn.init.PIDR.DerFilt = single(0);
    simIn.init.PIDR.SROut = single(0);
end

%% Pitch Rate Controller
iVec = (sid.PIDP.TimeS >= tStart & sid.PIDP.TimeS <= tEnd);
simIn.signals.PIDP.Time = single(sid.PIDP.TimeS(iVec)-tStart);
% Inputs
simIn.signals.PIDP.Tar = single(sid.PIDP.Tar(iVec)); % target values filtered
simIn.signals.PIDP.Act = single(sid.PIDP.Act(iVec)); % actual values
simIn.signals.PIDP.Err = single(sid.PIDP.Err(iVec)); % error values                                          
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    % Clock of Slew Limiter
    % Use tracked logging time of PID message in validation modes
    simIn.signals.PIDP.ClockDMod = single(1000 * simIn.signals.PIDP.Time); 
else
    simIn.signals.PIDP.ClockDMod = single(1000 * simIn.param.dt * (0:length(simIn.signals.PIDP.Time)-1)');
end

% Outputs
simIn.signals.PIDP.FF = single(sid.PIDP.FF(iVec));
simIn.signals.PIDP.P = single(sid.PIDP.P(iVec));
simIn.signals.PIDP.I = single(sid.PIDP.I(iVec));
simIn.signals.PIDP.D = single(sid.PIDP.D(iVec));
simIn.signals.PIDP.DMod = single(sid.PIDP.Dmod(iVec));
simIn.signals.PIDP.ILimit = single(sid.PIDP.Limit(iVec));
simIn.signals.PIDP.SRate = single(sid.PIDP.SRate(iVec)); % Output slew rate of the slew limiter, stored in _output_slew_rate in AC_PID.cpp, line 161

% Parameters - Read from PARM.Value cell array with logical indexing
simIn.param.PIDP.TC = getParamVal(sid, 'ATC_INPUT_TC');
simIn.param.PIDP.FLTT_f = getParamVal(sid, 'ATC_RAT_PIT_FLTT');
simIn.param.PIDP.FLTE_f = getParamVal(sid, 'ATC_RAT_PIT_FLTE');
simIn.param.PIDP.FLTD_f = getParamVal(sid, 'ATC_RAT_PIT_FLTD');
simIn.param.PIDP.P = getParamVal(sid, 'ATC_RAT_PIT_P');
simIn.param.PIDP.I = getParamVal(sid, 'ATC_RAT_PIT_I');
simIn.param.PIDP.D = getParamVal(sid, 'ATC_RAT_PIT_D');
simIn.param.PIDP.IMAX = getParamVal(sid, 'ATC_RAT_PIT_IMAX');
simIn.param.PIDP.FF = getParamVal(sid, 'ATC_RAT_PIT_FF');
simIn.param.PIDP.SR_MAX = getParamVal(sid, 'ATC_RAT_PIT_SMAX');
simIn.param.PIDP.SR_TAU = single(1.0); % Slew Rate Tau - not yet available as a parameter of the copter. Set to 1.0 by default (AC_PID.h, line 24).
simIn.param.PIDP.SR_FLT_f = single(25.0);  % Slew Rate lowpass filter cutoff frequency. Defined in SlewLimiter.h, line 14.

% Inital inputs
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    simIn.init.PIDP.P = single(sid.PIDP.P(1));
    simIn.init.PIDP.I = single(sid.PIDP.I(1));
    simIn.init.PIDP.D = single(sid.PIDP.D(1));
    simIn.init.PIDP.Tar = single(sid.RATE.PDes(1)*pi/180); % Convert to radian due to conversion to degree for logging (AP_AHRS_View::Write_Rate)
    simIn.init.PIDP.TarFilt = single(sid.PIDP.Tar(1));
    simIn.init.PIDP.ErrFilt = single(sid.PIDP.Err(1));
    simIn.init.PIDP.SrOut = single(sid.PIDP.SRate(1)); % Initial value of the slew rate determined by the slew limiter. Used for both modifier_slew_rate and output_slew_rate.
    if (simIn.param.PIDP.D ~= 0)
        simIn.init.PIDP.DerFilt = single(sid.PIDP.D(1))/simIn.param.PIDP.D;
    else
        simIn.init.PIDP.DerFilt = single(0);
    end
else
    simIn.init.PIDP.P = single(0);
    simIn.init.PIDP.I = single(0);
    simIn.init.PIDP.D = single(0);
    simIn.init.PIDP.Tar = single(0);
    simIn.init.PIDP.TarFilt = single(0);
    simIn.init.PIDP.ErrFilt = single(0);
    simIn.init.PIDP.DerFilt = single(0);
    simIn.init.PIDP.SrOut = single(0);
end

%% Yaw Rate Controller
iVec = (sid.PIDY.TimeS >= tStart & sid.PIDY.TimeS <= tEnd);
simIn.signals.PIDY.Time = single(sid.PIDY.TimeS(iVec)-tStart);
% Inputs
simIn.signals.PIDY.Tar = single(sid.PIDY.Tar(iVec)); % target values filtered
simIn.signals.PIDY.Act = single(sid.PIDY.Act(iVec)); % actual values
simIn.signals.PIDY.Err = single(sid.PIDY.Err(iVec)); % error values                                          
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    % Clock of Slew Limiter
    % Use tracked logging time of PID message in validation modes
    simIn.signals.PIDY.ClockDMod = single(1000 * simIn.signals.PIDY.Time); 
else
    simIn.signals.PIDY.ClockDMod = single(1000 * simIn.param.dt * (0:length(simIn.signals.PIDY.Time)-1)');
end

% Outputs
simIn.signals.PIDY.FF = single(sid.PIDY.FF(iVec));
simIn.signals.PIDY.P = single(sid.PIDY.P(iVec));
simIn.signals.PIDY.I = single(sid.PIDY.I(iVec));
simIn.signals.PIDY.D = single(sid.PIDY.D(iVec));
simIn.signals.PIDY.DMod = single(sid.PIDY.Dmod(iVec));
simIn.signals.PIDY.ILimit = single(sid.PIDY.Limit(iVec));
simIn.signals.PIDY.SRate = single(sid.PIDY.SRate(iVec)); % Output slew rate of the slew limiter, stored in _output_slew_rate in AC_PID.cpp, line 161

% Parameters - Read from PARM.Value cell array with logical indexing
simIn.param.PIDY.TC = getParamVal(sid, 'ATC_INPUT_TC');
simIn.param.PIDY.FLTT_f = getParamVal(sid, 'ATC_RAT_YAW_FLTT');
simIn.param.PIDY.FLTE_f = getParamVal(sid, 'ATC_RAT_YAW_FLTE');
simIn.param.PIDY.FLTD_f = getParamVal(sid, 'ATC_RAT_YAW_FLTD');
simIn.param.PIDY.P = getParamVal(sid, 'ATC_RAT_YAW_P');
simIn.param.PIDY.I = getParamVal(sid, 'ATC_RAT_YAW_I');
simIn.param.PIDY.D = getParamVal(sid, 'ATC_RAT_YAW_D');
simIn.param.PIDY.IMAX = getParamVal(sid, 'ATC_RAT_YAW_IMAX');
simIn.param.PIDY.FF = getParamVal(sid, 'ATC_RAT_YAW_FF');
simIn.param.PIDY.SR_MAX = getParamVal(sid, 'ATC_RAT_YAW_SMAX');
simIn.param.PIDY.SR_TAU = single(1.0); % Slew Rate Tau - not yet available as a parameter of the copter. Set to 1.0 by default (AC_PID.h, line 24).
simIn.param.PIDY.SR_FLT_f = single(25.0);  % Slew Rate lowpass filter cutoff frequency. Defined in SlewLimiter.h, line 14.

% Inital inputs
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    simIn.init.PIDY.P = single(sid.PIDY.P(1));
    simIn.init.PIDY.I = single(sid.PIDY.I(1));
    simIn.init.PIDY.D = single(sid.PIDY.D(1));
    simIn.init.PIDY.Tar = single(sid.RATE.PDes(1)*pi/180); % Convert to radian due to conversion to degree for logging (AP_AHRS_View::Write_Rate)
    simIn.init.PIDY.TarFilt = single(sid.PIDY.Tar(1));
    simIn.init.PIDY.ErrFilt = single(sid.PIDY.Err(1));
    simIn.init.PIDY.SrOut = single(sid.PIDY.SRate(1)); % Initial value of the slew rate determined by the slew limiter. Used for both modifier_slew_rate and output_slew_rate.
    if (simIn.param.PIDY.D ~= 0)
        simIn.init.PIDY.DerFilt = single(sid.PIDY.D(1))/simIn.param.PIDY.D;
    else
        simIn.init.PIDY.DerFilt = single(0);
    end
else
    simIn.init.PIDY.P = single(0);
    simIn.init.PIDY.I = single(0);
    simIn.init.PIDY.D = single(0);
    simIn.init.PIDY.Tar = single(0);
    simIn.init.PIDY.TarFilt = single(0);
    simIn.init.PIDY.ErrFilt = single(0);
    simIn.init.PIDY.DerFilt = single(0);
    simIn.init.PIDY.SrOut = single(0);
end

%% RATE message
iVec = (sid.RATE.TimeS >= tStart & sid.RATE.TimeS <= tEnd);
simIn.signals.RATE.Time = single(sid.RATE.TimeS(iVec)-tStart);

% Roll
simIn.signals.RATE.ROut = single(sid.RATE.ROut(iVec)); % Rate Controller total output (except FF)
% Rate Controller target (deg/s): Logged signal in val phase, step in opt phase.
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    simIn.signals.RATE.RDes = single(sid.RATE.RDes(iVec));
elseif simIn.param.optCtrl && axis == 10
    simIn.signals.RATE.RDes = single((simIn.signals.RATE.Time > OPT.stepTime)'*OPT.stepMag); % Step at 5s
else
    simIn.signals.RATE.RDes = single(zeros(numel(simIn.signals.RATE.Time), 1)); % Set to zero if axis not optimized
end

% Pitch
simIn.signals.RATE.POut = single(sid.RATE.POut(iVec));
% Rate Controller target (deg/s): Logged signal in val phase, step in opt phase
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    simIn.signals.RATE.PDes = single(sid.RATE.PDes(iVec));
elseif simIn.param.optCtrl && axis == 11
    simIn.signals.RATE.PDes = single((simIn.signals.RATE.Time > OPT.stepTime)'*OPT.stepMag); % Step at 5s
else
    simIn.signals.RATE.PDes = single(zeros(numel(simIn.signals.RATE.Time), 1)); % Set to zero if axis not optimized
end

% Yaw 
simIn.signals.RATE.YOut = single(sid.RATE.YOut(iVec));
% Rate Controller target (deg/s): Logged signal in val phase, step in opt phase
if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal
    simIn.signals.RATE.YDes = single(sid.RATE.YDes(iVec));
elseif simIn.param.optCtrl && axis == 12
    simIn.signals.RATE.YDes = single((simIn.signals.RATE.Time > OPT.stepTime)'*OPT.stepMag); % Step at 5s
else
    simIn.signals.RATE.YDes = single(zeros(numel(simIn.signals.RATE.Time), 1)); % Set to zero if axis not optimized
end

% Throttle
simIn.signals.RATE.AOut = single(sid.RATE.AOut(iVec)); % Throttle command to motors
simIn.signals.RATE.A = single(sid.RATE.A(iVec)); % Vertical acceleration in earth frame

%% Position Controller
% Default definitions
simIn.param.PSCD.THR_DZ = getParamVal(sid, 'THR_DZ'); % Deadzone above and below mid throttle in PWM microseconds. Used in Althold, Loiter, PosHold. Defined in Parameters.cpp
simIn.param.PSCD.VEL_MAX_DOWN_DFLT = single(-150); % Default descent rate in cm/s, defined in AC_PosControl.h, line 27
simIn.param.PSCD.VEL_MAX_UP_DFLT = single(250); % Default climb rate in cm/s, defined in AC_PosControl.h, line 28
simIn.param.PSCD.ACC_MAX_Z_DFLT = single(250); % Default vertical acceleration cm/s/s, defined in AC_PosControl.h, line 30
simIn.param.PSCD.JERK_MAX_Z_DFLT = single(500); % Default vertical jerk, defined as m/s/s/s in AC_PosControl.h, line 31. Converted to cm/s/s/s in AC_PosControl.cpp, line 318.

simIn.param.PSCD.VEL_MAX_DN = getParamVal(sid, 'PILOT_SPEED_DN'); % Maximum vertical descending velocity in cm/s, defined in parameters.cpp, line 906
simIn.param.PSCD.VEL_MAX_UP = getParamVal(sid, 'PILOT_SPEED_UP'); % Maximum vertical ascending velocity in cm/s, defined in parameters.cpp, line 223
simIn.param.PSCD.ACC_MAX_Z = getParamVal(sid, 'PILOT_ACCEL_Z'); % Maximum vertical acceleration used when pilot is controlling the altitude, parameters.cpp, line 232
simIn.param.PSCD.JERK_MAX_Z = getParamVal(sid, 'PSC_JERK_Z'); % Jerk limit of vertical kinematic path generation in m/s^3. Determines how quickly aircraft changes acceleration target
simIn.param.PSCD.OVERSPEED_GAIN_Z = single(2); % gain controlling rate at which z-axis speed is brought back within SPEED_UP and SPEED_DOWN range, defined in AC_PosControl.h 
%% Position Controller z

% Parameters z position controller
simIn.param.PSCD.P_POS = getParamVal(sid, 'PSC_POSZ_P'); % P gain of the vertical position controller

% Parameters z velocity controller
simIn.param.PSCD.P_VEL = getParamVal(sid, 'PSC_VELZ_P'); % P gain of the vertical velocity controller
simIn.param.PSCD.I_VEL = getParamVal(sid, 'PSC_VELZ_I'); % I gain of the vertical velocity controller
simIn.param.PSCD.D_VEL = getParamVal(sid, 'PSC_VELZ_D'); % D gain of the vertical velocity controller
simIn.param.PSCD.IMAX_VEL = getParamVal(sid, 'PSC_VELZ_IMAX'); % I gain maximu vertical velocity controller
simIn.param.PSCD.FF_VEL = getParamVal(sid, 'PSC_VELZ_FF'); % Feed Forward gain of the vertical velocity controller
simIn.param.PSCD.FLTE_VEL = getParamVal(sid, 'PSC_VELZ_FLTE'); % Cutoff frequency of the error filter (in Hz)
simIn.param.PSCD.FLTD_VEL = getParamVal(sid, 'PSC_VELZ_FLTD'); % Cutoff frequency of the D term filter (in Hz)
simIn.param.PSCD.CTRL_RATIO_INIT = single(2.0); % Initial value of _vel_z_control_ratio, defined in PosControl.h, line 456

% Parameters z acceleration controller
simIn.param.PIDA.P = getParamVal(sid, 'PSC_ACCZ_P'); % P gain of the vertical acceleration controller
simIn.param.PIDA.I = getParamVal(sid, 'PSC_ACCZ_I'); % I gain of the vertical acceleration controller
simIn.param.PIDA.D = getParamVal(sid, 'PSC_ACCZ_D'); % D gain of the vertical acceleration controller
simIn.param.PIDA.IMAX = getParamVal(sid, 'PSC_ACCZ_IMAX'); % I gain maximu vertical acceleration controller
simIn.param.PIDA.FF = getParamVal(sid, 'PSC_ACCZ_FF'); % Feed Forward gain of the vertical acceleration controller
simIn.param.PIDA.FLTE_f = getParamVal(sid, 'PSC_ACCZ_FLTE'); % Cutoff frequency of the error filter (in Hz)
simIn.param.PIDA.FLTD_f = getParamVal(sid, 'PSC_ACCZ_FLTD'); % Cutoff frequency of the D term filter (in Hz)
simIn.param.PIDA.FLTT_f = getParamVal(sid, 'PSC_ACCZ_FLTT'); % Cutoff frequency of the target filter (in Hz)
simIn.param.PIDA.SR_MAX = getParamVal(sid, 'PSC_ACCZ_SMAX'); % Upper limit of the slew rate produced by combined P and D gains
simIn.param.PIDA.SR_TAU = single(1.0); % Slew Rate Tau - not yet available as a parameter of the copter. Set to 1.0 by default (AC_PID.h, line 24).
simIn.param.PIDA.SR_FLT_f = single(25.0);  % Slew Rate lowpass filter cutoff frequency. Defined in SlewLimiter.h, line 14.

% Read signals for z position controller if PSCD message was logged
if isfield(sid, 'PSCD')
    % Signals of PSCD message
    iVec = (sid.PSCD.TimeS >= tStart & sid.PSCD.TimeS <= tEnd);
    simIn.signals.PSCD.Time = single(sid.PSCD.TimeS(iVec)-tStart);
    % Inputs - Multiply by 100 to get cm as used in controller
    % Inversion of signals necessary due to inverted logging
    simIn.signals.PSCD.zPosTar = single(-sid.PSCD.TPD(iVec)*100); % Target z-Position, calculated in set_pos_target_z_from_climb_rate_cm()
    simIn.signals.PSCD.zPosAct = single(-sid.PSCD.PD(iVec)*100); % Actual z-Position, obtained from INAV
    simIn.signals.PSCD.zVelDes = single(-sid.PSCD.DVD(iVec)*100); % Desired z-Velocity, calculated in set_pos_target_z_from_climb_rate_cm()
    simIn.signals.PSCD.zVelAct = single(-sid.PSCD.VD(iVec)*100); % Actual z-Velocity, obtained from INAV
    simIn.signals.PSCD.zVelTar = single(-sid.PSCD.TVD(iVec)*100); % Target z-Velocity, calculated in update_z_controller()
    simIn.signals.PSCD.zAccDes = single(-sid.PSCD.DAD(iVec)*100); % Desired z-Acceleration, calculated in set_pos_target_z_from_climb_rate_cm()
    simIn.signals.PSCD.zAccAct = single(-sid.PSCD.AD(iVec)*100); % Actual z-Acceleration, obtained from AHRS  
    simIn.signals.PSCD.zAccTar = single(-sid.PSCD.TAD(iVec)*100); % Target z-Acceleration, calculated in update_z_controller()  
    
    % Signals of PIDA message
    iVec = (sid.PIDA.TimeS >= tStart & sid.PIDA.TimeS <= tEnd);
    simIn.signals.PIDA.Time = single(sid.PIDA.TimeS(iVec)-tStart);
    % Inputs PID z acceleration
    simIn.signals.PIDA.Err = single(sid.PIDA.Err(iVec)); % Error between target and actual z-acceleration
    simIn.signals.PIDA.Tar = single(sid.PIDA.Tar(iVec));
    simIn.signals.PIDA.Act = single(sid.PIDA.Act(iVec));
    if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal || simIn.param.altCtrlVal
        % Clock of Slew Limiter
        % Use tracked logging time of PID message in validation modes
        simIn.signals.PIDA.ClockDMod = single(1000 * simIn.signals.PIDA.Time);
    else
        simIn.signals.PIDA.ClockDMod = single(1000 * simIn.param.dt * (0:length(simIn.signals.PIDA.Time)-1)');
    end
    % Outputs PID z acceleration
    simIn.signals.PIDA.FF = single(sid.PIDA.FF(iVec));
    simIn.signals.PIDA.P = single(sid.PIDA.P(iVec));
    simIn.signals.PIDA.I = single(sid.PIDA.I(iVec));
    simIn.signals.PIDA.D = single(sid.PIDA.D(iVec));
    simIn.signals.PIDA.DMod = single(sid.PIDA.Dmod(iVec));
    simIn.signals.PIDA.ILimit = single(sid.PIDA.Limit(iVec));
    simIn.signals.PIDA.SRate = single(sid.PIDA.SRate(iVec));    
    
    % Inital inputs
    if simIn.param.rateCtrlVal || simIn.param.attCtrlVal || simIn.param.mdlVal || simIn.param.altCtrlVal
        % Initial actual values - in meter
        simIn.init.PSCD.posAct = single(-sid.PSCD.PD(1));
        simIn.init.PSCD.velAct = single(-sid.PSCD.VD(1));
        simIn.init.PSCD.accAct = single(-sid.PSCD.AD(1));
        
        % Init the position controller to the current position, velocity
        % and acceleration (from AC_PosControl::init_z()) in cm
        simIn.init.PIDA.posTar = single(-sid.PSCD.TPD(1)*100);
        simIn.init.PIDA.velDes = single(-sid.PSCD.DVD(1)*100);
        simIn.init.PIDA.accDes = single(min(max(-sid.PSCD.DAD(1)*100, -simIn.param.PSCD.ACC_MAX_Z_DFLT) ,simIn.param.PSCD.ACC_MAX_Z_DFLT));
        
        simIn.init.PIDA.P = single(sid.PIDA.P(1));
        simIn.init.PIDA.I = single(simIn.signals.CTUN.ThrIn(1) - simIn.signals.CTUN.ThrHov(1)) * 1000 ...
            - simIn.param.PIDA.P * (-sid.PSCD.TAD(1)*100 - (-sid.PSCD.AD(1))*100) ...
            - simIn.param.PIDA.FF * sid.PIDA.Tar(1); % Integrator init. according to AC_PosControl::init_z_controller()
        simIn.init.PIDA.D = single(sid.PIDA.D(1));
        simIn.init.PIDA.TarFilt = single(sid.PIDA.Tar(1));
        simIn.init.PIDA.ErrFilt = single(sid.PIDA.Err(1));
        simIn.init.PIDA.SROut = single(sid.PIDA.SRate(1)); % Initial value of the slew rate determined by the slew limiter. Used for both modifier_slew_rate and output_slew_rate.
        if (simIn.param.PIDA.D ~= 0)
            simIn.init.PIDA.DerFilt = single(sid.PIDA.D(1))/PSC_ACC_Z.param.D; % Divide through D gain to obtain inital D input
        else
            simIn.init.PIDA.DerFilt = single(0);
        end
    else
        % Initial actual values
        simIn.init.PSCD.posAct = single(0);
        simIn.init.PSCD.velAct = single(0);
        simIn.init.PSCD.accAct = single(0);
        
        simIn.init.PIDA.posTar = single(0);
        simIn.init.PIDA.velDes = single(0);
        simIn.init.PIDA.accDes = single(0);
        
        simIn.init.PIDA.P = single(0);
        simIn.init.PIDA.I = single(0);
        simIn.init.PIDA.D = single(0);
        simIn.init.PIDA.TarFilt = single(0);
        simIn.init.PIDA.ErrFilt = single(0);
        simIn.init.PIDA.DerFilt = single(0);
        simIn.init.PIDA.SROut = single(0);
    end
else
    % Set all signals to zero, if PSCD message was not logged and z
    % controller has been deactivated
    iVec = (sid.RATE.TimeS >= tStart & sid.RATE.TimeS <= tEnd);
    simIn.signals.PSCD.Time = single(sid.RATE.TimeS(iVec)-tStart);
    simIn.signals.PSCD.zPosTar = single(zeros(length(simIn.signals.PSCD.Time),1)); % Target z-Position, calculated in set_pos_target_z_from_climb_rate_cm()
    simIn.signals.PSCD.zPosAct = single(zeros(length(simIn.signals.PSCD.Time),1)); % Actual z-Position, obtained from INAV
    simIn.signals.PSCD.zVelDes = single(zeros(length(simIn.signals.PSCD.Time),1)); % Desired z-Velocity, calculated in set_pos_target_z_from_climb_rate_cm()
    simIn.signals.PSCD.zVelAct = single(zeros(length(simIn.signals.PSCD.Time),1)); % Actual z-Velocity, obtained from INAV
    simIn.signals.PSCD.zVelTar = single(zeros(length(simIn.signals.PSCD.Time),1)); % Target z-Velocity, calculated in update_z_controller()
    simIn.signals.PSCD.zAccDes = single(zeros(length(simIn.signals.PSCD.Time),1)); % Desired z-Acceleration, calculated in set_pos_target_z_from_climb_rate_cm()
    simIn.signals.PSCD.zAccAct = single(zeros(length(simIn.signals.PSCD.Time),1)); % Actual z-Acceleration, obtained from AHRS  
    simIn.signals.PSCD.zAccTar = single(zeros(length(simIn.signals.PSCD.Time),1)); % Target z-Acceleration, calculated in update_z_controller()
    
    % Signals of PIDA message
    simIn.signals.PIDA.Time = single(simIn.signals.RATE.Time)-tStart;
    % Inputs PID z acceleration
    simIn.signals.PIDA.Err = single(zeros(length(simIn.signals.PIDA.Time),1)); % Error between target and actual z-acceleration
    simIn.signals.PIDA.ClockDMod = single(zeros(length(simIn.signals.PIDA.Time),1)); % Clock of Slew Limiter

    % Outputs PID z acceleration
    simIn.signals.PIDA.FF = single(zeros(length(simIn.signals.PIDA.Time),1));
    simIn.signals.PIDA.P = single(zeros(length(simIn.signals.PIDA.Time),1));
    simIn.signals.PIDA.I = single(zeros(length(simIn.signals.PIDA.Time),1));
    simIn.signals.PIDA.D = single(zeros(length(simIn.signals.PIDA.Time),1));
    simIn.signals.PIDA.DMod = single(zeros(length(simIn.signals.PIDA.Time),1));
    simIn.signals.PIDA.ILimit = single(zeros(length(simIn.signals.PIDA.Time),1));
    simIn.signals.PIDA.SRate = single(zeros(length(simIn.signals.PIDA.Time),1)); 
    simIn.signals.PIDA.Act = single(zeros(length(simIn.signals.PIDA.Time),1));
    simIn.signals.PIDA.Tar = single(zeros(length(simIn.signals.PIDA.Time),1));
    
    % Initial inputs
    % Initial actual values
    simIn.init.PSCD.posAct = single(0);
    simIn.init.PSCD.velAct = single(0);
    simIn.init.PSCD.accAct = single(0);
    
    simIn.init.PIDA.posTar = single(0);
    simIn.init.PIDA.velDes = single(0);
    simIn.init.PIDA.accDes = single(0);
    
    simIn.init.PIDA.P = single(0);
    simIn.init.PIDA.I = single(0);
    simIn.init.PIDA.D = single(0);
    simIn.init.PIDA.TarFilt = single(0);
    simIn.init.PIDA.ErrFilt = single(0);
    simIn.init.PIDA.DerFilt = single(0);
    simIn.init.PIDA.SROut = single(0);
end


% Real Flight Signals of PID z velocity - only load if log messages
% exists
if isfield(sid, 'PCVZ')
    iVec = (sid.PCVZ.TimeS >= tStart & sid.PCVZ.TimeS <= tEnd);
    simIn.signals.PCVZ.Time = single(sid.PCVZ.TimeS(iVec)-tStart);
    simIn.signals.PCVZ.Err = single(sid.PCVZ.Err(iVec));
    simIn.signals.PCVZ.P = single(sid.PCVZ.P(iVec));
    simIn.signals.PCVZ.I = single(sid.PCVZ.I(iVec));
    simIn.signals.PCVZ.D = single(sid.PCVZ.D(iVec));
    simIn.signals.PCVZ.FF = single(sid.PCVZ.FF(iVec));
else % Set to zero otherwise
    iVec = (sid.RATE.TimeS >= tStart & sid.RATE.TimeS <= tEnd);
    simIn.signals.PCVZ.Time = single(sid.RATE.TimeS(iVec)-tStart);
    simIn.signals.PCVZ.Err = single(zeros(length(simIn.signals.PCVZ.Time),1));
    simIn.signals.PCVZ.P = single(zeros(length(simIn.signals.PCVZ.Time),1));
    simIn.signals.PCVZ.I = single(zeros(length(simIn.signals.PCVZ.Time),1));
    simIn.signals.PCVZ.D = single(zeros(length(simIn.signals.PCVZ.Time),1));
    simIn.signals.PCVZ.FF = single(zeros(length(simIn.signals.PCVZ.Time),1));
end

%% Load identified plant models data
if simIn.param.optCtrl || simIn.param.mdlVal && exist('idModel', 'var')
    % Roll
    if isempty(findobj(idModel, 'axis', 'RLL'))
        simIn.models.RollAxis = idpoly(tf(1, [1 1], simIn.param.dt));
    else
        model = findobj(idModel, 'axis', 'RLL');
        % If more than one model is found, let user choose
        if length(model) > 1 
            mdlIdx = input('More than one model found for the roll axis. Input index of desired model: ');
        else
            mdlIdx = 1;
        end
        simIn.models.RollAxis = idss(c2d(model(mdlIdx).tfModel, simIn.param.dt));
        % Decrease delay to account for the additional feedback delay
        simIn.models.RollAxis.InputDelay = simIn.models.RollAxis.InputDelay-1;
    end
    % Pitch
    if isempty(findobj(idModel, 'axis', 'PIT'))
        simIn.models.PitchAxis = idpoly(tf(1, [1 1], simIn.param.dt));
    else
        model = findobj(idModel, 'axis', 'PIT');
        % If more than one model is found, let user choose
        if length(model) > 1 
            mdlIdx = input('More than one model found for the pitch axis. Input index of desired model: ');
        else
            mdlIdx = 1;
        end
        simIn.models.PitchAxis = idss(c2d(model(mdlIdx).tfModel, simIn.param.dt));
        simIn.models.PitchAxis.InputDelay = simIn.models.PitchAxis.InputDelay-1;
    end
    % Yaw
    if isempty(findobj(idModel, 'axis', 'YAW'))
        simIn.models.YawAxis = idpoly(tf(1, [1 1], simIn.param.dt));
    else
        model = findobj(idModel, 'axis', 'YAW');
        % If more than one model is found, let user choose
        if length(model) > 1 
            mdlIdx = input('More than one model found for the yaw axis. Input index of desired model: ');
        else
            mdlIdx = 1;
        end
        simIn.models.YawAxis = idss(c2d(model(mdlIdx).tfModel, simIn.param.dt));
        simIn.models.YawAxis.InputDelay = simIn.models.YawAxis.InputDelay-1;
    end
    % Vertical Motion - Throttle
    if isempty(findobj(idModel, 'axis', 'THR'))
        simIn.models.RollAxis = idpoly(tf(1, [1 1], simIn.param.dt));
    else
        model = findobj(idModel, 'axis', 'THR');
        % If more than one model is found, let user choose
        if length(model) > 1
            mdlIdx = input('More than one model found for the vertical axis. Input index of desired model: ');
        else
            mdlIdx = 1;
        end
        simIn.models.VerticalAxis = idss(c2d(model(mdlIdx).tfModel, simIn.param.dt));
        % Decrease delay to account for the additional feedback delay
        simIn.models.VerticalAxis.InputDelay = simIn.models.VerticalAxis.InputDelay-1;
    end
else % Create dummy models of optimization is not active
    simIn.models.RollAxis = idpoly(tf(1, [1 1], simIn.param.dt));
    simIn.models.PitchAxis = idpoly(tf(1, [1 1], simIn.param.dt));
    simIn.models.YawAxis = idpoly(tf(1, [1 1], simIn.param.dt));
    simIn.models.VerticalAxis = idpoly(tf(1, [1 1], simIn.param.dt));
end

% Initial values
if simIn.param.mdlVal && mode == 25
    simIn.models.RollInit = single(sid.SIDD.Gx(1)*pi/180);
    simIn.models.PitchInit = single(sid.SIDD.Gy(1)*pi/180);
    simIn.models.YawInit = single(sid.SIDD.Gz(1)*pi/180);
    simIn.models.VerticalOutInit = single(sid.SIDD.Az(1));
    simIn.models.VerticalInInit = single(sid.RATE.AOut(1));
else
    simIn.models.RollInit = single(0);
    simIn.models.PitchInit = single(0);
    simIn.models.YawInit = single(0);
    simIn.models.VerticalOutInit = single(-9.81);
    simIn.models.VerticalInInit = single(sid.RATE.AOut(1));
end

%% Delete all irrelevant data
clear i dist start thr_man rel_msg mode param_names msgs tEnd tStart simMode
clear phase_des runOutsideLoop mdlIdx pitchCh rollCh yawCh iVec Fs axis dt
clear simMode thrCh