pwlint(x,y)

function [Apos, Aneg, Anet] = pwlint(x,y)
% PWLINT.m
%  Calculate the positive and negative area under a curve
%  using trapezoids, assuming a piecewise linear function
%
% Requires insertrows:
%   http://www.mathworks.com/matlabcentral/fileexchange/9984
%
% Inspired by this thread:
%   http://www.mathworks.com/matlabcentral/newsreader/view_thread/252131
%
% Usage:
%   [Apos, Aneg, Anet] = pwlint(x,y)
%        x: vector of independent variable
%        y: vector of dependent variable
%     Apos: area between y=0 and positive y values
%     Aneg: area between y=0 and negative y values
%     Anet: Apos+Aneg; equal to trapz(x,y)
%
% See also TRAPZ, SUM

% v0.9 (Jan 2012) by Andrew Davis -- [e-mail]

% Check inputs. Always, always, always check inputs! (Sometimes)
if ~isvector(x) || ~isvector(y),
   whos x y;
   error('x and y must be 1D vectors');    
end
if size(x) ~= size(y),
   whos x y;
   error('x and y must be the same size');
end
if isrow(x), x=x'; end;             % column vectors
if isrow(y), y=y'; end;

% Find zero-crossings and insert roots into x
i = 1:(length(x)-1);                % leave off index of last point
dx = diff(x);
icross = find(((y(i) > 0) & (y(i+1) < 0)) | ((y(i) < 0) & (y(i+1) > 0)));
x0 = y(icross).*dx(icross)./(y(icross)-y(icross+1)) + x(icross);

xplus = insertrows(x, x0, icross);
yplus = insertrows(y, 0, icross);
iplus = 1:(length(xplus)-1);
dxplus = diff(xplus);

% Find and sum positive, negative trapezoids
ipos = find((yplus(iplus) >= 0) & (yplus(iplus+1) >= 0));
Apos = sum(0.5.*dxplus(ipos).*(yplus(ipos) + yplus(ipos+1)));

ineg = find((yplus(iplus) <= 0) & (yplus(iplus+1) <= 0));
Aneg = sum(0.5.*dxplus(ineg).*(yplus(ineg) + yplus(ineg+1)));

Anet = Apos + Aneg;