hudson_map.m

% !-*- coding:UTF-8 -*-
% @Time : 2025/04/29
% @File : hudson_map
% @Software: Matlab 
% Author: Jiefang Song  jfang97@foxmail.com 

% Main function
clear;
clc;
[M11, M22, M33, M23, M13, M12] = textread('E:\....\MTs.dat', '%f%f%f%f%f%f', 'commentstyle', 'matlab');
A = [M11 M22 M33 M23 M13 M12];
drawhudsonnet_density(A);   


% main function
function drawhudsonnet_density(A)

% 计算所有 (u,v)
u = nan(size(A,1),1);
v = nan(size(A,1),1);

for i = 1:size(A,1)
    M_MATRIX = [A(i,1) A(i,2) A(i,3);
                A(i,2) A(i,4) A(i,5);
                A(i,3) A(i,5) A(i,6)];
    try
        [k,t] = ae_srctype(M_MATRIX);
        [u(i),v(i)] = ae_uvt(t,k);
    catch
        u(i) = nan;
        v(i) = nan;
    end
end

% generate grid
xlin = linspace(-4/3, 4/3, 300);
ylin = linspace(-1, 1, 300);
[X,Y] = meshgrid(xlin, ylin);

% density
valid = ~isnan(u) & ~isnan(v);
[f,xi] = ksdensity([u(valid), v(valid)], [X(:), Y(:)]);

F = reshape(f, size(X));  % metraix


% Mask processing
bx = [0 4/3 0 -4/3 0];
by = [1 1/3 -1 -1/3 1];
mask = inpolygon(X, Y, bx, by);
F(~mask) = NaN; 


% 画密度云图
p = pcolor(X, Y, F);
set(p, 'EdgeColor', 'none');  % delete grid
hold on;
p = pcolor(X, Y, F);
set(p, 'EdgeColor', 'none');  % delete grid
colormap(parula);
cb = colorbar;
ylabel(cb, 'AE Density', 'FontSize', 15, 'FontName', 'Times New Roman');
axis equal;
axis tight;
set(gca, 'XLim', [-4/3 4/3], 'YLim', [-1 1]);
hold on;

VV = linspace(-1, 1, 200)';  
levels = [-0.5, 0, 0.5];     

lightGray = [0.7 0.7 0.7];   % LightGrey

for i = 1:length(levels)
    k = levels(i);
    [x1, y1] = ae_uvt(VV, k * ones(size(VV)));  
    [x2, y2] = ae_uvt(k * ones(size(VV)), VV); 
    
    darkerGray = [0.5 0.5 0.5];  % Grey
    line(x1, y1, 'Color', darkerGray, 'LineStyle', '--', 'LineWidth', 1);
    line(x2, y2, 'Color', darkerGray, 'LineStyle', '--', 'LineWidth', 1);
end

drawhudsonframe();

% Add caption（DC, CLVD, Explosion ...）
offset = 0.02;
ha = 'HorizontalAlignment';
va = 'VerticalAlignment';
fontopt = {'FontSize', 13, 'FontName', 'Times New Roman'};

[x,y] = ae_uvt(1,1);     text(x,y+offset,     'Explosion',      ha,'Center', va,'Bottom', fontopt{:});
[x,y] = ae_uvt(1,-1);    text(x,y-offset,     'Implosion',      ha,'Center', va,'Top',    fontopt{:});
[x,y] = ae_uvt(1,0);     text(x,y-offset,     'CLVD (-)',       ha,'Center', va,'Top',    fontopt{:});
[x,y] = ae_uvt(1,-5/9);  text(x,y-offset,     'Anticrack',      ha,'Left',   va,'Top',    fontopt{:});
[x,y] = ae_uvt(-1,0);    text(x,y-offset,     'CLVD (+)',       ha,'Center', va,'Top',    fontopt{:});
[x,y] = ae_uvt(-1,5/9);  text(x,y+offset,     'Tensile Crack',  ha,'Right',  va,'Bottom', fontopt{:});
[x,y] = ae_uvt(0,0);     text(x,y+offset,     'DC',             ha,'Center', va,'Bottom', fontopt{:});
[x,y] = ae_uvt(-1,1/3);  text(x,y,            'LVD (+)',        ha,'Right',  va,'Bottom', fontopt{:});
[x,y] = ae_uvt(1,-1/3);  text(x,y,            'LVD (-)',        ha,'Left',   va,'Top',    fontopt{:});


hold off;
title('Hudson Diagram Density Plot');
xlabel('u');
ylabel('v');
grid off;

set(gca, 'FontSize', 16, 'FontName', 'Times New Roman'); 
set(gca,'Visible','off'); 

set(gca,'Visible','off');

end

%% Hudson frame
function drawhudsonframe()
    line([0 4/3 0 -4/3 0]', [1 1/3 -1 -1/3 1]', 'Color','k', 'LineWidth',1.5);
    line([-1 1]', [0 0]', 'Color','k');
    line([0 0]', [-1 1]', 'Color','k');
    
end

function [k,T] = ae_srctype(M)
[dummy,DD] = eig(M); %#ok<ASGLU>
EGV = DD([1 5 9]);
EGV = sort(EGV,'descend');
MM = mean(EGV);
EGV = EGV - MM;
if (abs(MM) + max(abs(EGV(1)),abs(EGV(3)))) == 0
  error('Tensor cannot be translated into [k,T] space.');
end
k = MM / (abs(MM) + max(abs(EGV(1)),abs(EGV(3))));
if max(abs(EGV(1)),abs(EGV(3))) == 0
  T = 0;
else
  T = 2 * EGV(2) / max(abs(EGV(1)),abs(EGV(3)));
end
end

function [UU,VV] = ae_uvt(TT,KK)

TAU = TT.*(1-abs(KK));
UU = NaN*ones(size(TT));
VV = UU;

% 2nd and 4th quadrants
II = (TAU > 0 & KK < 0) | (TAU < 0 & KK > 0);
UU(II) = TAU(II);
VV(II) = KK(II);

% First quadrant, Region A
II = (TAU < 4*KK) & (TAU >= 0 & KK >= 0);
UU(II) = TAU(II)./(1-TAU(II)/2);
VV(II) = KK(II)./(1-TAU(II)/2);

% First quadrant, Region B
II = (TAU >= 4*KK) & (TAU >= 0 & KK >= 0);
UU(II) = TAU(II)./(1-2*KK(II));
VV(II) = KK(II)./(1-2*KK(II));

% Third quadrant
II = (TAU >= 4*KK) & (TAU <= 0 & KK <= 0);
UU(II) = TAU(II)./(1+TAU(II)/2);
VV(II) = KK(II)./(1+TAU(II)/2);

II = (TAU < 4*KK) & (TAU <= 0 & KK <= 0);
UU(II) = TAU(II)./(1+2*KK(II));
VV(II) = KK(II)./(1+2*KK(II));

if sum(isnan(VV))
  error('Error plotting point in [uu,vv] space.');
end
end