function [D1, D2, D3, D4] = sorting(D1, D2, D3, D4)
%**************************************************************************
% CLUSTER project: Magnetometerdata sorter         Autor: Christian Maszl * 
% Unit System = [SI]                  Mail: e0025754@student.tuwien.ac.at *
%**************************************************************************
%
% This script sort the magnetometer data of the four satellites according
% to the time. "Accuracy" set the maximum time difference between the 
% measurements of one "event".
%
accuracy = 8e-6; searched1 = 1; searched2 = 1; searched3 = 1; 
k = 1; l = 0; setmatch1=0; setmatch2=0;setmatch3=0; startnext=1;
a = 0; b = 0; c = 0; d = 0; n=0; j1=1; j2=1; j3=1;
%Find the minimum length of the input data
mindim=min([length(D1) length(D2) length(D3) length(D4)]);
% Preallocate output arrays
    X1 =zeros(mindim,7);
    X2 =zeros(mindim,7);
    X3 =zeros(mindim,7);
    X4 =zeros(mindim,7);
    X  =zeros(mindim,4);
% [Ok] Find the data list with the earliest start value
% Permutation of indices 1 to 4
P = perms([1 2 3 4]);
for i=1:24
    if (eval(['D',int2str(P(i,1)),'(1,1)']) <= ...
        eval(['D',int2str(P(i,2)),'(1,1)']) &  ...
        eval(['D',int2str(P(i,1)),'(1,1)']) <= ... 
        eval(['D',int2str(P(i,3)),'(1,1)']) &  ...
        eval(['D',int2str(P(i,1)),'(1,1)']) <= ...
        eval(['D',int2str(P(i,4)),'(1,1)']) )
        % Name of the data file with the earliest start value = a
        a = P(i,1); b = P(i,2); c = P(i,3); d = P(i,4); break;
    else
    end
end
P2 = perms([b c d]);
for j=1:6
    if (eval(['D',int2str(P2(j,1)),'(1,1)']) <= ...
        eval(['D',int2str(P2(j,2)),'(1,1)']) &  ...
        eval(['D',int2str(P2(j,1)),'(1,1)']) <= ... 
        eval(['D',int2str(P2(j,3)),'(1,1)']) )
        b = P2(j,1); c = P2(j,2); d = P2(j,3); break;
    else    
    end
end
if (eval(['D',int2str(c),'(1,1)']) <= ...
    eval(['D',int2str(d),'(1,1)']) )
else
    x = c; c = d; d = x; 
end
% Construct arrays for search
% Array with the latest start time X4
% Row X3<X1<X2<X4
% Write data into X1-4
    eval(['X4=D',int2str(d),';']);
    eval(['X3=D',int2str(a),';']);
    eval(['X2=D',int2str(c),';']);
    eval(['X1=D',int2str(b),';']);
% Length of the input list's
    enda=eval(['size(D',int2str(a),',1);']);
    endb=eval(['size(D',int2str(b),',1);']);
    endc=eval(['size(D',int2str(c),',1);']);
    endd=eval(['size(D',int2str(d),',1);']);
    clear P P2 x D1 D2 D3 D4;
% Perform search with the latest search time
for i=startnext:endd
     for j1=searched1:endc
         if (abs(X4(i,1)-X2(j1,1))<=accuracy)            
            % First hit in the if statement 
                firsthit=n+1;
            % Set index to which all others will be compared
                if (firsthit==1 & startnext==1)
                    startnext=i;
                else
                end
            %Decrease search length    
                searched1 = j1;
            % Found a match
                setmatch1 = 1;
            % Break match found
                break;
         else              
            % Found no match
                setmatch1 = 0;
         end 
         % Prevent scanning the whole list when no match can be found
         if (X4(i,1)+accuracy<X2(j1,1))
             break;
         else
         end
     end
     for j2=searched2:endb
         if (abs(X4(i,1)-X1(j2,1))<=accuracy)           
            %Decrease search length    
                searched2=j2;
            % Found a match
                setmatch2=1;
            % Break match found
                break;
        else
            % Found no match
                setmatch2=0;  
        end
        % Prevent scanning the whole list when no match can be found
        if (X4(i,1)+accuracy<X1(j2,1))
             break;
        else
        end
    end
    for j3=searched3:enda
         if (abs(X4(i,1)-X3(j3,1))<=accuracy)  
            %Decrease search length    
                searched3=j3;
            % Found a match
                setmatch3=1;
            % Break match found
                break;
        else
            % Found no match
                setmatch3=0;  
        end
        % Prevent scanning the whole list when no match can be found
        if (X4(i,1)+accuracy<X3(j3,1))
            break;
        else
        end
    end
    if (setmatch1 & setmatch2 & setmatch3)
          % Store matching lines
            X(k,:)=[i j1 j2 j3];
          % Increment line pointer
            k = k+1;
    elseif (j1==endc |  j2==endb | j3==enda)
          break;
    else
    end
end
%Delete dispensable variables
    clear j1 j2 j3 i setmatch1 setmatch2 setmatch3;
    clear searched1 searched2 searched3; 
    clear con1 con2 con3 con4;
    clear D1 D2 D3 D4 j k l mindim;
    clear k1 k2 k3 firsthit startnext;
    clear enda endb endc endd;
% Find "nonzero" length of the X array
    n=size(find(X(:,1)),1);
% Store these values in an data array
    DATA(:,:)=X(1:n,:); 
% Use data array to extract correct measurements
% Row X3<X1<X2<X4
% Write data into D1-4
    eval(['D',int2str(d),'=X4(DATA(:,1),:);']);
    eval(['D',int2str(a),'=X3(DATA(:,4),:);']);
    eval(['D',int2str(c),'=X2(DATA(:,2),:);']);
    eval(['D',int2str(b),'=X1(DATA(:,3),:);']);
% Clear all except the sortet data for the main programm
    clear DATA a b c d con1 con2 con3 con4 ;
    clear X X1 X2 X3 X4 n accuracy;
% Keep accuracy for a possible GUI, if you like.
