% Code to plot the bending modes from the sdl_flexible_shaft_transverse model: flexShaft_eigen block
%% Plot Description:
%
% This figure shows the bending mode shapes of a shaft
% in its two transverse directions with parameters specified in the 
% flexShaft_eigen block of the sdl_flexible_shaft_transverse model.

% Copyright 2018 The MathWorks, Inc.

% Inputs for mode calculation
blockHandle= getSimulinkBlockHandle([system_name '/' eigenmodes_block_name]);

% Reuse figure if it exists, else create new figure
if ~exist('hModes', 'var') || ...
        ~isgraphics(hModes, 'figure')
    hModes = figure('Name', 'Eigenmodes Shaft Bending Vibration');
end
figure(hModes)
clf(hModes)

hX= subplot(2,2,1);
hX.Position= [0.1 0.6 0.6 0.32];
box on
hold on
title('Mode Shapes in X Direction');
ylabel('U_X');
xlabel('Shaft position (m)');

hY= subplot(2,2,3);
hY.Position= [0.1 0.1 0.6 0.32];
box on
hold on
title('Mode Shapes in Y Direction');
ylabel('U_Y');
xlabel('Shaft position (m)');
L= legend();
L.Position= [ 0.66 0.22 0.137 0.043 ];

hCalc= annotation('textbox',...
       'Units', 'normalized', 'Position', [0.38,0.56,0.15,0.06],...
       'FontSize', 16, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle',...
       'String','Calculating...', 'backgroundcolor', 'none', 'LineStyle', 'none',...
        'visible', 'off');

% Choose initial shaft rotation speed [RPM]:
OmegaNom_RPM= 50;

if ~( isvalid(hModes) && isgraphics(hModes) )
    hModes= figure;
end

% User interface

annotation('textbox',...
   'Units', 'normalized', 'Position',[0.7 0.7 0.15 0.06],...
   'FontSize', 8, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle',...
   'String','Shaft rotation speed (RPM):', 'backgroundcolor', 'none', 'LineStyle', 'none');

hText = uicontrol('Style', 'edit', 'String', num2str(OmegaNom_RPM, '%.0f'), ...
    'Units', 'normalized', 'Position', [0.86 0.7 0.1 0.05]);

annotation('textbox',...
   'Units', 'normalized', 'Position',[0.76 0.64 0.16 0.045],...
   'FontSize', 8, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle',...
   'String','Modes to plot:', 'backgroundcolor', 'none', 'LineStyle', 'none');

hText_mode = uicontrol('Style', 'listbox', ...
    'Min', 1, 'Units', 'normalized', 'Position', [0.76 0.38 0.17 0.26]);

hPlot= uicontrol('Style', 'pushbutton', 'String', 'Plot', ...
    'Units', 'normalized', 'Position', [0.73 0.79 0.1 0.05], ... 
    'Callback', {@plot_modes, hX, hY, hText, hText_mode} );

uicontrol('Style', 'pushbutton', 'String', 'Clear', ...
    'Units', 'normalized', 'Position', [0.86 0.79 0.1 0.05], ...
    'Callback', {@clear_axes, hModes});

h_calculate= uicontrol('Style', 'pushbutton', 'String', 'Calculate modes', ...
    'Units', 'normalized', 'Position', [0.73 0.87 0.23 0.05], ...
    'Callback', {@calculate_modes, blockHandle, hCalc, hX, hY, hModes, hText_mode});

if ~exist( 'flex_shaft_mode_props', 'var' )
    calculate_modes([],[],blockHandle, hCalc, hX, hY, hModes, hText_mode);
end


% --- Subfunctions ---
function calculate_modes(~, ~, blockHandle, hCalc, hX, hY, hModes, hText_mode)

    set(hCalc, 'visible', 'on')
    set(hX, 'visible', 'off')
    set(hY, 'visible', 'off')
    clear_axes([], [], hModes)
    drawnow
    
    % Get mode shape arrays with the format [zshaft, mode, shaft_speed]
    [WX, WY, zshaft, eigenfrequencies, shaft_speeds]=...
    sdl_internal.sdl_flexible_shaft_transverse_get_modes_from_block(blockHandle, [], []);

    set(hCalc, 'visible', 'off')
    set(hX, 'visible', 'on')
    set(hY, 'visible', 'on')
    
    % Store properties in a structure
    flex_shaft_mode_props.WX= WX;
    flex_shaft_mode_props.WY= WY;
    flex_shaft_mode_props.zshaft= zshaft;
    flex_shaft_mode_props.modeFreqs= eigenfrequencies;
    flex_shaft_mode_props.shaft_speeds_RPM= shaft_speeds * 60 / (2*pi);

    assignin('base', 'flex_shaft_mode_props', flex_shaft_mode_props);

    % Update listbox of mode #'s
    hText_mode.String= num2str( (1: size(flex_shaft_mode_props.modeFreqs,2) )' );
    hText_mode.Max= size(flex_shaft_mode_props.modeFreqs,2);
    hText_mode.Value= hText_mode.Value( hText_mode.Value <=  hText_mode.Max); % saturate selected values
    
end

function plot_modes( ~, ~, hX, hY, hText, hText_mode)

    OmegaNom_RPM= str2double(hText.String);
    props= evalin('base', 'flex_shaft_mode_props');

    if length( props.shaft_speeds_RPM ) > 1
        % saturate mode properties at input data bounds, as done during simulation
        shaft_speed_lookup= min( max( props.shaft_speeds_RPM(1), OmegaNom_RPM ), props.shaft_speeds_RPM(end) ); 

        % Determine mode shape at the specified rotation speed
        % WX has format [zshaft, mode, shaft_speed]. WX_plot has format [zshaft, mode]
        % interp1 interpolates along columns
        WX_plot= permute( interp1(props.shaft_speeds_RPM, permute(props.WX, [3,1,2]), shaft_speed_lookup), [2,3,1]);
        WY_plot= permute( interp1(props.shaft_speeds_RPM, permute(props.WY, [3,1,2]), shaft_speed_lookup), [2,3,1]);
    else
        % Calculated mode shapes do not depend on rotation speed
        WX_plot= props.WX;
        WY_plot= props.WY;
    end

    if ~isempty(props.modeFreqs)

        if length ( props.shaft_speeds_RPM ) > 1
            mode_frequency_plot= interp1(props.shaft_speeds_RPM, props.modeFreqs, shaft_speed_lookup);
        else
            mode_frequency_plot = props.modeFreqs;
        end

        for i= hText_mode.Value
            line_label = [ 'w' num2str(i) '= ',...
            num2str( mode_frequency_plot(i)*60/(2*pi), '%.0f'  ) ' RPM. ',...
            'Shaft speed= ' num2str( OmegaNom_RPM ) ' RPM' ];
            plot(hX, props.zshaft, WX_plot(:,i), 'DisplayName', line_label );
        end
        for i= hText_mode.Value
            line_label = [ 'w' num2str(i) '= ',...
            num2str( mode_frequency_plot(i)*60/(2*pi), '%.0f'  ) ' RPM. ',...
            'Shaft speed= ' num2str( OmegaNom_RPM ) ' RPM' ];
            plot(hY, props.zshaft, WY_plot(:,i), 'DisplayName', line_label );
        end
    end
   
end

function clear_axes(~, ~, hModes)
    allaxes = findall(hModes, 'type', 'axes');
    hX= allaxes(1);
    hY= allaxes(2);

    delete(hX.Children)
    delete(hY.Children)
    
end

