clear all
clc
set(0,'defaultfigurecolor',[1 1 1])

%% INFO

% Read and process friction data

% Figure concering effect of heating time for ESAFORM 2022 conference paper

% Matlab version: R2021a


%% initialisation
mod_coeff = 0.90; %threshold for modulus fit, end; max difference local and average modulus
serie = 9; % measurement series
range_meas = [25 26 27 22 23 24 20 21];  % measurement within series, tool-ply and ply-ply tests
t_wait = [10 5 2 10 10 10 5 5]; % heating times in minutes
t_wait_prev = 0; % switch for legend and plotting

%% Go through the measurements

figure(1), clf(1), hold on

% go through measurements within series
for meas_i = 1:length(range_meas)

    meas = range_meas(meas_i);

    %% Load data       
    file_find = dir(['S_c_pekk_','*_',num2str(serie),'-',num2str(meas),'.*']);
    file_name = file_find.name;        
    data = readmatrix(file_name,'NumHeaderLines',1);
    fprintf(['File: ', file_find.name,'\n'])

    t = data(:,1); % time [s]
    d = data(:,6) - data(1,6);  % displacement [mm]
    force_corr = mean(data(end-200:end,5)); % end value (end of experiment)
    f_pull = data(:,5) - force_corr; % pull force, corrected for end value [N}   

    %% Analysis

    %% index start/end using displacement data, mean velocity, shear stress

    % velocity from file name
    name_u = find(file_find.name=='u');
    name_u_end = name_u + 1;
    while file_find.name(name_u_end) ~= '_'
        name_u_end = name_u_end + 1;
    end
    vel = str2double(file_find.name(name_u+1:name_u_end-1)); % velocity from file name

    % determine window for smoothening based on velocity form file name
    if vel <= 1
        window = 525; %1500
    elseif vel <= 3
        window = 75;
    elseif vel <= 5
        window = 100;
    elseif vel <= 11
        window = 25;
    elseif vel <= 15
        window = 20;
    elseif vel <= 30
        window = 20;
    else
        window = 15;
    end

    % smoothening data
    d_filt = smoothdata(d,'gaussian',window); % filter data
    
    % compute velocity and acceleration; required for start/end friction response
    v_vec = [0]; % velocity vector
    a_vec = [0]; % acceleration vector
    for i = 2:length(t)-1
        dt = t(i+1) - t(i);
        v = (d_filt(i+1)-d_filt(i-1))/(2*dt);
        a = (d_filt(i+1)-2*d_filt(i)+d_filt(i-1))/(dt^2);
        v_vec(i) = v;
        a_vec(i) = a;
    end
    v_vec(length(t)) = v;
    a_vec(length(t)) = a;

    % determine global start and end based on displacement
    i_glob = 1; % global end range of start point
    while d(i_glob) < 1 % deformation below 1 mm
        i_glob = i_glob + 1;
    end

    i_glob_start = i_glob; % find around i_glob also global start range for start point
    while t(i_glob_start) > t(i_glob)-(2/vel*60)  % move 2 mm in time earlier based on velocity in experiment
        i_glob_start = i_glob_start - 1;
    end

    i_glob_end = i_glob;
    while t(i_glob_end) < t(i_glob) + d(end)/vel*60 % find index for time after full displacement
        i_glob_end = i_glob_end + 1;
    end
    
    % compute start/end friction response based on max/min acceleration around global indices for start/end
    idx_start = find(a_vec(i_glob_start:i_glob)==max(a_vec(i_glob_start:i_glob)))+i_glob_start-1;
    idx_end = find(a_vec(i_glob:i_glob_end)==min(a_vec(i_glob:i_glob_end)))+i_glob-1;

    % shear stress
    tau = f_pull/2500/2*1e3; % [N/mm] -> MPa, times 1e3 results in [kPa] (division by 2, as 2 slip interfaces are present)      

    %% modulus fit
    
    % fit modulus at start-up; correction for t and d data
    [t,d,mod_i] = mod_fit_v1(idx_start,idx_end,t,d,tau,mod_coeff,window);

    %% Show data

    figure(1), hold on
    % color
    if t_wait(meas_i) == 10
        i_color = 1;
    elseif t_wait(meas_i) == 5
        i_color = 5;
    else 
        i_color = 7;
    end

    if range_meas(meas_i) >= 25
        plot(d(idx_start:idx_end),tau(idx_start:idx_end),'-k','color',[0.1 0.1 0.1]*i_color,'linewidth',2.0,'displayname',['Tool-ply, ',num2str(t_wait(meas_i)),' min']);
    else
        if t_wait(meas_i) ~= t_wait_prev
            plot(d(idx_start:idx_end),tau(idx_start:idx_end),'--k','color',[0.1 0.1 0.1]*i_color,'linewidth',2.0,'displayname',['Ply-ply, ',num2str(t_wait(meas_i)),' min']);
        else
            plt = plot(d(idx_start:idx_end),tau(idx_start:idx_end),'--k','color',[0.1 0.1 0.1]*i_color,'linewidth',2.0,'displayname',['Ply-ply, ',num2str(t_wait(meas_i)),' min']);
            set( get( get( plt, 'Annotation'), 'LegendInformation' ), 'IconDisplayStyle', 'off' );
        end
    end

    t_wait_prev = t_wait(meas_i);
end

% plot controls
legend('location','NorthEast','interpreter','latex','box','off')
xlabel('Displacement, d [mm]','interpreter','Latex','fontsize',18)
ylabel('Shear stress, $\tau$ [kPa]','interpreter','Latex','fontsize',18)
set(gca,'fontsize',18)
box on
grid on
axis([0 10 0 200])


