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

%% INFO

% Read and process friction data

% Figure concering effect of normal pressure and release agent for ESAFORM 2022 conference paper

% Matlab version: R2021a

%% initialisation
max_slip_length = 10; % mm, maximum slip length for evaluation of long-time shear stress
mod_coeff = 0.4; %threshold for modulus fit, end; max difference local and average modulus
meas_range = [33 34 35 20 21 22 41 42 43 74 75 76 62 63 64 85 86 87]; % unique labels of specimens
pressure = [5 5 5 15 15 15 45 45 45 5 5 5 15 15 15 45 45 45]; % corresponding pressures
vel = 25; % corresponding rate

figure(1), clf(1)

%% Go through the measurements

V_mean_vec = [];  % store mean rate in loop
tau_mean_vec = []; % store mean long-time shear stress in loop
curves_store = []; % store: disp, tau
    
for meas_i = 1:length(meas_range)
    
    %% Load data
    file_find = dir(['TP','*',num2str(meas_range(meas_i)),'*']);
    file_name = file_find.name; 
    fid = fopen(file_name);
    C = textscan(fid,'%f');
    data = zeros(int32(length(C{1})/13),13);
    for i = 0:length(C{1})-1
        i_row = 1+floor(i/13);
        i_col = 1+mod(i,13);
        data(int32(i_row),int32(i_col)) = C{1}(i+1);
    end
    fprintf(['File: ', file_find.name,'\n'])
    fclose(fid);

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

    %% Analysis

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

    % smoothening data
    window = 100;
    d_filt = smoothdata(d,'gaussian',window); % filter data
    
    % compute velocity and acceleration 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
    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; % find around i_glob also global start range for start point
    while d(i_glob_end) < min(max(d),12)  % move 12 mm in time based on velocity in experiment
        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;

    % mean velocity
    V_mean = (d(idx_end)-d(idx_start))/(t(idx_end)-t(idx_start)); % mm/s
    V_mean_vec = [V_mean_vec, V_mean]; % store values

    % 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 for correction on t and d data

    % filter displacement data (multiple inputs: averaging)
    d_red = []; % reduced displacement vector
    tau_red = []; % reduced shear stress vector
    idx_avg_end = idx_start; % find index for range for modulus fit
    while d(idx_avg_end) < min(3,max(d))  % stop after 3 mm
        idx_avg_end = idx_avg_end + 1;
    end
    d_prev = d(idx_start);
    i_count = 1;
    i_red = 1;
    tau_tot = tau(idx_start);
    for i = 1:(idx_avg_end-idx_start)  % data for averaging (because same displacmeent input)
        if d(i+idx_start) <= d_prev
            tau_tot = tau_tot + tau(i+idx_start);
            i_count = i_count + 1;
        else
            d_red(i_red) = d_prev;   % include data point in reduced vector
            tau_red(i_red) = tau_tot/i_count;  % use average shear stress value         
            d_prev = d(i+idx_start);
            i_count = 1;  % start again with counting 
            tau_tot = tau(i+idx_start); % tau_tot equals subsequent data point
            i_red = i_red + 1; % index for reduced vector
        end
    end
    
    % derivative tau wrt displacement
    dtau_dd_first = [0 0]; % first derivative of tua to d (d tau / d d)
    for i_der = 3:length(tau_red)-2
        dtau_dd_first(i_der) = (-0.5*tau_red(i_der-1)+0.5*tau_red(i_der+1))/(0.5*(d_red(i_der)-d_red(i_der-1))+0.5*(d_red(i_der+1)-d_red(i_der)));
    end

    idx_mod_start = 1;
    while dtau_dd_first(idx_mod_start) < mod_coeff*max(dtau_dd_first) % use fraction of max derivative for start range fit modulus
        idx_mod_start = idx_mod_start + 1;
    end
    idx_mod_start = idx_mod_start - 1;
    if idx_mod_start < 1
        idx_mod_start = 1;
    end
    idx_mod_end = find(dtau_dd_first==max(dtau_dd_first)); % start after the peak
    while dtau_dd_first(idx_mod_end) > mod_coeff*max(dtau_dd_first) & idx_mod_end < length(dtau_dd_first) % use fraction of max derivative for start range fit modulus
        idx_mod_end = idx_mod_end + 1;
    end

    % modulus fit in selected range
    tau_mod_fit = tau_red(idx_mod_start:idx_mod_end);
    d_mod_fit = d_red(idx_mod_start:idx_mod_end);
    p = polyfit(d_mod_fit,tau_mod_fit,1);
    tau_mod = [p(1)*(min(d_mod_fit)-0.1)+p(2) p(1)*(max(d_mod_fit)+0.1)+p(2)];  % vector to plot straight line fit
    d_mod = [min(d_mod_fit)-0.1 max(d_mod_fit)+0.1]; % vector to plot straight line fit

    % correction for displacement and time
    d_corr = -p(2)/p(1);
    idx_corr = idx_start;
    while d(idx_corr) < d_corr
        idx_corr = idx_corr + 1;
    end
    d = d - d_corr; 
    t_corr = t(idx_corr);
    t = t - t_corr; 
    
    % store average curves
    curves_store{length(V_mean_vec)} = [d(idx_start:idx_end) tau(idx_start:idx_end)];
    
    % plot data
    figure(1), hold on
    plot(d(idx_start:idx_end),tau(idx_start:idx_end),'displayname',['Meas: ', num2str(meas_range(meas_i))])
    plt = plot(d_mod-d_corr,tau_mod,'-k');
    set( get( get( plt, 'Annotation'), 'LegendInformation' ), 'IconDisplayStyle', 'off' );
    plt = plot(d_mod_fit-d_corr,tau_mod_fit,'ko');
    set( get( get( plt, 'Annotation'), 'LegendInformation' ), 'IconDisplayStyle', 'off' );
    
