%% Simulate Model
sim(model)

%% Get simulation results
simlog01 = simlog_ShaftWithTorsionalFlexibility;
simlog_t = simlog01.Clutch.P.series.time;
simlog_pClutch = simlog01.Clutch.P.series.values('Pa');
simlog_wMotor = simlog01.Shaft_Motor.F.w.series.values('rad/s');
simlog_wLoad = simlog01.Shaft_Load.F.w.series.values('rad/s');

%% Plot results
X1 = simlog_t;
YMatrix1 = [simlog_wMotor simlog_wLoad];
Y1 = simlog_pClutch;

% Create figure
figure1 = figure('Name','ShaftWithTorsionalFlexibility',...
    'OuterPosition',[107 336 733 822]);

% Create axes
axes1 = axes('Parent',figure1,...
    'Position',[0.13 0.709 0.775 0.216]);
hold(axes1,'on');

% Activate the left side of the axes
yyaxis(axes1,'left');

% Create multiple lines using matrix input to plot
plot1 = plot(X1,YMatrix1,'LineWidth',2);
set(plot1(1),'DisplayName','Motor Shaft', 'Color', 'k');
set(plot1(2),'DisplayName','Load Shaft', 'Color', 'b');

% Create ylabel
ylabel('Speed (rad/s)');

% Comment/uncomment the following line to remove/
%  preserve the Y-limits of the axes
ylim(axes1,[-100 250]);

% Set the remaining axes properties
set(axes1,'YColor',[0 0.447 0.741]);

% Activate the right side of the axes
yyaxis(axes1,'right');

% Create plot
plot(X1,Y1,'DisplayName','Clutch Pressure','LineWidth',2, 'Color', 'r');

% Create ylabel
ylabel('Pressure (Pa)');

% Comment/uncomment the following line to remove/
%  preserve the Y-limits of the axes
ylim(axes1,[0 2000000]);

% Set the remaining axes properties
set(axes1,'YColor',[0.85 0.325 0.098]);
% Create xlabel
xlabel('Time (s)');

% Create title
title('Shaft Speeds and Clutch Pressure');

% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[0 0.2);

% Set the remaining axes properties
set(axes1,'LineStyleOrderIndex',2);
box(axes1,'on');
grid(axes1,'on');

% Create legend
legend1 = legend(axes1,'show');
% set(legend1,'Location','best');
set(legend1,...
    'Position',[0.6944 0.9125 0.1982 0.07270]);

%% Replace Rotational Damper from Foundation Library with
%  Faultable Rotational Damper from the Simscape Driveline Library

% Define Unfaultable Damper Block 
foundationDamper = [model,'/Damping Bearing'];

% Get Unfaultable Damper Block Damping Coefficient
dampingCoefficient = get_param(foundationDamper,'D');

% Get Unfaultable Damper Block Position
damperPosition = get_param(foundationDamper, 'Position');

% Delete Unfaultable Damper Block
delete_block(foundationDamper)

% Add Faultable Damper Block
faultableDamperLib = ...
    'sdl_lib/Couplings & Drives/Springs & Dampers/Rotational Damper';
faultableDamperPath = [model,'/Faultable Damper'];
add_block(faultableDamperLib,faultableDamperPath,...
    'Position',damperPosition,...
    'NamePlacement', 'alternate',...
    'D',dampingCoefficient)

%% Define Faultable Damper Block Parameters
underDamp = '10';
overDamp = '150';

faultTime = '.06';

aMaxDef = '100';
aMaxHigh = '400';
aMaxLow = '50';

shockNMaxDef = '1';
shockNMaxHigh = '5';
shockNMaxLow = '2';

disabled = '0';
enabled = '1';

none = '1';
warning = '2';
error = '3';

%% Parameterize a Timed Fault
set_param(faultableDamperPath,...
     'enable_faults', enabled,...       % Enable faults
     'b_fault', underDamp,...           % Faulted damping coefficient 
     'temporal_fault', enabled,...      % Enable temporal fault trigger
     'fault_time', faultTime,....       % Simulation time for fault event
     'shock_fault', disabled,...        % Disable behavioral fault trigger
     'acceleration_limit', aMaxDef,...  % Maximum permissible acceleration
     'shock_limit', shockNMaxDef ,...   % Maximum number of shocks
     'report_fault', warning)           % Reporting when fault occurs

%% Simulate
sim(model)

%% Get Simulation Results
simlog02 = simlog_ShaftWithTorsionalFlexibility;
simlog_t02 = simlog02.Clutch.P.series.time;
simlog_pClutch02 = simlog02.Clutch.P.series.values('Pa');
simlog_wMotor02 = simlog02.Shaft_Motor.F.w.series.values('rad/s');
simlog_wLoad02 = simlog02.Shaft_Load.F.w.series.values('rad/s');

