function [t,d,mod_i] = mod_fit(idx_start,idx_end,t,d,tau,mod_coeff,window)
%mod_fit: fit line on first part of friction response, also used to correct
%time and displacement data
%   Version: 11-11-2020

% Cleaned for use in ESAFORM 2022 conference paper figures   


% 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)-d(idx_start) < 3 && idx_avg_end < idx_end   % 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)  
    if d(i+idx_start) <= d_prev
        tau_tot = tau_tot + tau(i+idx_start); % data for averaging (because same displacmeent input)
        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
dtau_dd_first = smoothdata([dtau_dd_first, 0, 0],'gaussian',window);

% determine start and end for fit range
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) % use fraction of max derivative for start range fit modulus
    idx_mod_end = idx_mod_end + 1;
end

% mod 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);
mod_i = p(1);
tau_mod = [p(1)*(min(d_mod_fit)-0.1)+p(2) p(1)*(max(d_mod_fit)+0.1)+p(2)];
d_mod = [min(d_mod_fit)-0.1 max(d_mod_fit)+0.1];

% 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; 

end

