%% Analyzing a Ball Impacting a Lens
% This example demonstrates basic approaches using various DSP
% functionality to correlate impact results from a physical experiment 
% and a FEA simulation for a steel ball-bearing hitting the plastic 
% lens of a cell phone. 

%% Overview
%{
Key concepts and features demonstrated include: 
 o Utilizing Units with all the calculations.
 o Using a combination of DSP techniques to characterize frequency
   content and then filter out noise to achieve significant 
   improvement in data interpretation and correlation. 

Tips:
 o Consider running the example "By Section". It will help you walk 
   through the example and its results.
%}

%% Set-up
%{
This section makes a few initial settings for the example.
 
Additional Details: 
The settings below include: 
 o Popup asking to clean MATLAB session before start of example.
 o Define a handy "adv" variable so that we can later use it to 
   quickly access any Kornucopia function's ADV options with ease.
 o Setting a specific Units Preference for this example. 
 o Define some units variables that will be used throughout the 
   example.
 o Turn off a warning regarding Props not being inheritted during
   certain math operations.
 o Create a tabbed figure with the main tab titles that will be used for
   the entire analysis
%}

k_cleanup();
adv = k_adv();

k_unitsPreferenceActivate('mm_N_ms'); 
k_unitsVariables('mm, ms, N, G, kHz');
k_set('inheritPropsAutoWarning', 'off');

figH = figure('pos', k_figZoom(2));
tabTitles = {
    'Description and Raw Data'
    'Initial Comparision of Accel'
    'Data Clean-up and Fourier Assessments of Accel'
    'After LP Filtering Accel'
    'Assessing Displacements and Adjusting Stiffness'   
    'Evaluating Shock Response'
    };
tabStuffMain = k_figTabsCreate(tabTitles, figH);

%% Loading data into MATLAB session
%{
This section loads and displays the raw data to be analyzed along with a 
png image file that describes the impact problem being analyzed. 
%}

dataDir = '..\data';

% Load the raw FEA and experimental data 
fileNames = {
    'experiment.txt'
    'fea.txt'
    }; 
delim = ',';
startStr = [];
rs1 = [];
colHeadStr = 'Time';
rs2 = 0;
endStr = [];
rs3 = [];

fName = fullfile(dataDir, fileNames{1});
rawExper = k_readText(fName,delim,startStr,rs1,colHeadStr,rs2,endStr,rs3);

fName = fullfile(dataDir, fileNames{2});
rawFea = k_readText(fName,delim,startStr,rs1,colHeadStr,rs2,endStr,rs3);

% Define the image file and fully qualified file name (with path).
pngImage_fileName = 'ball_lens_impact.png';
png_fullFileName = fullfile(dataDir, pngImage_fileName);

% Present the data
currentTabH = tabStuffMain.TabsH(1);
k_figTabsClear(currentTabH);

% Create "subtabs"
titles = {'Picture of Problem', 'Raw Data'};
tabStuffLocal = k_figTabsCreate(titles, currentTabH);

% Picture of problem
currentTabH = tabStuffLocal.TabsH(1);
k_figTabsClear(currentTabH);

k_picDisplay(png_fullFileName, 'center', ...
    'parent', currentTabH);

k_figTabsDisplay(currentTabH);

% The raw data
currentTabH = tabStuffLocal.TabsH(2);
k_figTabsClear(currentTabH);

aH = subplot(1,2,1, 'parent', currentTabH);
k_varViewer(rawExper, 'parent', aH)

aH = subplot(1,2,2);
k_varViewer(rawFea, 'parent', aH)
k_figTabsDisplay(currentTabH);

%% Plot the raw acceleration data
%{
This section plots the acceleration data that was loaded. 

Before we plot the data, we need to do a few house-keeping items and
then we plot the data.

Additional Details: 
 o Both data sets are converted to the active Units Preference so that
   quantities appear in our preferred units. 
 o Additionally, the Accel cols are converted to units of 1000*G
 o We will be doing DSP on the data, so the FEA data is regularized
   since we suspect it has non-constant time increments due to the
   Explicit Dynamics solution technique used to generate the data. 
 o The sampling frequency of the regularized FEA data is much higher 
   than the experimental data, so we decimate the FEA data to a
   sampling rate that is similar to the experimental data. 
    - Because the data starts with a sharp pulse or spike very early
      in time, we need to "zero pad" the beginning of the FEA data via 
      the function k_extend. This is to avoid DSP related distortions 
      that will occur when bi-directionally filtering a spike that 
      occurs at the very beginning of a data set. 
%}

