Bi-linear (piecewise) curve fitting

T1 = readtable('testdata.xlsx');
x = T1.Displacement;
y = T1.BaseForce;

% top flat curve (apex)
[m,idy] = max(y);
ym = m*ones(size(y));

plot(x,y,'*-',x,ym,'g','Linewidth',3);
hold on

% red curve ( "yellow" net area must be zero)
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);

% let's find out which second point of the curve let us find the min area
xx = x(1:idy);
yy = y(1:idy);

for ck = 2:idy
    x2 = x(ck);
    y2 = y(ck);
    slope = (y2-y1)/(x2-x1);
    yl = y1 + slope*(xx - x1);
    % compute yellow area A
    A(ck) = trapz(xx,(yy-yl));
    
end
 

% find zero crossing value for A (with linear interpolation)
threshold = 0; 
xx_zc = find_zc(xx,A,threshold);
yy_zc = interp1(xx,yy,xx_zc);
slope = (yy_zc-y1)/(xx_zc-x1);
yl = y1 + slope*(xx - x1);
plot(xx,yl,'r','Linewidth',3);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end