end

%% Compute average curves

figure(1), clf(1), hold on
color_i = 1;
color_i_RA = 1;

% determine ranges
data_ranges = [1 4 7 10 13 16 19]; %start meas, up to meas (not included), up to meas (not included) etc

for i = 1:length(data_ranges)-1

    % data range for mean curve
    j_start = data_ranges(i);
    j_end = data_ranges(i+1)-1;
    j_max = 1;
    length_max = 1;
    d_std_thres = 0.75; % for showing std bars
    d_std = 2; % distance between std bars

    % find index with longest dataset
    for j = j_start:j_end
        if length(curves_store{j}) > length_max
            j_max = j;
            length_max = length(curves_store{j});
        end
    end
    
    if pressure(j_start) ~= pressure(j_end)
        fprintf('Unequal variables \n')
    end       

    % find start and end indices
    d_base = curves_store{j_max}(:,1);
    k_start = min(find(d_base>0));
    k_end = min(find(d_base>10));
    if isempty(k_end)
         k_end = length(d_base);
    end

    % determine mean curves
    d_curve = zeros(1,k_end-k_start+1);
    tau_curve = zeros(1,k_end-k_start+1);
    std_curve = zeros(1,k_end-k_start+1);
    k = k_start;
    while k <= k_end
        d_k = d_base(k);
        d_curve(k-k_start+1) = d_k;

        % find clostes tau values; interpolation
        tau_k = zeros(1,j_end-j_start+1);
        idx_zero = []; % account for end of curves; no data
        for j = j_start:j_end
            if j == j_max
                tau_k(j-j_start+1) = curves_store{j}(k,2);
            else
                idx_left = max(find(curves_store{j}(:,1)<d_k));
                idx_right = min(find(curves_store{j}(:,1)>d_k));
                if ~isempty(idx_right) && ~isempty(idx_left)

                    % interpolation to find tau_k
                    tau_L = curves_store{j}(idx_left,2); 
                    tau_R = curves_store{j}(idx_right,2);
                    d_L = curves_store{j}(idx_left,1);
                    d_R = curves_store{j}(idx_right,1);
                    tau_k(j-j_start+1) = tau_L + (tau_R-tau_L)*((d_k-d_L)/(d_R-d_L));
                else
                    idx_zero = [idx_zero, j-j_start+1];
                end
            end
        end

        tau_k(idx_zero) = [];
        tau_curve(k-k_start+1) = mean(tau_k);
        std_curve(k-k_start+1) = std(tau_k);
        if d_k > d_std_thres
            d_std_thres = d_std_thres + d_std;
            plt = plot([d_k d_k],[mean(tau_k)-std(tau_k) mean(tau_k)+std(tau_k)],'-k_');
            set( get( get( plt, 'Annotation'), 'LegendInformation' ), 'IconDisplayStyle', 'off' );
        end 
        k = k + 1;
    end
    
    if i <= 3
        plot(d_curve,tau_curve,'-k','color',[0.1 0.1 0.1]*color_i,'linewidth',2.0,'displayname',[num2str(pressure(j_max)),' kPa']);
        color_i = color_i + 3;
    elseif i == 4
        plot(d_curve,tau_curve,'--k','color',[0.1 0.1 0.1]*color_i_RA,'linewidth',2.0,'displayname','Release Agent');
        color_i_RA = color_i_RA + 3;
    else
        plt = plot(d_curve,tau_curve,'--k','color',[0.1 0.1 0.1]*color_i_RA,'linewidth',2.0,'displayname','Release Agent');
        set( get( get( plt, 'Annotation'), 'LegendInformation' ), 'IconDisplayStyle', 'off' );
        color_i_RA = color_i_RA + 3;
    end
end

% plot controls
figure(1), hold on
axis([0 10 0 150])
box on
grid on
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)