rawFea = rawFea.convert();
rawExper = rawExper.convert();

rawFea('Accel') = rawFea('Accel').convert('1000*G');
rawExper('Accel') = rawExper('Accel').convert('1000*G');

% Regularizing the FEA data
regFea = k_regularize(rawFea);

% Extend the beginning of the FEA data with zeros.
extendedFea = k_extend(regFea, [], 1/6, 'zero', [], []);
extendedFea.Props.Description = 'Extended FEA Data';

% Assessing sampling frequency
fs_Exper = k_samplingFreq(rawExper('Time'));
fs_regFea = k_samplingFreq(extendedFea('Time'));

% Decimating the FEA data
decFactor = fs_regFea / fs_Exper;
decFea = k_decimate(extendedFea, [], decFactor);
decFea.Props.Description = 'Extended & Decimated FEA Data';
fs_decFea = k_samplingFreq(decFea('Time'));

% Plotting the results
currentTabH = tabStuffMain.TabsH(2);
k_figTabsClear(currentTabH);

txt = {'No experimental data available that is'
       'comparable to high sampling rate FEA data'}; 
k_displayOnFigure(0.25, 0.75, txt, ...
    'parent', currentTabH, ...
    'box', 'on', ...
    'alignHor', 'center')

subplot(2,2,2, 'parent', currentTabH);
k_plot(regFea, [], ...
    'depend', 'Accel', ...
    'title', {'user', 'Regularized FEA data'})
k_displayOnFigure(0.5, 0.2, fs_regFea, 'removeVar', 'off')


pltADVs = {'xLim', [-2, 10]*ms, ...
           'yLim', [-30, 10]*1000*G};

subplot(2,2,3)
k_plot(rawExper, [], ...
    pltADVs, ...
    'depend', 'Accel', ...
    'title', {'user', 'Raw Experimental data'})    
k_displayOnFigure(0.5, 0.2, fs_Exper, 'removeVar', 'off')


subplot(2,2,4)
k_plot(decFea, [], ...
    pltADVs, ...
    'depend', 'Accel', ...
    'title', {'user', 'Extended & Decimated FEA data'})
k_displayOnFigure(0.5, 0.2, fs_decFea, 'removeVar', 'off')

k_figTabsDisplay(currentTabH);

%% Data Clean-up, then perform DFS Analysis
%{
In this section we do a touch more house-keeping to get the experimental
data set and FEA data sets ready for analysis, and then we use the function
k_DFS to perform Fourier analysis of the data.

Additional Details: 
 o We trim the experiment to end at a similar time as the FEA
 o For ease of computation, the two data sets are placed into a cell
   array and then a simple for-loop is used for the calculations.
 o In the DFS calculation, a transient Hann window is specified. The
   transient ADV option is used so that the window does not squelch
   the impact spike at the beginning. You are encouraged to explore 
   the results using other supported Kornucopia windows as well as 
   other setting values for the transient ADV option for k_DFS. 
 o The second DFS subplot is used to make a zoomed-in plot of the Spectra.

%}

% Trim the experimental data to be similar time duration as FEA
idxMax = k_indexFind(rawExper('Time'),'cutBefore',decFea(end,'Time'));
trimmedExper = rawExper(1:idxMax, :);
trimmedExper.Props.Description = 'Trimmed Experimental Data';

% Stack the data sets into cell array and then loop to compute DFS 
dataSets = {trimmedExper; decFea};
descriptions = {'Experiment', 'FEA'};
nSets = numel(dataSets);

dfs = cell(nSets, 1);
for z = 1 : nSets
    tmp = dataSets{z};
    tmp = k_DFS(tmp('Time, Accel'), [], 'hann', 'mag', ...
                'transient', 0.5);
    tmp.Props.Description = descriptions{z};
    dfs{z,1} = tmp;
end


% Plot the results
currentTabH = tabStuffMain.TabsH(3);
k_figTabsClear(currentTabH);
      
