function [tm,ip]=meshmerg(t,tt)
%                  Merging two compatible uniform meshes
% =======================================================================
% function [tm,ip]=meshmerg(t,tt)
% INPUT:
% t  ... vector of the first mesh of length n which MUST be uniform with
%        spacing dt.
% tt ... vector of the second mesh of length N which should be compatible
%        with mesh t in the following sense:
%        (1) tt is uniform having the same spacing dt,
%        (2) both t(1) and t(n) are knots of the infinite continuation of tt
%            (t is a submesh of infinite continuation of tt).
% OUTPUT:
% tm ... = t  if tt==[]
%        = tt if t ==[], otherwise:
%        = mesh on interval [min(t(1),tt(1),max(t(n),tt(N)] with spacing dt
%          into which t an tt will merge if they are compatible (column)
%          (minimal mesh comprising both t and tt).
%        = [] ... if t and tt are not compatible (cannot be merged),
%                 error is reported only if t is not a uniform mesh
% ip ... = [] if t=[] or tt=[], otherwise:
%        = [i1,i2,ii1,ii2] ... pointers to t and tt in tm if tm is not empty:
%                             all(tm(i1:i2)==t)=1, all(tm(ii1:ii2)==tt)=1.

% (C) 1998, V. Vesely, Masaryk University of Brno (Czech Republic)
%
error(nargchk(2,2,nargin));
tol = 5*eps; % tolerance relaxing strict equality of two numbers
tm = []; ip = [];
if isempty(t), tm = tt; return; end
if isempty(tt),tm = t ; return; end
n = length(t); N = length(tt);
if size(t,1)>1 & size(t,2)>1
   error('First mesh t is not a vector.')
end
if size(tt,1)>1 & size(tt,2)>1
   error('Second mesh tt is not a vector.');
end
t = reshape(t,n,1); tt = reshape(tt,N,1);
if n>1, dt = (t(n)-t(1))./(n-1); else dt = 0; end
if any(abs(t.'-linspace(t(1),t(n),n))>tol)
   error('First mesh t is not uniform');
end
if N>1, dtt = (tt(N)-tt(1))./(N-1); else dtt = 0; end
if dt<0 | dtt<0 | (dt~=0 & dtt~=0 & abs(dt-dtt)>tol)
                           % t and tt are incompatible due to bad spacing
   return;
end
if any(abs(tt.'-linspace(tt(1),tt(N),N)>tol))  % tt is not uniform
   return;
end
tmin= min(t(1),tt(1)); tmax=max(t(n),tt(N));
dtm = dt +  dtt*(dt==0);
if dtm>0, Nm = (tmax-tmin)./dtm+1; else Nm = 1+((tmax-tmin)>0); end
tm  = linspace(tmin,tmax,Nm).';
i1  = find(abs(tm- t(1))<=tol); i2  = find(abs(tm- t(n))<=tol);
ii1 = find(abs(tm-tt(1))<=tol); ii2 = find(abs(tm-tt(N))<=tol);
ip = [i1,i2,ii1,ii2];
if length(ip)<4, tm=[]; return; end  % t and tt are not compatible
if (i2-i1+1)~=n | (ii2-ii1+1)~=N, tm=[]; return; end
                                     % another test to be safe
