%% run06e_compare_Himawari_SST_CHLA_PIV_student.m
% Compare and overlay MCC/PIV vectors estimated from Himawari SST and CHLA.
%
% Run order:
%   run06c_Himawari_SST_PIV_student
%   run06d_Himawari_CHLA_PIV_student
%   run06e_compare_Himawari_SST_CHLA_PIV_student
%
% This version reads the legacy output names, which are overwritten by the
% current PIV_MODE in run06c/d.  Therefore, if run06c/d were run in 'fine'
% mode, run06e compares the fine-mode vectors.

clear; close all; clc;

OUTPUT_ROOT = fullfile(pwd, 'out06_Himawari_PIV');
FIG_DIR = fullfile(OUTPUT_ROOT, 'figures');
DATA_DIR = fullfile(OUTPUT_ROOT, 'data');
ensure_dir(FIG_DIR);
ensure_dir(DATA_DIR);

sst_file  = fullfile(DATA_DIR, 'run06c_SST_PIV_vectors.mat');
chla_file = fullfile(DATA_DIR, 'run06d_CHLA_PIV_vectors.mat');

if ~exist(sst_file, 'file')
    error('Cannot find %s. Run run06c first.', sst_file);
end
if ~exist(chla_file, 'file')
    error('Cannot find %s. Run run06d first.', chla_file);
end

S = load(sst_file);
C = load(chla_file);

pivS = S.piv;
pivC = C.piv;

fprintf('=== run06e SST/CHLA PIV comparison ===\n');
fprintf('SST mode : %s, final vectors = %d\n', pivS.settings.piv_mode, nnz(pivS.mask_final));
fprintf('CHLA mode: %s, final vectors = %d\n', pivC.settings.piv_mode, nnz(pivC.mask_final));

ROI_LON = [147.0 149.8];
ROI_LAT = [37.2 39.9];
PLOT_ARROW_SCALE = 3.0;

%% Overlay figure
fig1 = figure('Color','w','Position',[80 80 1400 850]);
hold on;
plot_vectors_struct(pivS, pivS.mask_final, PLOT_ARROW_SCALE, 'k');
plot_vectors_struct(pivC, pivC.mask_final, PLOT_ARROW_SCALE, 'r');
axis tight;
set(gca,'YDir','normal');
xlabel('Longitude'); ylabel('Latitude');
title('run06e: SST-PIV black and CHLA-PIV red');
legend({'SST PIV','CHLA PIV'}, 'Location','best');
grid on;

%% ROI overlay figure
fig2 = figure('Color','w','Position',[80 80 1100 850]);
hold on;
maskSroi = pivS.mask_final & pivS.lon >= ROI_LON(1) & pivS.lon <= ROI_LON(2) & pivS.lat >= ROI_LAT(1) & pivS.lat <= ROI_LAT(2);
maskCroi = pivC.mask_final & pivC.lon >= ROI_LON(1) & pivC.lon <= ROI_LON(2) & pivC.lat >= ROI_LAT(1) & pivC.lat <= ROI_LAT(2);
plot_vectors_struct(pivS, maskSroi, PLOT_ARROW_SCALE, 'k');
plot_vectors_struct(pivC, maskCroi, PLOT_ARROW_SCALE, 'r');
xlim(ROI_LON); ylim(ROI_LAT);
set(gca,'YDir','normal');
xlabel('Longitude'); ylabel('Latitude');
title(sprintf('run06e ROI: SST black (%d), CHLA red (%d)', nnz(maskSroi), nnz(maskCroi)));
legend({'SST PIV','CHLA PIV'}, 'Location','best');
grid on;

%% Nearest-neighbor comparison at SST vector points
% For each SST vector, find the nearest CHLA vector within a radius.
MAX_DIST_DEG = 0.05;
[lonS, latS, US, VS, spS] = extract_vectors(pivS, pivS.mask_final);
[lonC, latC, UC, VC, spC] = extract_vectors(pivC, pivC.mask_final);

n = numel(lonS);
U_C_on_S = nan(n,1);
V_C_on_S = nan(n,1);
sp_C_on_S = nan(n,1);
dist_near = nan(n,1);

for i = 1:n
    d2 = (lonC - lonS(i)).^2 + (latC - latS(i)).^2;
    [dmin2, j] = min(d2);
    if isfinite(dmin2) && sqrt(dmin2) <= MAX_DIST_DEG
        U_C_on_S(i) = UC(j);
        V_C_on_S(i) = VC(j);
        sp_C_on_S(i) = spC(j);
        dist_near(i) = sqrt(dmin2);
    end
end