tmpH = cell(2,1);       
subplot(2,2,1, 'parent', currentTabH)
tmpH{1} = k_plot(trimmedExper, [], ...
    'depend', 'Accel', ...
    'title', 'description', ...
    'lineColor', 'blue');

subplot(2,2,3)
tmpH{2} = k_plot(decFea, [], ...
    'depend', 'Accel', ...
    'title', 'description', ...
    'lineColor', 'red');

k_plotMatchLimits(tmpH, 'xy'); % This makes all the plots in tmpH have same limits

subplot(2,2,2)
k_plot(dfs, [], ...
    'title', {'user', 'DFS (full spectrum)'}, ...
    'lineColor', 'blue, red', ...
    'legendLoc', 'NorthEast')    

subplot(2,2,4)
k_plot(dfs, [], ...
    'title', {'user', 'DFS (lower freq''s of spectrum)'}, ...
    'lineColor', 'blue, red', ...    
    'legendLoc', 'NorthEast', ...
    'xlim', [0, 25]*kHz)


k_figTabsDisplay(currentTabH);

%% Apply LP filtering
%{
From the previous DFS analysis, it would appear that below 5*kHz there
is very good correlation between the FEA and the experiment.

You are encouraged to try other cutoff frequencies (higher and lower)
than the one demonstrated below. 

The LP filter applied is a bi-directional filter, so there are no time
delays in the filtered results.
%}

% Cutoff freq for the filter
fc = 5*kHz;

% Sampling frequency for each data set is the same since the FEA was 
% decimated to same sampling rate as the Experimental data
fs = fs_Exper;

% Compute Filter Coefficients and apply to the data
C_filt = k_cIIR(fs, fc, 'lp', 'butter');  % Filter coefficients

dataSets_f = cell(nSets,1);
for z = 1 : nSets            
    % Applying the filter coefficients to the data sets
    dataSets_f{z} = k_filter(dataSets{z}, [], C_filt, 'bi-dir');
    dataSets_f{z}.Props.Description = descriptions{z};
end

% Calculating the filter response is optional. It is not required 
% to perform filtering. It is helpful to see what the filter will
% do.
% A vector of frequencies used to make a FRF plot to show the 
% filter's response. 
f = k_vectorByPoints(0*kHz, 2*fc, 100); % A vector of freqs for FRF plot
f.Props.ColNames = 'Frequency';
frf = k_FRF(fs, f, C_filt, 'bi-dir', 'mag'); % The actual FRF calc
frf.Props.Description = 'Filter''s FRF';
    

% Plotting the results
currentTabH = tabStuffMain.TabsH(4);
k_figTabsClear(currentTabH);

subplot(2,2,1, 'parent', currentTabH)
k_plot(f, frf, ...
    'yLim', [0,1.2], ...
    'legendSingleCurve', 'on', ...
    'legendLoc', 'NorthEast')

txt = {'FRF for the two data'
       'sets essentially identical'};
k_displayOnFigure(0.75, 0.7, txt)

subplot(2,2,2)
k_plot(dataSets_f, [], ...
    'depend', 'Accel', ...
    'legendLoc', 'SouthEast',...
    'title', {'user', 'Excellent Correlation via Filtering'})

k_displayOnFigure(0.5, 0.98, fc, ...
    'parent', currentTabH, ...
    'alignVert', 'top', ...
    'box', 'on', ...
    'removeVarName', 'off');

k_figTabsDisplay(currentTabH);

%% Evaluating transient displacement results for the lens
%{
This next portion of the example assesses the displacement response of
the lens due to the impact. 

While the FEA data already has transient displacement results, the 
experimental data does not. We need to compute it from the 
experiment's acceleration vs. time transient data. 

Note that a warning may appear related to Props Inheritance. The
warning can be ignored in this particular case. 
%}

% Computing transient displacements via integration of Accel data.
tmp = dataSets_f{1}; % This is the filtered experimental data
v_exper = k_integrate(tmp('Time'), tmp('Accel'));
u_exper = k_integrate(tmp('Time'), v_exper);
u_exper.Props.ColNames = 'Disp';

% Including the calculated displacement to previous filtered 
% experimental data.
dataSets_f{1}(:,3) = u_exper;

