function merge(ts1,ts2,mergemethod,varargin)
%MERGE Merge two time series objects onto a common time vector
%
%   MERGE(TS1, TS2, 'MERGEMETHOD') creates a new time series object by 
%   merging the time series objects TS1 and TS2 onto a common time vector.
%   The string 'MERGEMETHOD' defines how the time series objects are
%   merged. MERGE replaces both the original time series objects TS1 and  
%   TS2 with the new merged time series object. The options for 
%   'MERGEMETHOD' are
%
%   'Union' - Takes the union of the time vectors over intersecting
%   intervals. The resulting time vectors are equal to the union of the
%   time vectors of TS1 and TS2 on the time interval where the two
%   time vectors overlap
%
%   'Intersection' - Take the intersection of the time
%   vectors. The resulting time vectors  are equal to the intersection of
%   the time vectors of TS1 and TS2.
%
%   'Uniform' - this option requires an additional argument. See below for
%   syntax and description.
%
%   MERGE(TS1, TS2, 'UNIFORM', INCREMENT) creates a new time series object
%   by resampling both time series objects so that they are uniformly 
%   sampled with the interval specified by INCREMENT
%   on the time interval where the time vectors of TS1 and TS2 overlap. 
%   The units of the increment are assumed to be the smaller time units of 
%   the two input timeseries. 
%
%   See also TSDATA.TIMESERIES/TIMESERIES

%   Author(s): James G. Owen
%   Copyright 1986-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.1 $  $Date: 2004/12/26 21:35:08 $

%TO DO: Merge quality 

nargchk(3,inf,nargin);
% Extract time
t1 = ts1.time;
t2 = ts2.time;
if isempty(t1) || isempty(t2)
    error('timeseries:merge:noemptytime',...
        'One or more of the time vectors is empty')
end

% Merge time vectors onto a common basis
[ts1timevec, ts2timevec, outprops,outtrans] = ...
    timemerge(ts1.timeInfo, ts2.timeInfo,t1,t2);

% Find overlapping time interval
interval = [max(ts1timevec(1),ts2timevec(1)) min(ts1timevec(end),ts2timevec(end))];
if interval(2)<=interval(1)
    error('timeseries:merge:strictpos', ...
        'Overlapping time interval must have strictly positive length')
end

% Find output time interval
ts1timevec = ts1timevec(ts1timevec>=interval(1) & ts1timevec<=interval(end));
ts2timevec = ts2timevec(ts2timevec>=interval(1) & ts2timevec<=interval(end));

% Merge the time vectors
switch lower(mergemethod)
case 'union'
    tout = union(ts1timevec,ts2timevec);
    % The output of the 'union.m' function is a vector of time, in which
    % two successive points with even the slightest difference, e.g. 1e-15,
    % will be treated as two different points. In real cases we'd like to
    % treat them as a single point. A threshold has been used to check if
    % any two succesive points are essentially the same after the merge.
    threshold = 1e-10;
    db = diff(tout);
    db = db/mean(db);
    tout = tout([db>threshold;true]); 

    % Parse extra args
    [interpobj1,interpobj2,modcode] = localParseExtras(nargin-3,varargin,ts1,ts2);
    
case 'intersection'
    tout = intersect(ts1timevec,ts2timevec);
    if isempty(tout)
        error('timeseries:merge:emptyset',...
            'There are no intersecting time instants')
    end
    
    % Parse extra args
    [interpobj1,interpobj2,modcode] = localParseExtras(nargin-3,varargin,ts1,ts2);
    
case 'uniform'
   if nargin>=4
      tout = interval(1):varargin{1}:interval(end);
   else
      error('timeseries:merge:syntax', ...
          'A sampling interval must be specified for merging on uniform time vectors')
   end
   
    % Parse extra args
    [interpobj1,interpobj2,modcode] = localParseExtras(nargin-4,varargin,ts1,ts2);
    
otherwise
      error('timeseries:merge:interpsyntax', ...
          'Invalid merge method specification')
end

%% Map the output time interval back into the units and offsets of the input
%% time series
if outtrans.deltaTS==1 % deltaTS==1 means the first time series has been shifted 
    ts1.resample((tout-outtrans.delta)/outtrans.scale{1},interpobj1,modcode);
    ts2.resample(tout/outtrans.scale{2},interpobj2,modcode);
elseif outtrans.deltaTS==2 % deltaTS==1 means the second time series has been shifted 
    ts1.resample(tout/outtrans.scale{1},interpobj1,modcode);
    ts2.resample((tout-outtrans.delta)/outtrans.scale{2},interpobj2,modcode);
else
    ts1.resample(tout/outtrans.scale{1},interpobj1,modcode);
    ts2.resample(tout/outtrans.scale{2},interpobj2,modcode);
end


function [interpobj1,interpobj2,modcode] = localParseExtras(nargs,argin,ts1,ts2)

%% Has a custom interpolation method been defined?
interpobj1 = ts1.DataInfo.Interpolation;
interpobj2 = ts2.DataInfo.Interpolation;
if nargs>=1 && ~isempty(argin{1})
    if ischar(argin{1}) % Interpolation method specified by a string
        interpobj1 = tsdata.interpolation(argin{4});
        interpobj2 = tsdata.interpolation(argin{4});
    elseif isa(argin{1},'tsdata.interpolation')
        interpobj1 = argin{1};
        interpobj2 = argin{1};
    else
        error('timeseries:merge:invinterp',...
            'Invalid specification of interpolation method')
    end
end

%% Has a modifed quality code been defined?
if nargs>=2 && ~isempty(argin{2})
    if isempty(ts1.Quality) || isempty(ts2.Quality)
        error('timeseries:merge:noqual','Quality vector must be initialized')
    end
    if isnumeric(argin{2}) && isscalar(argin{2}) && argin{2}>=0
        if floor(argin{2})-argin{2}<0
            warning('Data status code has been converted to an integer')
        end
        modcode = floor(argin{2});
    else
        error('timeseries:merge:invcode','Invalid data status code')
    end
else
    modcode = [];
end