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