%{
The plot that will be created below will show that the max
displacement of the FEA data is larger than the experimental
results. We assume this difference is caused by a stiffness 
difference in the FEA model relative to the experiment, likely due to
strain-rate effects causing a higher physical stiffness during impact.

Note: 
dataSets_f{1} is the LP filtered experimental data
dataSets_f{2} is the LP filtered FEA data
%}

K_ratio = min(dataSets_f{2}('Disp')) / min(dataSets_f{1}('Disp'));
K_ratio.Props = []; % Remove any Props

% Based on simple spring/mass model, the frequency should scale via:
f_ratio = sqrt(K_ratio);

% Apply the stiffness and frequency scaling factors to the FEA to
% adjust the displacement results. 
u_FEA_adjusted = (1/K_ratio) * dataSets_f{2}('Disp');
t_FEA_adjusted = (1/f_ratio) * dataSets_f{2}('Time');
adjustedFea_Disp = [t_FEA_adjusted, u_FEA_adjusted];
adjustedFea_Disp.Props.ColNames = 'Time, Disp';
adjustedFea_Disp.Props.Description = 'Stiffness & time warp adjusted FEA';


% Plotting the results
currentTabH = tabStuffMain.TabsH(5);
k_figTabsClear(currentTabH);

subplot(3,2,[1,3], 'parent', currentTabH)
toPlot = dataSets_f;
k_plot(toPlot, [], ...
    'depend', 'Disp', ...
    'legendLoc', 'SouthOutside')

k_displayOnFigure(0.4, 0.9, K_ratio, 'removeVar', 'off')
k_displayOnFigure(0.4, 0.8, f_ratio, 'removeVar', 'off')

subplot(3,2,[2,4])
toPlot = {
    dataSets_f{1}
    adjustedFea_Disp
    };
k_plot(toPlot, [], ...
    'depend', 'Disp', ...
    'legendLoc', 'SouthOutside',...
    'title', {'user', 'Adjusted FEA Matches Well'})

txt = {'The excellent match suggests'
       'the FEA model''s material law'
       'needs to be stiffened by'
       k_stringCreate({'K_ratio = ', K_ratio, '.'})};
   
k_displayOnFigure(0.75, 0.25, txt, ...
    'parent', currentTabH, ...
    'box', 'on');

k_displayOnFigure(0.5, 0.98, fc, ...
    'parent', currentTabH, ...
    'alignVert', 'top', ...
    'box', 'on', ...
    'removeVarName', 'off');

k_figTabsDisplay(currentTabH);

%% Assessing Shock Response via PVSS
%{
Below we compute shock response for both the experiment and the FEA (for
both before and after making the stiffness adjustement correction).
%}

% First we re-calculate the accel of the adjustedFea.
adjustedFea_vel = k_derivative(adjustedFea_Disp);
adjustedFea_vel.Props.ColNames{2} = 'Vel';

adjustedFea_accel = k_derivative(adjustedFea_vel);
adjustedFea_accel.Props.ColNames{2} = 'Accel';

toProcess = [
    dataSets_f
    {adjustedFea_accel}
    ];

% For SRS calcs via PVSS
damp = 0.05; 

PVSS = cell(3,1);
for z = 1 : 3
    PVSS{z} = k_SRS(toProcess{z}('Time, Accel'), [], '10*Hz', 'fnMax', damp, 'PVSS', ...
        adv.k_SRS.modifyColNames.replace);
end

% Plot results
currentTabH = tabStuffMain.TabsH(6);
k_figTabsClear(currentTabH);

subplot(2,1,1, 'parent', currentTabH)
k_plot(toProcess, [], ...
    'depend', 'Accel', ...
    'legendLoc', 'best')

k_displayOnFigure(0.98, 0.98, fc, ...
    'parent', currentTabH, ...
    'alignVert', 'top', ...
    'alignHor', 'right', ...
    'box', 'on', ...
    'removeVarName', 'off');

subplot(2,1,2, 'parent', currentTabH)
k_plot(PVSS, [], ...
    'legendLoc', 'best', ...
    'scale', 'fourCP')

k_figTabsDisplay(currentTabH);

%% Options to save figure windows

k_exampleSetup('figSavingOptions');




