% Code to plot simulation results from sh_testrig_prv_direct_acting
%% Plot Description:
%
% The plots below show the behavior of the direct-acting pressure relief
% valve.  As the load orifice closes, pressure in the system rises.  When
% the system pressure reaches the pressure setting for the pressure relief
% valve, the valve opens which diverts flow away from the orifice and
% allows the system to maintain the set pressure of the pressure relief
% valve.

% Copyright 2015 The MathWorks, Inc.

% Generate simulation results if they don't exist
if ~exist('simlog_sh_testrig_prv_direct_acting', 'var')
    sim('sh_testrig_prv_direct_acting')
end

% Reuse figure if it exists, else create new figure
if ~exist('h1_sh_testrig_prv_direct_acting', 'var') || ...
        ~isgraphics(h1_sh_testrig_prv_direct_acting, 'figure')
    h1_sh_testrig_prv_direct_acting = figure('Name', 'sh_testrig_prv_direct_acting');
end
figure(h1_sh_testrig_prv_direct_acting)
clf(h1_sh_testrig_prv_direct_acting)

% Get simulation results
simlog_t = simlog_sh_testrig_prv_direct_acting.Orifice.flow_rate.series.time;
simlog_orifq = simlog_sh_testrig_prv_direct_acting.Orifice.flow_rate.series.values('m^3/s');
simlog_prvq = simlog_sh_testrig_prv_direct_acting.Pressure_Relief_Valve_Direct_Acting.Orifice.flow_rate.series.values('m^3/s');
simlog_orifS = simlog_sh_testrig_prv_direct_acting.Orifice.S.series.values('m');
simlog_prvp = simlog_sh_testrig_prv_direct_acting.Pressure_Relief_Valve_Direct_Acting.Orifice.pressure_drop.series.values('Pa');

temp_orifMaxO = str2num(get_param([bdroot '/Orifice'],'opening_max'));
temp_orifMaxA = str2num(get_param([bdroot '/Orifice'],'area_max'));
temp_orifIniO = str2num(get_param([bdroot '/Orifice'],'x_0'));
temp_orifori = str2num(get_param([bdroot '/Orifice'],'or'));
if(temp_orifori==2)
    temp_orifori = -1;
end
    
% Plot results
simlog_handles(1) = subplot(2, 1, 1);
plot(simlog_t, simlog_orifq, 'LineWidth', 1)
hold on
plot(simlog_t, simlog_prvq, 'LineWidth', 1)
hold off
grid on
title('Flow Rates')
ylabel('Flow Rates (m^3)')
legend({'Orifice','Pressure Relief Valve'},'Location','Best');

simlog_handles(2) = subplot(2, 1, 2);
[handles_Ax, handles_Line1, handles_Line2] = plotyy(...
    simlog_t, max((temp_orifIniO+simlog_orifS*temp_orifori)/temp_orifMaxO*temp_orifMaxA,0), ...
    simlog_t, simlog_prvp);
grid on
handles_Line1.LineWidth = 1;
handles_Line2.LineWidth = 1;
ylabel(handles_Ax(1), 'Opening (m^2)') % left y-axis
ylabel(handles_Ax(2), 'Pressure (Pa)') % right y-axis

title('Orifice Opening and Pressure Drop')
xlabel('Time (s)')

linkaxes([handles_Ax simlog_handles(1)], 'x')

% Remove temporary variables
clear simlog_t simlog_handles
clear simlog_orifS simlog_orifp simlog_orifq simlog_prvp simlog_prvq 
clear temp_orifIniO temp_orifori temp_orifMaxA temp_orifMaxO
clear handles_Ax handles_Line1 handles_Line2 