function kDI = kDIfromtensor(tensorelements)
kDI = zeros(length(tensorelements(:,1)),9);

for jj = 1:length(tensorelements(:,1))
    
    % build tensor
    kxx = tensorelements(jj,1);
    kyy = tensorelements(jj,2);
    kzz = tensorelements(jj,3);
    kxy = tensorelements(jj,4);
    kyz = tensorelements(jj,5);
    kzx = tensorelements(jj,6);
    
    if kxx >= kyy -1e-16 && kxx <= kyy + 1e-16 && kyy >= kzz-1e-16 && kyy <= kzz+ 1e-16 && kxx >= kzz -1e-16 && kxx <= kzz + 1e-16
        kmean = (kxx(jj) + kyy(jj) + kzz(jj))/3;
        kDI(jj,:) = [kmean 0 0 kmean 90 0 kmean 0 90];
    else
    
    tensor = [kxx kxy kzx; kxy kyy kyz; kzx kyz kzz];
    
    % compute eigenvalues and eigenvectors
    [v,u] = eig(tensor);
    
   u(1,2) = NaN;
   u(2,3) = NaN;
   u(1,3) = NaN;
   u(2,1) = NaN;
   u(3,2) = NaN;
   u(3,1) = NaN;
    
    % sort into (k1, v1) (k2, v2) and (k3, v3)
    a1 = find(u == max(max(u)));
    if length(a1) ==2
        a1 = a1(1);
    end
    if a1 == 1
        k1 = u(1,1);
        u(1,1) = NaN;
        v1 = v(:,1);
    elseif a1 == 5
        k1 = u(2,2);
        u(2,2) = NaN;
        v1 = v(:,2);
    elseif a1 == 9
        k1 = u(3,3);
        u(3,3) = NaN;
        v1 = v(:,3);
    end
    
    a2 = find(u == max(max(u)));
    if length(a2) ==2
        a2 = a2(1);
    end
    if a2 == 1
        k2 = u(1,1);
        u(1,1) = NaN;
        v2 = v(:,1);
    elseif a2 == 5
        k2 = u(2,2);
        u(2,2) = NaN;
        v2 = v(:,2);
    elseif a2 == 9
        k2 = u(3,3);
        u(3,3) = NaN;
        v2 = v(:,3);
    end
    a3 = find(u == max(max(u)));
    if a3 == 1
        k3 = u(1,1);
        u(1,1) = NaN;
        v3 = v(:,1);
    elseif a3 == 5
        k3 = u(2,2);
        u(2,2) = NaN;
        v3 = v(:,2);
    elseif a3 == 9
        k3 = u(3,3);
        u(3,3) = NaN;
        v3 = v(:,3);
    end
    
    % compute D and I from each eigenvector
    % v1
    
    vlength = vectorlength(v1);
    
    x = v1(1)/vlength;
    y = v1(2)/vlength;
    z = v1(3)/vlength;
    
    if x>=0 && y>=0
        D1 = atand(y/x);
    elseif x<0 && y>=0
        D1 = 180 + atand(y/x);
    elseif x<0 && y<0
        D1 = 180 + atand(y/x);
    elseif x>=0 && y<0
        D1 =  360 + atand(y/x);
    end
    
    if isnan(D1)
        D1 = 0;
    end
    
    I1 = asind(z);
    
    if I1 < 0
        I1 = - I1;
        D1 = D1 + 180;
    end
    if D1 > 360
        D1 = D1 -360;
    end
    
    % v2
    
    vlength = vectorlength(v2);
    
    x = v2(1)/vlength;
    y = v2(2)/vlength;
    z = v2(3)/vlength;
    
    if x>=0 && y>=0
        D2 = atand(y/x);
    elseif x<0 && y>=0
        D2 = 180 + atand(y/x);
    elseif x<0 && y<0
        D2 = 180 + atand(y/x);
    elseif x>=0 && y<0
        D2 =  360 + atand(y/x);
    end
    
    if isnan(D2)
        D2 = 0;
    end
    
    I2 = asind(z);
    
    if I2 < 0
        I2 = - I2;
        D2 = D2 + 180;
    end
    if D2 > 360
        D2 = D2 -360;
    end
    
    
    % v3
    
    vlength = vectorlength(v3);
    
    x = v3(1)/vlength;
    y = v3(2)/vlength;
    z = v3(3)/vlength;
    
    if x>=0 && y>=0
        D3 = atand(y/x);
    elseif x<0 && y>=0
        D3 = 180 + atand(y/x);
    elseif x<0 && y<0
        D3 = 180 + atand(y/x);
    elseif x>=0 && y<0
        D3 =  360 + atand(y/x);
    end
    
    if isnan(D3)
        D3 = 0;
    end
    
    I3 = asind(z);
    
    if I3 < 0
        I3 = - I3;
        D3 = D3 + 180;
    end
    if D3 > 360
        D3 = D3 -360;
    end
    
    kDI(jj,:) = [k1 D1 I1 k2 D2 I2 k3 D3 I3];
    
    clear v u v1 v2 v3 k1 D1 I1 k2 D2 I2 k3 D3 I3 vlength x y z
    end
    
end
