function [D,tdt,tphi]=split(uE,uN,smp)

% A function for shear-wave splitting
% input = east-component, north-component and sampling period (s)
% output = 2D map of determinant of split-rotated correlation matrix
% J.A. Olive, July 2010, after Silver and Chan, 1991


utotE=uE;
utotN=uN;
N=length(utotE);
tmax=(N-1)*smp;

% phi
Nphi=50;
phimin=-90;
phimax=90;
% dt
Ndt=50;
dtmin=0;
dtmax=3;
%
tphi=linspace(phimin,phimax,Nphi);
tdt=linspace(dtmin,dtmax,Ndt);


% cross-correlating the components
lag=linspace(-tmax,tmax,2*N-1);
CEN=xcorr(utotE,utotN);
CNE=xcorr(utotN,utotE);
CEE=xcorr(utotE,utotE);
CNN=xcorr(utotN,utotN);

D=zeros(Ndt,Nphi);
R=zeros(2,2);
C0=zeros(2,2);
Cdt=C0;

for i=1:Ndt
    for j=1:Nphi
        
        temp=90-tphi(j);
        R(1,1)=cosd(temp);
        R(2,2)=cosd(temp);
        R(1,2)=-sind(temp);
        R(2,1)=sind(temp);
        
        C0(1,1)=interp1(lag,CEE,0);
        C0(2,2)=interp1(lag,CNN,0);
        C0(1,2)=interp1(lag,CEN,0);
        C0(2,1)=interp1(lag,CNE,0);
        
        Cdt(1,1)=interp1(lag,CEE,tdt(i));
        Cdt(2,2)=interp1(lag,CNN,tdt(i));
        Cdt(1,2)=interp1(lag,CEN,tdt(i));
        Cdt(2,1)=interp1(lag,CNE,tdt(i));
        
        a=0;
        b=0;
        c=0;
        for u=1:2
            for v=1:2
                a=a+R(1,u)*C0(u,v)*R(1,v);
                b=b+R(2,u)*C0(u,v)*R(2,v);
                c=c+R(1,u)*Cdt(u,v)*R(2,v); 
            end
        end
        
        Ct(1,1)=a;
        Ct(2,2)=b;
        Ct(2,1)=c;
        Ct(1,2)=c;  
        [~,valp]=eig(Ct);
        l2=valp(1,1); % smallest eigenvalue
        D(i,j)=l2;
 
    end
end

tphi=90-tphi;


figure

contourf(tdt,tphi,D'./min(min(D')))
colorbar
hold on
ylabel('\phi','fontsize',20)
xlabel('\delta t','fontsize',20)
title('Map of \lambda_2 / min(\lambda_2)','fontsize',20)


