% 原连续域(浅灰) + 新增连续域(浅红) + 原始光锥(蓝虚线) + 折叠段(黑实线)
clear; clc;

% 周期/折射率
ax = 0.5; ay = 1.0; n = 1.4;
sz = 10;
% 高对称点（绝对 k）
Gpt = [0,     0    ];
Xpt = [pi/ax, 0    ];
Ypt = [0,     pi/ay];
Mpt = [pi/ax, pi/ay];

% 路径按长度比例采样：M–X, X–Γ, Γ–Y → 总 Ntot 点
L = @(A,B) norm(A-B);
L_MX = L(Mpt,Xpt);  L_XG = L(Xpt,Gpt);  L_GY = L(Gpt,Ypt);
L_tot = L_MX + L_XG + L_GY;
Ntot = 800;
n_MX = round(Ntot * L_MX / L_tot);
n_XG = round(Ntot * L_XG / L_tot);
n_GY = Ntot - n_MX - n_XG;

seg = @(A,B,N) [linspace(A(1),B(1),N)', linspace(A(2),B(2),N)'];
path_MX = seg(Mpt, Xpt, n_MX);
path_XG = seg(Xpt, Gpt, n_XG);
path_GY = seg(Gpt, Ypt, n_GY);
path = [path_MX; path_XG; path_GY];   % (Ntot x 2)

% 原始 light line（纵轴归一到 a_y）
kabs = vecnorm(path,2,2);
f = (ay/(2*pi*n)) * kabs;

% -------- 折叠段1：X–M 折到 Γ–Y --------
idxGY   = (n_MX + n_XG + 1) : (n_MX + n_XG + n_GY);   % 目标 Γ→Y
f_XM    = f(1:n_MX);                                   % 源 X–M（数据里顺序是 M→X）
f_newGY = flipud(f_XM);                                 % 折叠到 Γ→Y

% -------- 折叠段2：X/2–X 折到 Γ–X/2 --------
len_half   = floor(n_XG/2);
idx_srcGX2 = (n_MX + 1) : (n_MX + len_half);           % 源：X→X/2
idx_dstGX2 = (n_MX + len_half + 1) : (n_MX + n_XG);    % 目标：X/2→Γ
f_srcGX2   = f(idx_srcGX2);
f_newGX2   = flipud(f_srcGX2);                          % 折叠到 Γ→X/2（升序）

% -------- 画图 --------
figure; hold on;

% y 上限（统一）
ytop = max([f; f_newGY; f_newGX2]) * 1.05;