ok = isfinite(U_C_on_S) & isfinite(V_C_on_S);
angle_diff_deg = nan(n,1);
for i = 1:n
    if ok(i)
        dotuv = US(i)*U_C_on_S(i) + VS(i)*V_C_on_S(i);
        den = hypot(US(i),VS(i)) * hypot(U_C_on_S(i),V_C_on_S(i));
        if den > 0
            angle_diff_deg(i) = acosd(max(-1,min(1,dotuv/den)));
        end
    end
end
speed_diff = sp_C_on_S - spS;

fig3 = figure('Color','w','Position',[80 80 1400 600]);
tiledlayout(1,3,'TileSpacing','compact','Padding','compact');

nexttile;
scatter(spS(ok), sp_C_on_S(ok), 15, 'filled');
hold on; plot([0 1.5],[0 1.5],'k--');
xlabel('SST-PIV speed [m s^{-1}]'); ylabel('CHLA-PIV speed [m s^{-1}]');
title(sprintf('Nearest vectors within %.2f deg: n=%d', MAX_DIST_DEG, nnz(ok)));
grid on; axis square;

nexttile;
histogram(speed_diff(ok), 30);
xlabel('CHLA speed - SST speed [m s^{-1}]'); ylabel('Count');
title('Speed difference');
grid on;

nexttile;
histogram(angle_diff_deg(ok), 30);
xlabel('Angle difference [deg]'); ylabel('Count');
title('Direction difference');
grid on;

fig_bases = { ...
    'run06e_SST_CHLA_overlay', ...
    'run06e_SST_CHLA_overlay_ROI_148E_385N', ...
    'run06e_SST_CHLA_nearest_vector_comparison'};
figure_files = save_figure_set([fig1 fig2 fig3], fig_bases, FIG_DIR);

comparison = struct();
comparison.sst_file = sst_file;
comparison.chla_file = chla_file;
comparison.max_dist_deg = MAX_DIST_DEG;
comparison.lonS = lonS;
comparison.latS = latS;
comparison.US = US;
comparison.VS = VS;
comparison.speedS = spS;
comparison.UC_on_S = U_C_on_S;
comparison.VC_on_S = V_C_on_S;
comparison.speedC_on_S = sp_C_on_S;
comparison.dist_near = dist_near;
comparison.angle_diff_deg = angle_diff_deg;
comparison.speed_diff = speed_diff;
comparison.ok = ok;
comparison.figure_files = figure_files;

save(fullfile(DATA_DIR, 'run06e_SST_CHLA_PIV_comparison.mat'), 'comparison', 'pivS', 'pivC');
fprintf('Saved comparison MAT: %s\n', fullfile(DATA_DIR, 'run06e_SST_CHLA_PIV_comparison.mat'));
fprintf('Saved figures in    : %s\n', FIG_DIR);

%% local functions
function ensure_dir(dirname)
    if ~exist(dirname, 'dir')
        mkdir(dirname);
    end
end

function figure_files = save_figure_set(figHandles, baseNames, figDir)
    if numel(figHandles) ~= numel(baseNames)
        error('figHandles and baseNames must have the same length.');
    end
    ensure_dir(figDir);
    figure_files = struct('fig', cell(numel(figHandles),1), 'png', cell(numel(figHandles),1));
    for k = 1:numel(figHandles)
        fig_file = fullfile(figDir, [baseNames{k} '.fig']);
        png_file = fullfile(figDir, [baseNames{k} '.png']);
        if ~ishandle(figHandles(k)) || ~isvalid(figHandles(k))
            continue;
        end
        set(figHandles(k), 'PaperPositionMode', 'auto');
        savefig(figHandles(k), fig_file);
        print(figHandles(k), png_file, '-dpng', '-r200');
        figure_files(k).fig = fig_file;
        figure_files(k).png = png_file;
        fprintf('Saved figure        : %s\n', png_file);
    end
end

function plot_vectors_struct(piv, mask, scaleFactor, colorSpec)
    qx = piv.dx_grid .* piv.grid.dlon .* scaleFactor;
    qy = piv.dy_grid .* piv.grid.dlat .* scaleFactor;
    quiver(piv.lon(mask), piv.lat(mask), qx(mask), qy(mask), 0, colorSpec, 'LineWidth', 1.0, 'MaxHeadSize', 1.8);
end

function [lon, lat, U, V, sp] = extract_vectors(piv, mask)
    lon = piv.lon(mask);
    lat = piv.lat(mask);
    U = piv.U(mask);
    V = piv.V(mask);
    sp = piv.speed(mask);
    lon = lon(:); lat = lat(:); U = U(:); V = V(:); sp = sp(:);
end
