
% This function performs union with a range tolerance, the merged array
% using input A as coordinates. See Matlab function 'union' for a detail
% introduction on exact array union.

function [c,idA,idB, ampsA_all, ampsB_all]...
    = union_peaks(A, B, ampsA, ampsB, dt)
pA=1;
pB=1;
lenA=length(A);
lenB=length(B);

[sA,idxA]=sort(A,'ascend');
[sB,idxB]=sort(B,'ascend');

flag_last=0; % =0 if last point is from A
             % =1 if last point is from B
if sA(1)<sB(1)
    tc=sA(1);
    idA(idxA(1))=1;
    pA=pA+1;
else
    tc=sB(1);
    idB(idxB(1))=1;
    pB=pB+1;
    flag_last=1;
end
count=1;

% [nA,mA]=size(A);
% if nA > mA
%     flag_column=1;
% else
%     flag_column=0;
% end

ub=max(sA(end),sB(end))+1;
while pA<=lenA | pB<=lenB    
   
    if sA(pA)<=sB(pB) 
%        pA=pA+1;
        if sA(pA)-tc(count)>=dt
            count=count+1;
            tc(count)=sA(pA);
        elseif flag_last==1 % here make sure using input A as coordinates
             tc(count)=sA(pA);  
        end
        flag_last=0;
        idA(idxA(pA))=count;
        pA=pA+1;
        
        if pA>lenA
            sA(end+1)=ub;
        end
    elseif sA(pA)>sB(pB) 
        if sB(pB)-tc(count)>=dt
            count=count+1;
        	tc(count)=sB(pB);            
            flag_last=1;
        end
        
        idB(idxB(pB))=count;
        pB=pB+1;
        
        if pB>lenB
            sB(end+1)=ub;
%            sB=[sB ub];
        end
    end       
end
c=tc;

lnall=length(c);
ampsA_all=zeros(lnall,1);
ampsA_all(idA)=ampsA;
%ampsNB_all(invIdxNB(find(invIdxNB>0)))=ampsNB;
ampsB_all=zeros(lnall,1);
ampsB_all(idB)=ampsB;
%%%
%%%%%%%% Make sure the largest amptitude in duplicate positions  
%%%
% check if there is duplicated point in idA
for k=1:length(idA)-1
   if  idA(k)~=idA(k+1)
       continue;
   end
   pos=find(idA==idA(k));
   tmaxamp=max(ampsA(pos));

   ampsA_all(idA(pos(1)))=tmaxamp;
end

% check if there is duplicated point in idA
for k=1:length(idB)-1
   if  idB(k)~=idB(k+1)
       continue;
   end
   pos=find(idB==idB(k));
   tmaxamp=max(ampsB(pos));

   ampsB_all(idB(pos(1)))=tmaxamp;
end