% (A) 原连续域：整条路径上方（原光线 f 以上）涂浅灰
fill([1:Ntot, fliplr(1:Ntot)], ...
     [f',      ytop*ones(1,Ntot)], ...
     [1.00 0.85 0.85], 'EdgeColor','none', 'FaceAlpha',0.55);

% (B) 新增连续域：在两段上，把“新光线 f_new”与“原光线 f”之间的带状区域涂浅红
%   段1：Γ–Y（x = idxGY）
x = idxGY; y_low = f_newGY'; y_up = f(idxGY)';           % 新增 = [f_new, f_orig] 之间
fill([x, fliplr(x)], [y_low, fliplr(y_up)], ...
     [0.90 0.90 0.90], 'EdgeColor','none', 'FaceAlpha',0.70);

%   段2：Γ–X/2（x = idx_dstGX2）
x = idx_dstGX2; y_low = f_newGX2'; y_up = f(idx_dstGX2)';
fill([x, fliplr(x)], [y_low, fliplr(y_up)], ...
     [0.90 0.90 0.90], 'EdgeColor','none', 'FaceAlpha',0.70);

% (C) 线条：原光锥（蓝虚线）+ 两段折叠（黑实线）
plot(f, 'k-', 'LineWidth', 1.6);                  % 原始 light line（全程虚线）
plot(idxGY,      f_newGY,  'k--', 'LineWidth', 2);  % 折叠段1：Γ–Y（黑实）
plot(idx_dstGX2, f_newGX2, 'k--', 'LineWidth', 2);  % 折叠段2：Γ–X/2（黑实）

% 分段竖线 & 标签
xline(n_MX,'k--'); xline(n_MX+n_XG,'k--');xline(n_MX+n_XG/2,'k--');
xticks([1, n_MX, n_MX+n_XG, Ntot]);
xticklabels({'M','X','X1','\Gamma','Y'});

% xlabel('Path index (M–X–\Gamma–Y)');
% ylabel('\omega a_y / (2\pi c)');
% title('原连续域(浅灰) + 新增连续域(浅红) + 原始光锥(虚线) + 折叠段(实线)');
grid off; box on;
% legend({'原连续域','新增连续域(折叠)','新增连续域(折叠)', ...
        % '原光锥','折叠：\Gamma–Y','折叠：\Gamma–X/2'}, 'Location','northwest');
ylim([0 ytop]);
hold on


%% 1. 载入能带和Q数据
a = load('bands_original.mat');
a1 = load('bands_folding.mat');
k1 = a1.kpoints/pi;
k = a.kpoints/pi;   % k 向量 (单位 π/a_x)
freqs = a.freqs;
Qs = a.Qs;

[nk, nm] = size(freqs);

%% 2. 设置晶格 & 高对称点采样
ax = 1; ay = 2*ax;      % y方向周期是x方向的2倍
bx = 2*pi/ax; by = 2*pi/ay;

% 高对称点（绝对 k）
M = [pi/ax, pi/ay];
X = [pi/ax, 0];
X1 = [pi/2*ax,0];
G = [0, 0];
Y = [0, pi/ay];

% 计算各段长度
len_MX = norm(M - X);
len_XG = norm(X - G);
len_GY = norm(G - Y);
len_total = len_MX + len_XG + len_GY;

% 按长度比例分配点数，总点数= nk
n_MX = round(nk * len_MX / len_total);
n_XG = round(nk * len_XG / len_total);
n_GY = nk - n_MX - n_XG;   % 剩下的点数给最后一段

% 构造路径（绝对k）
seg = @(A,B,N) [linspace(A(1),B(1),N)', linspace(A(2),B(2),N)'];
path_MX = seg(M, X, n_MX);
path_XG = seg(X, G, n_XG);
path_GY = seg(G, Y, n_GY);
path = [path_MX; path_XG; path_GY];

%% 3. 折叠后的 light line 计算（可选仅x折叠）
norm_by = 'ay';           % 纵轴归一：'ax' 或 'ay'
n_use   = 1.4;              % 折射率（空气）
scale   = ax; if strcmp(norm_by,'ay'), scale = ay; end

only_x_fold = false;      % <<<<<< 只做 x 方向折叠就设 true

if only_x_fold
    ms = -1:1; ns = 0;    % 只沿 x 取相邻平移
else
    ms = -1:1; ns = -1:1; % x,y 都取相邻平移
end

Gs = [];
for m = ms
    for n = ns
        Gs = [Gs; m*bx, n*by];
    end
end

% 原始 light line 边界（未折叠）
kabs = vecnorm(path, 2, 2);
f_orig = (scale/(2*pi*n_use)) * kabs;      % ω a_norm /(2πc)

% 折叠后包络（取 min |k+G|）
f_all = zeros(nk, size(Gs,1));
for i = 1:nk
    kabs_list = vecnorm((path(i,:) + Gs), 2, 2);
    f_all(i,:) = (scale/(2*pi*n_use)) * kabs_list;
end
f_fold = min(f_all, [], 2);                % 折叠后连续域边界

%% 4. Q 值处理（和你原来一样）
Qs_clipped = Qs;
Qs_clipped(Qs > 1e15) = NaN;

%% 5. 绘图（原连续域=浅灰；新增连续域=浅红；散点=log10(Q)）
 hold on;

% 先定 Y 上限（不然填充会被后面更改）
f_max = max([freqs(:); f_orig; f_fold])*1.05;
ylim([0.1 0.5]);  % 你原本的限制，如需自适应可注释掉
yl = ylim;      % 取最终 Y 轴范围

% (A) 原连续域：f_orig 以上整体浅灰
% fill([1:nk, fliplr(1:nk)], [f_orig', yl(2)*ones(1,nk)], ...
%      [0.90 0.90 0.90], 'FaceAlpha', 0.55, 'EdgeColor', 'none');

% (B) 新增连续域：在 f_fold < f_orig 的区域，填 f_fold 与 f_orig 之间的“带”（浅红）
mask = (f_fold < f_orig);
% 找出 mask 的连续片段
dmask = diff([false; mask; false]);
run_starts = find(dmask == 1);
run_ends   = find(dmask == -1) - 1;

for r = 1:numel(run_starts)
    i1 = run_starts(r); i2 = run_ends(r);
    x = i1:i2;
    y_low = f_fold(x)'; y_up = f_orig(x)';
    fill([x, fliplr(x)], [y_low, fliplr(y_up)], ...
         [1.00 0.85 0.85], 'FaceAlpha', 0.70, 'EdgeColor', 'none'); % 浅红
end

% (C) Q 色散点（log10 颜色）
for m = 1:nm
    scatter(1:nk, freqs(:,m), sz, log10(Qs_clipped(:,m)), 'filled');
    mask_inf = Qs(:,m) > 1e15;
    scatter(find(mask_inf), freqs(mask_inf,m), sz, 'y', 'filled');
end

% (D) 画边界线：原始（蓝虚线） + 折叠后包络（黑实线）
plot(1:nk, f_orig, 'b--', 'LineWidth', 1.2);
plot(1:nk, f_fold, 'k-',  'LineWidth', 1.6);

% 轴与标签
xline(n_MX, 'k--');
xline(n_MX + n_XG, 'k--');
xline(n_MX + n_XG/2, 'b--');
xticks([1, n_MX, n_MX + n_XG/2,  n_MX + n_XG, nk]);
xticklabels({'M','X','X1','\Gamma','Y'});
% ylabel(sprintf('\\omega a_{%s} / (2\\pi c)', norm_by));
xlim([0 nk+1]);
grid off; box on;

% 颜色映射（你原本使用的）
colormap(slanCM('plasma'));
cb = colorbar; cb.Label.String = 'log_{10}(Q)';

hold on
%% 
%% 1. 载入能带和Q数据

a= load('bands_folding.mat');

k = a.kpoints/pi;   % k 向量 (单位 π/a_x)
freqs = a.freqs;
Qs = a.Qs;

[nk, nm] = size(freqs);

%% 2. 设置晶格 & 高对称点采样
ax = 1; ay = 2*ax;      % y方向周期是x方向的2倍
bx = 2*pi/ax; by = 2*pi/ay;

% 高对称点（绝对 k）
M = [pi/ax, pi/ay];
X = [pi/ax, 0];
X1 = [pi/2*ax,0];
G = [0, 0];
Y = [0, pi/ay];

% 计算各段长度
len_MX = 0;
len_XG = norm(X/2 - G);
len_GY = norm(G - Y);
len_total = len_MX + len_XG + len_GY;

% 按长度比例分配点数，总点数= nk
n_MX = round(nk * len_MX / len_total);
n_XG = round(nk * len_XG / len_total);
n_GY = nk - n_MX - n_XG;   % 剩下的点数给最后一段

% 构造路径（绝对k）
seg = @(A,B,N) [linspace(A(1),B(1),N)', linspace(A(2),B(2),N)'];
path_MX = seg(M, X, n_MX);
path_XG = seg(X, G, n_XG);
path_GY = seg(G, Y, n_GY);
path = [path_MX; path_XG; path_GY];

%% 3. 折叠后的 light line 计算（可选仅x折叠）
norm_by = 'ay';           % 纵轴归一：'ax' 或 'ay'
n_use   = 1;              % 折射率（空气）
scale   = ax; if strcmp(norm_by,'ay'), scale = ay; end

only_x_fold = false;      % <<<<<< 只做 x 方向折叠就设 true

if only_x_fold
    ms = -1:1; ns = 0;    % 只沿 x 取相邻平移
else
    ms = -1:1; ns = -1:1; % x,y 都取相邻平移
end

Gs = [];
for m = ms
    for n = ns
        Gs = [Gs; m*bx, n*by];
    end
end

% 原始 light line 边界（未折叠）
kabs = vecnorm(path, 2, 2);
f_orig = (scale/(2*pi*n_use)) * kabs;      % ω a_norm /(2πc)

% 折叠后包络（取 min |k+G|）
f_all = zeros(nk, size(Gs,1));
for i = 1:nk
    kabs_list = vecnorm((path(i,:) + Gs), 2, 2);
    f_all(i,:) = (scale/(2*pi*n_use)) * kabs_list;
end
f_fold = min(f_all, [], 2);                % 折叠后连续域边界

%% 4. Q 值处理（和你原来一样）
Qs_clipped = Qs;
Qs_clipped(Qs > 1e15) = NaN;

%% 5. 绘图（原连续域=浅灰；新增连续域=浅红；散点=log10(Q)）
 hold on;

% 先定 Y 上限（不然填充会被后面更改）
f_max = max([freqs(:); f_orig; f_fold])*1.05;
ylim([0 0.5]);  % 你原本的限制，如需自适应可注释掉
yl = ylim;      % 取最终 Y 轴范围

% (A) 原连续域：f_orig 以上整体浅灰
% fill([1:nk, fliplr(1:nk)], [f_orig', yl(2)*ones(1,nk)], ...
%      [0.90 0.90 0.90], 'FaceAlpha', 0.55, 'EdgeColor', 'none');

% (B) 新增连续域：在 f_fold < f_orig 的区域，填 f_fold 与 f_orig 之间的“带”（浅红）
mask = (f_fold < f_orig);
% 找出 mask 的连续片段
dmask = diff([false; mask; false]);
run_starts = find(dmask == 1);
run_ends   = find(dmask == -1) - 1;

for r = 1:numel(run_starts)
    i1 = run_starts(r); i2 = run_ends(r);
    x = i1:i2;
    y_low = f_fold(x)'; y_up = f_orig(x)';
    fill([x+399, fliplr(x)], [y_low, fliplr(y_up)], ...
         [1.00 0.85 0.85], 'FaceAlpha', 0.70, 'EdgeColor', 'none'); % 浅红
end

% (C) Q 色散点（log10 颜色）
for m = 1:nm
    scatter(400:(399+nk), freqs(:,m), sz, log10(Qs_clipped(:,m)), 'filled');
    mask_inf = Qs(:,m) > 1e15;
    scatter(find(mask_inf)+399, freqs(mask_inf,m), sz, 'y', 'filled');
end

% % (D) 画边界线：原始（蓝虚线） + 折叠后包络（黑实线）
% plot(1:nk, f_orig, 'b--', 'LineWidth', 1.2);
% plot(1:nk, f_fold, 'k-',  'LineWidth', 1.6);

% 轴与标签


% xlim([0 nk+1]);
grid off; box on;

% 颜色映射（你原本使用的）
colormap(slanCM('plasma'));
cb = colorbar; cb.Label.String = 'log_{10}(Q)';


