% Code to animate transverse vibrations from sdl_flexible_shaft_transverse model: flexShaft_eigen block
%% Plot Description:
%
% This figure shows an animation of the bending deflection of a rotating shaft
% in its two transverse directions over time. The transverse deflection is modeled by eigenmodes.
% Positive shaft rotation in this animation corresponds to counter-clockwise rotation about the
% Z-axis, which points from the shaft B node toward the F node. 
% The red circles correspond to concentrated inertias on the shaft.

% Copyright 2018-2020 The MathWorks, Inc.

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

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

eigenmodes_blockHandle= getSimulinkBlockHandle([system_name '/' eigenmodes_block_name]);

shaft_length= sdl_internal.sdl_flexible_shaft_transverse_get_dialog_box_value(eigenmodes_blockHandle, 'length', 'm');
numPts= 100; %number of points to plot along shaft axis
z= linspace(0, shaft_length, numPts); %shaft locations to plot in animation

% Get mode shapes at each reference shaft speed.
% mode shape arrays have the format [shaft_speed, z_shaft, mode]
[UX_tab, UY_tab,~, ~, shaft_speed_tab]= sdl_internal.sdl_flexible_shaft_transverse_get_modes_from_block(eigenmodes_blockHandle, z, []);
numSpeeds_computed= size(shaft_speed_tab, 2);
numModes_computed= size(UX_tab,2);

UX_tab= UX_tab(:,1:numModes_computed,:); 
UY_tab= UY_tab(:,1:numModes_computed,:);

% Simulation data
t = eval(['simlog_sdl_flexible_shaft_transverse.' eigenmodes_block_name '.fsBendEigen.eta.series.time']);
shaft_speed= eval(['simlog_sdl_flexible_shaft_transverse.' eigenmodes_block_name '.B.w.series.values(''rad/s'')']);

%  Deta, mode time-dependent coefficient displacement of real and imaginary components of the modes
% format: [real components; imaginary components]
Deta= eval(['simlog_sdl_flexible_shaft_transverse.' eigenmodes_block_name '.fsBendEigen.eta.series.values(''1'')']);


% Animation frame data based on geometry and simulation data
Dx= zeros(numel(t), numel(z));
Dy= zeros(numel(t), numel(z));

for j= 1:numel(t)
    
  % time components of the modes
  eta= Deta(j,:); 

  %expand time component to every position along shaft
  etaMat= repmat(eta, numel(z),1 ); 

  if numSpeeds_computed >1
     shaft_speed_lookup= min( max( shaft_speed_tab(1),shaft_speed(j) ), shaft_speed_tab(end) ); % saturate mode properties at input data bounds
     
     % 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
     UX_inst= permute( interp1(shaft_speed_tab, permute(UX_tab, [3,1,2]), shaft_speed_lookup), [2,3,1]);
     UY_inst= permute( interp1(shaft_speed_tab, permute(UY_tab, [3,1,2]), shaft_speed_lookup), [2,3,1]);
  else % constant mode properties
     UX_inst= UX_tab;
     UY_inst= UY_tab; 
  end
  
  %Sum column-wise to get final instantaneous displacement of shaft at each position.
  Dx(j,:)=  sum( real( UX_inst'.*etaMat') , 1); 
  Dy(j,:)=  sum( real( UY_inst'.*etaMat') , 1);

end

D= Dx + 1i.*Dy;

% Rigid masses or "disks" that are attached to shaft
no_Disk= str2double(get_param(eigenmodes_blockHandle, 'noDisk')); % Flag indicating if shaft has disks
z_Disk= sdl_internal.sdl_flexible_shaft_transverse_get_dialog_box_value(eigenmodes_blockHandle, 'zDisk', 'm');

if no_Disk ~= 1
    %Disk deflection versus time
    d_Disk= interp1(z', D', z_Disk);
    d_Disk= d_Disk';
else
    d_Disk= [];
    z_Disk= [];
end

sdl_flexible_shaft_transverse_animate(t, z, D, z_Disk, d_Disk, h_sdl_flexible_shaft_transverse_eigenmodes);

% Remove temporary variables
clear eigenmodes_blockHandle numPts z
clear t Vx Vy DetaX DetaY
clear UX_tab UY_tab UX_inst UY_inst numModes_computed
clear shaft_speed_tab shaft_speed numSpeeds_computed
clear Dx Dy D Deta etaX eta etaMat
clear no_Disk d_Disk z_Disk j
clear shaft_length number_modes

