% Function: PlotUnevenDiff(x1,y1,x2,y2,extrap)
% Purpose:  Plot y2-y1 vs. x.  x values are taken to be x1(union)x2.
%           If extrap == 1, perform extrapolation. 
%
% Inputs:
%     x1 : list of x-values (e.g. time)
%     x2 : another list of x-values
%     y1 : dependent variable values corresponding to x1
%     y2 : dependent variable values corresponding to x2
% extrap : flag determining whether extrapolation is performed. 1 ==> true
%   axes : figure axis handle (optional)
%
% Outputs
%   Plots graph of (x1(union)x2,y2-y1) to screen.

function [x,y1_new,y2_new] = PlotUnevenDiff(x1,y1,x2,y2,extrap,axes, linespec)

    if (nargin <= 5)
        figure;
        axes = gca;
    end
    
    if (nargin <= 6)
        linespec = ''
    end
    
    % Collapse arrays to contain unique x values.  If repeated x values
    % occur, all but the first are thrown out (with their corresponding
    % y values)
    [x1_unique, indices1, null] = unique(x1);
    [x2_unique, indices2, null] = unique(x2);
    y1_unique = y1(indices1);
    y2_unique = y2(indices2);
    
    % Get new array of x values: x1(union)x2
    x = unique([x1_unique; x2_unique]); % x1,x2 assumed to be colum vectors
    sort(x);
    
    % Interpolate to get new y values
    if (1 == extrap)
        y1_new = interp1(x1_unique,y1_unique,x,'linear','extrap');
        y2_new = interp1(x2_unique,y2_unique,x,'linear','extrap');
    else
        y1_new = interp1(x1_unique,y1_unique,x,'linear');
        y2_new = interp1(x2_unique,y2_unique,x,'linear');
    end
    
    % Plot diff
    plot(axes,x,y2_new-y1_new, linespec);
end