%% Plot Results
X2 = simlog_t02;
YMatrix2 = [simlog_wMotor02 simlog_wLoad02];
Y2 = simlog_pClutch02;

% Create axes
axes2 = axes('Parent',figure1,...
    'Position',[0.125 0.4 0.775 0.216]);
hold(axes2,'on');

% Activate the left side of the axes
yyaxis(axes2,'left');

% Create multiple lines using matrix input to plot
plot1 = plot(X2,YMatrix2,'LineWidth',2);
set(plot1(1),'DisplayName','Motor Shaft', 'Color', 'k');
set(plot1(2),'DisplayName','Load Shaft', 'Color', 'b');

% Create ylabel
ylabel('Speed (rad/s)');

% Comment/uncomment the following line to remove/
%  preserve the Y-limits of the axes
ylim(axes2,[-100 250]);

% Set the remaining axes properties
set(axes2,'YColor',[0 0.447 0.741]);

% Activate the right side of the axes
yyaxis(axes2,'right');
% Create plot
plot(X2,Y2,'DisplayName','Clutch Pressure','LineWidth',2, 'Color', 'r');

% Create ylabel
ylabel('Pressure (Pa)');

% Comment/uncomment the following line to remove/
%  preserve the Y-limits of the axes
ylim(axes2,[0 2000000]);

% Set the remaining axes properties
set(axes2,'YColor',[0.85 0.325 0.098]);

% Create title
title('Time-Faulted Speeds and Clutch Pressure');

% Create xlabel
xlabel('Time (s)');

% Uncomment the following line to 
% preserve the X-limits of the axes
% xlim(axes1,[0 0.2]);

% Set the remaining axes properties
box(axes2,'on');
grid(axes2,'on');
set(axes2,'LineStyleOrderIndex',2);

%% Paramaterize a Shock Fault
set_param(faultableDamperPath,...
     'enable_faults', enabled,...       % Enable faults
     'b_fault', overDamp,...           % Faulted damping coefficient 
     'temporal_fault', disabled,...     % Disable temporal fault trigger
     'fault_time', faultTime,....       % Simulation time for fault event
     'shock_fault', enabled,...         % Enable behavioral fault trigger
     'acceleration_limit', aMaxHigh,...  % Maximum permissible acceleration
     'shock_limit', shockNMaxLow ,...   % Maximum number of shocks
     'report_fault', warning)           % Reporting when fault occurs

%% Simulate
sim(model)

%% Get Simulation Results
simlog03 = simlog_ShaftWithTorsionalFlexibility;
simlog_t03 = simlog03.Clutch.P.series.time;
simlog_pClutch03 = simlog03.Clutch.P.series.values('Pa');
simlog_wMotor03 = simlog03.Shaft_Motor.F.w.series.values('rad/s');
simlog_wLoad03 = simlog03.Shaft_Load.F.w.series.values('rad/s');
 
%% Plot Results
X3 = simlog_t03;
YMatrix3 = [simlog_wMotor03 simlog_wLoad03];
Y3 = simlog_pClutch03;

% Create axes
axes3 = axes('Parent',figure1,...
    'Position',[0.125 0.08 0.775 0.216]);
hold(axes3,'on');

% Activate the left side of the axes
yyaxis(axes3,'left');

% Create multiple lines using matrix input to plot
plot1 = plot(X3,YMatrix3,'LineWidth',2);
set(plot1(1),'DisplayName','Motor Shaft', 'Color', 'k');
set(plot1(2),'DisplayName','Load Shaft', 'Color', 'b');

% Create ylabel
ylabel('Speed (rad/s)');

% Comment/uncomment the following line to 
% remove/preserve the Y-limits of the axes
ylim(axes3,[-100 250]);

% Set the remaining axes properties
set(axes3,'YColor',[0 0.447 0.741]);

% Activate the right side of the axes
yyaxis(axes3,'right');

% Create plot
plot(X3,Y3,'DisplayName','Clutch Pressure','LineWidth',2, 'Color', 'r');

% Create ylabel
ylabel('Pressure (Pa)');

% Comment/uncomment the following line to 
% remove/preserve the Y-limits of the axes
ylim(axes3,[0 2000000]);

% Set the remaining axes properties
set(axes3,'YColor',[0.85 0.325 0.098]);

% Create title
title('Shock-Faulted Shaft Speeds and Clutch Pressure');

% Create xlabel
xlabel('Time (s)');

% Uncomment the following line to 
% preserve the X-limits of the axes
% xlim(axes1,[0 0.2);

% Set the remaining axes properties
box(axes3,'on');
grid(axes3,'on');
set(axes3,'LineStyleOrderIndex',2);