SLAから地衡流・渦を診断する

%% run05a_geostrophic_current_central_student.m
% ============================================================
% Ocean Remote Sensing Practical
% Step 1: Calculate geostrophic current from SLA/ADT
%         using central differences
%
% Purpose:
%   - Read SLA or ADT from satellite altimetry NetCDF
%   - Calculate horizontal gradients on lon-lat grid
%   - Calculate geostrophic velocity:
%
%       u_g = -(g/f) * d eta / dy
%       v_g =  (g/f) * d eta / dx
%
%   - Save results for the next step
%
% Notes:
%   eta = SLA : geostrophic velocity anomaly
%   eta = ADT : absolute geostrophic velocity
%
% Required:
%   MATLAB only
% ============================================================

clear; close all; clc;

%% ===================== User settings =====================

ncfile = 'nrt_global_allsat_phy_l4_20240602_20240608.nc';

% Display / analysis region
lonlim = [120 170];
latlim = [20 55];

% false: use SLA
% true : use ADT
useAbsolute = false;

% Output directory
outdir = 'run05a_geostrophic_output';
if ~exist(outdir, 'dir')
    mkdir(outdir);
end

%% ===================== Constants =====================

g = 9.81;                 % gravitational acceleration [m s^-2]
Omega = 7.2921159e-5;     % Earth's rotation rate [s^-1]
Re = 6371000;             % Earth radius [m]

%% ===================== Read data =====================

if ~isfile(ncfile)
    error('NetCDF file not found: %s', ncfile);
end

lon = double(ncread(ncfile, 'longitude'));
lat = double(ncread(ncfile, 'latitude'));

ilon = find(lon >= lonlim(1) & lon <= lonlim(2));
ilat = find(lat >= latlim(1) & lat <= latlim(2));

lonSub = lon(ilon);
latSub = lat(ilat);

[Lon, Lat] = meshgrid(lonSub, latSub);

if useAbsolute
    etaName = 'adt';
    etaLabel = 'ADT';
else
    etaName = 'sla';
    etaLabel = 'SLA';
end

eta = read_cmems_var_latlon(ncfile, etaName, lon, lat, ilon, ilat);

% Remove ice-covered points if flag_ice exists
try
    ice = read_cmems_var_latlon(ncfile, 'flag_ice', lon, lat, ilon, ilat);
    eta(ice == 1) = NaN;
catch
end

% Remove unrealistic values
eta(abs(eta) > 5) = NaN;

fprintf('Loaded %s from %s\n', etaName, ncfile);
fprintf('eta size = [%d %d] = [lat lon]\n', size(eta,1), size(eta,2));

%% ===================== Grid spacing =====================

% lon/lat are in degrees.
% Convert degree differences to meter distances.
%
% dx depends on latitude:
%   dx = R cos(lat) dlon
%
% dy is almost independent of latitude:
%   dy = R dlat

nlat = numel(latSub);
nlon = numel(lonSub);

deta_dx = NaN(nlat, nlon);
deta_dy = NaN(nlat, nlon);

%% ===================== Central difference =====================

% x-direction gradient: d eta / dx
for j = 2:nlon-1

    dlon_rad = deg2rad(__________________);

    for i = 1:nlat

        dx = ________ * cosd(________) * __________;

        deta_dx(i,j) = (____________ - ____________) / ____;
    end
end

% y-direction gradient: d eta / dy
for i = 2:nlat-1

    dlat_rad = deg2rad(__________________);
    dy = ________ * __________;

    for j = 1:nlon

        deta_dy(i,j) = (____________ - ____________) / ____;
    end
end

% Edges cannot be calculated by central difference.
% Keep them as NaN.
deta_dx(:,1)   = NaN;
deta_dx(:,end) = NaN;
deta_dy(1,:)   = NaN;
deta_dy(end,:) = NaN;

%% ===================== Geostrophic current =====================

f = 2 * ________ * sind(_____);

ug = -(____ ./ ____) .* __________;
vg =  (____ ./ ____) .* __________;

% Avoid near-equator singularity
ug(abs(f) < 1e-8) = NaN;
vg(abs(f) < 1e-8) = NaN;

speed = sqrt(______.^2 + ______.^2);

% Remove unrealistic speeds
ug(abs(ug) > 5) = NaN;
vg(abs(vg) > 5) = NaN;
speed(abs(speed) > 5) = NaN;

%% ===================== Scaling check =====================

fprintf('\nScaling check:\n');
fprintf('  If eta changes by 0.1 m over 100 km,\n');
fprintf('  slope is about 1e-6.\n');
fprintf('  g/f is about 1e5 at mid-latitudes.\n');
fprintf('  geostrophic speed is about 0.1 m/s.\n');

fprintf('\nComputed geostrophic current:\n');
fprintf('  max |ug| = %.3f m/s\n', max(abs(ug(:)), [], 'omitnan'));
fprintf('  max |vg| = %.3f m/s\n', max(abs(vg(:)), [], 'omitnan'));
fprintf('  max speed = %.3f m/s\n', max(speed(:), [], 'omitnan'));

%% ===================== Figure =====================

figure('Color','w','Position',[100 100 1500 650]);
tiledlayout(1,2,'TileSpacing','compact','Padding','compact');

nexttile;
imagesc(lonSub, latSub, eta);
set(gca,'YDir','normal');
axis tight;
grid on; box on;
colorbar;
caxis([-1 1]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
title([etaLabel ' and contours']);
xlabel('Longitude');
ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, speed);
set(gca,'YDir','normal');
axis tight;
grid on; box on;
colorbar;
caxis([0 1.5]);
hold on;
skip = 4;
quiver(Lon(1:skip:end,1:skip:end), Lat(1:skip:end,1:skip:end), ...
       ug(1:skip:end,1:skip:end), vg(1:skip:end,1:skip:end), ...
       1.5, 'k');
title('Geostrophic speed and vectors from central difference');
xlabel('Longitude');
ylabel('Latitude');

sgtitle('Step 1: Geostrophic current from SLA/ADT using central difference');

exportgraphics(gcf, fullfile(outdir, 'Fig01_geostrophic_current_central.png'), 'Resolution', 180);

%% ===================== Save =====================

outmat = fullfile(outdir, 'geostrophic_current_central.mat');

save(outmat, ...
    'ncfile', 'lonSub', 'latSub', 'Lon', 'Lat', ...
    'eta', 'etaName', 'etaLabel', ...
    'deta_dx', 'deta_dy', ...
    'ug', 'vg', 'speed', 'f', ...
    'useAbsolute', 'g', 'Omega', 'Re', ...
    '-v7.3');

fprintf('\nSaved: %s\n', outmat);
disp('Done.');

%% ===================== Local function =====================

function A = read_cmems_var_latlon(ncfile, varname, lonFull, latFull, ilon, ilat)

    raw = ncread(ncfile, varname);
    raw = double(squeeze(raw));

    info = ncinfo(ncfile, varname);

    fillValue = [];
    scaleFactor = [];
    addOffset = 0;

    for ia = 1:numel(info.Attributes)
        aname = info.Attributes(ia).Name;
        aval  = info.Attributes(ia).Value;

        if strcmp(aname, '_FillValue')
            fillValue = double(aval);
        elseif strcmp(aname, 'scale_factor')
            scaleFactor = double(aval);
        elseif strcmp(aname, 'add_offset')
            addOffset = double(aval);
        end
    end

    if ~isempty(fillValue)
        raw(raw == fillValue) = NaN;
    end

    finiteVals = raw(isfinite(raw));

    if ~isempty(scaleFactor) && ~isempty(finiteVals)
        if prctile(abs(finiteVals), 99) > 50
            raw = raw * scaleFactor + addOffset;
        end
    end

    nlon = numel(lonFull);
    nlat = numel(latFull);

    if size(raw,1) == nlon && size(raw,2) == nlat
        Aall = raw.';       % [lat lon]
    elseif size(raw,1) == nlat && size(raw,2) == nlon
        Aall = raw;         % already [lat lon]
    else
        error('Unexpected size for %s', varname);
    end

    A = Aall(ilat, ilon);
end

%% run05b_OW_EKE_vorticity_student.m
% ============================================================
% Ocean Remote Sensing Practical
% Step 2: Calculate relative vorticity, strain, OWrot, and KE/EKE-like
%
% Input:
%   run05a_geostrophic_output/geostrophic_current_central.mat
%
% Output:
%   run05b_diagnostics_output/OW_EKE_vorticity.mat
%
% Notes:
%   If eta is SLA, ug/vg are geostrophic velocity anomalies.
%   Then 0.5*(ug^2+vg^2) can be interpreted as EKE-like.
%
%   If eta is ADT, ug/vg are closer to absolute geostrophic velocity.
%   Then 0.5*(ug^2+vg^2) is KE, not EKE.
% ============================================================

clear; close all; clc;

%% ===================== Settings =====================

inmat = fullfile('run05a_geostrophic_output', 'geostrophic_current_central.mat');

outdir = 'run05b_diagnostics_output';
if ~exist(outdir, 'dir')
    mkdir(outdir);
end

% Smoothing window for velocity before differentiating
% 1 = no smoothing
% 3 or 5 is useful because derivatives amplify noise
smoothWindow = 5;

% Vector skip for plotting
vectorSkip = 4;

%% ===================== Load =====================

if ~isfile(inmat)
    error('Input MAT file not found: %s', inmat);
end

S = load(inmat);

lonSub = S.lonSub;
latSub = S.latSub;
Lon = S.Lon;
Lat = S.Lat;

eta = S.eta;
ug = S.ug;
vg = S.vg;
f = S.f;

etaLabel = S.etaLabel;
useAbsolute = S.useAbsolute;

fprintf('Loaded: %s\n', inmat);
fprintf('Field: %s\n', etaLabel);

%% ===================== Smooth velocity =====================

ug_s = smooth2_nan(ug, smoothWindow);
vg_s = smooth2_nan(vg, smoothWindow);

speed = sqrt(______.^2 + ______.^2);

%% ===================== Velocity gradients =====================

[du_dx, du_dy] = central_gradient_lonlat_m(ug_s, lonSub, latSub);
[dv_dx, dv_dy] = central_gradient_lonlat_m(vg_s, lonSub, latSub);

%% ===================== Relative vorticity and strain =====================

% Relative vorticity
zeta = ________ - ________;

% Strain components
sn = ________ - ________;     % normal strain
ss = ________ + ________;     % shear strain

strain2 = _____.^2 + _____.^2;
vort2   = ______.^2;

% Conventional Okubo-Weiss:
%   OW = strain^2 - vorticity^2
%   OW < 0 means rotation-dominant.
OW = ________ - ________;

% Teaching-friendly sign:
%   OWrot = vorticity^2 - strain^2
%   OWrot > 0 means rotation-dominant.
OWrot = ________ - ________;

% Relative vorticity normalized by planetary vorticity
zetaOverF = ______ ./ ____;
zetaOverF(abs(f) < 1e-8) = NaN;

%% ===================== KE / EKE-like =====================

KE = 0.5 * (______.^2 + ______.^2);

if useAbsolute
    KElabel = 'Geostrophic KE';
else
    KElabel = 'EKE-like from SLA-derived velocity anomaly';
end

%% ===================== Common mask =====================

valid = isfinite(eta) & isfinite(ug_s) & isfinite(vg_s);

speed(~valid) = NaN;
zeta(~valid) = NaN;
zetaOverF(~valid) = NaN;
sn(~valid) = NaN;
ss(~valid) = NaN;
strain2(~valid) = NaN;
vort2(~valid) = NaN;
OW(~valid) = NaN;
OWrot(~valid) = NaN;
KE(~valid) = NaN;

%% ===================== Print summary =====================

fprintf('\nDiagnostics summary:\n');
fprintf('  max speed = %.3f m/s\n', max(speed(:), [], 'omitnan'));
fprintf('  max KE/EKE-like = %.4f m^2/s^2\n', max(KE(:), [], 'omitnan'));
fprintf('  max |zeta/f| = %.3f\n', max(abs(zetaOverF(:)), [], 'omitnan'));
fprintf('  max OWrot = %.3e s^-2\n', max(OWrot(:), [], 'omitnan'));
fprintf('  min OWrot = %.3e s^-2\n', min(OWrot(:), [], 'omitnan'));

%% ===================== Figure =====================

figure('Color','w','Position',[80 80 1500 950]);
tiledlayout(2,2,'TileSpacing','compact','Padding','compact');

nexttile;
imagesc(lonSub, latSub, speed);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colorbar;
caxis([0 1.5]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
quiver(Lon(1:vectorSkip:end,1:vectorSkip:end), Lat(1:vectorSkip:end,1:vectorSkip:end), ...
       ug_s(1:vectorSkip:end,1:vectorSkip:end), vg_s(1:vectorSkip:end,1:vectorSkip:end), ...
       1.5, 'k');
title('Geostrophic speed');
xlabel('Longitude'); ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, KE);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colorbar;
caxis([0 0.3]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
title(KElabel);
xlabel('Longitude'); ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, zetaOverF);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colormap(gca, redblue_colormap(256));
colorbar;
caxis([-0.5 0.5]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
title('Relative vorticity: \zeta / f');
xlabel('Longitude'); ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, OWrot * 1e10);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colormap(gca, redblue_colormap(256));
colorbar;
caxis([-5 5]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
contour(Lon, Lat, OWrot, [0 0], 'm', 'LineWidth', 1.0);
title('OW_{rot} = \zeta^2 - strain^2; magenta: OW_{rot}=0');
xlabel('Longitude'); ylabel('Latitude');

sgtitle('Step 2: OW, KE/EKE-like, and relative vorticity');

exportgraphics(gcf, fullfile(outdir, 'Fig01_OW_EKE_vorticity.png'), 'Resolution', 180);

%% ===================== Save =====================

outmat = fullfile(outdir, 'OW_EKE_vorticity.mat');

save(outmat, ...
    'lonSub', 'latSub', 'Lon', 'Lat', ...
    'eta', 'etaLabel', 'useAbsolute', ...
    'ug', 'vg', 'ug_s', 'vg_s', 'speed', ...
    'du_dx', 'du_dy', 'dv_dx', 'dv_dy', ...
    'zeta', 'zetaOverF', 'sn', 'ss', ...
    'strain2', 'vort2', 'OW', 'OWrot', ...
    'KE', 'KElabel', 'smoothWindow', ...
    '-v7.3');

fprintf('\nSaved: %s\n', outmat);
disp('Done.');

%% ===================== Local functions =====================

function B = smooth2_nan(A, win)

    if win <= 1
        B = A;
        return
    end

    if mod(win,2) == 0
        win = win + 1;
    end

    K = ones(win, win);
    M = isfinite(A);

    A0 = A;
    A0(~M) = 0;

    num = conv2(A0, K, 'same');
    den = conv2(double(M), K, 'same');

    B = num ./ den;
    B(den == 0) = NaN;
end

function [dA_dx, dA_dy] = central_gradient_lonlat_m(A, lon, lat)

    Re = 6371000;

    nlat = numel(lat);
    nlon = numel(lon);

    dA_dx = NaN(nlat, nlon);
    dA_dy = NaN(nlat, nlon);

    for j = 2:nlon-1

        dlon_rad = deg2rad(lon(j+1) - lon(j-1));

        for i = 1:nlat

            dx = Re * cosd(lat(i)) * dlon_rad;

            dA_dx(i,j) = (A(i,j+1) - A(i,j-1)) / dx;
        end
    end

    for i = 2:nlat-1

        dlat_rad = deg2rad(lat(i+1) - lat(i-1));
        dy = Re * dlat_rad;

        for j = 1:nlon

            dA_dy(i,j) = (A(i+1,j) - A(i-1,j)) / dy;
        end
    end
end

function cmap = redblue_colormap(n)

    if nargin < 1
        n = 256;
    end

    n1 = floor(n/2);
    n2 = n - n1;

    blue = [linspace(0,1,n1)' linspace(0,1,n1)' ones(n1,1)];
    red  = [ones(n2,1) linspace(1,0,n2)' linspace(1,0,n2)'];

    cmap = [blue; red];
end

%% run05c_SLA_eddy_interpretation_student.m
% ============================================================
% Ocean Remote Sensing Practical
% Step 3: Eddy interpretation using SLA, vectors, zeta/f, OWrot
%
% Input:
%   run05b_diagnostics_output/OW_EKE_vorticity.mat
%
% Purpose:
%   - Compare SLA/ADT, velocity vectors, relative vorticity, OWrot
%   - Identify eddy-like structures by visual interpretation
%
% Important:
%   OWrot > 0 does NOT automatically mean "eddy".
%   Eddy-like structures should be judged using multiple diagnostics:
%
%     1. closed SLA/ADT anomaly
%     2. rotating velocity vectors
%     3. coherent zeta/f patch
%     4. rotation-dominant OWrot area
%
% ============================================================

clear; close all; clc;

%% ===================== Settings =====================

inmat = fullfile('run05b_diagnostics_output', 'OW_EKE_vorticity.mat');

outdir = 'run05c_eddy_interpretation_output';
if ~exist(outdir, 'dir')
    mkdir(outdir);
end

% Local anomaly background smoothing window
bgWindow = 21;

% Vector skip
vectorSkip = 4;

% Candidate support thresholds
% These are not universal. They are only for guided interpretation.
localSLA_thr = 0.05;       % m
zetaF_thr    = 0.05;       % nondimensional
OWrot_thr    = 0;          % rotation-dominant

%% ===================== Load =====================

if ~isfile(inmat)
    error('Input MAT file not found: %s', inmat);
end

S = load(inmat);

lonSub = S.lonSub;
latSub = S.latSub;
Lon = S.Lon;
Lat = S.Lat;

eta = S.eta;
etaLabel = S.etaLabel;

ug_s = S.ug_s;
vg_s = S.vg_s;
speed = S.speed;

zetaOverF = S.zetaOverF;
OWrot = S.OWrot;
KE = S.KE;
KElabel = S.KElabel;

fprintf('Loaded: %s\n', inmat);
fprintf('Field: %s\n', etaLabel);

%% ===================== Local SLA anomaly =====================

% A broad background is removed to emphasize mesoscale structures.
eta_bg = smooth2_nan(eta, bgWindow);
etaLocal = eta - eta_bg;

%% ===================== Candidate support masks =====================

% These masks are not final eddy detection.
% They only indicate places where eddy-like interpretation may be supported.

mask_localSLA = abs(__________) >= _____________;
mask_zeta     = abs(__________) >= _________;
mask_OWrot    = ________ > _________;

candidateSupport = ______________ & _________ & __________;

% Simple score: how many conditions are satisfied?
eddyScore = double(______________) + double(_________) + double(__________);

eddyScore(~isfinite(etaLocal)) = NaN;

%% ===================== Figure 1: Interpretation panels =====================

figure('Color','w','Position',[80 80 1500 950]);
tiledlayout(2,2,'TileSpacing','compact','Padding','compact');

nexttile;
imagesc(lonSub, latSub, eta);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colorbar;
caxis([-1 1]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
quiver(Lon(1:vectorSkip:end,1:vectorSkip:end), Lat(1:vectorSkip:end,1:vectorSkip:end), ...
       ug_s(1:vectorSkip:end,1:vectorSkip:end), vg_s(1:vectorSkip:end,1:vectorSkip:end), ...
       1.5, 'k');
title([etaLabel ' + geostrophic vectors']);
xlabel('Longitude'); ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, etaLocal);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colormap(gca, redblue_colormap(256));
colorbar;
caxis([-0.4 0.4]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
title(['Local ' etaLabel ' anomaly']);
xlabel('Longitude'); ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, zetaOverF);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colormap(gca, redblue_colormap(256));
colorbar;
caxis([-0.5 0.5]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
title('\zeta / f');
xlabel('Longitude'); ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, OWrot * 1e10);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colormap(gca, redblue_colormap(256));
colorbar;
caxis([-5 5]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
contour(Lon, Lat, OWrot, [0 0], 'm', 'LineWidth', 1.0);
title('OW_{rot}; magenta: OW_{rot}=0');
xlabel('Longitude'); ylabel('Latitude');

sgtitle('Step 3: Eddy interpretation using multiple diagnostics');

exportgraphics(gcf, fullfile(outdir, 'Fig01_eddy_interpretation_panels.png'), 'Resolution', 180);

%% ===================== Figure 2: Candidate support =====================

figure('Color','w','Position',[100 100 1500 700]);
tiledlayout(1,2,'TileSpacing','compact','Padding','compact');

nexttile;
imagesc(lonSub, latSub, eddyScore);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colorbar;
caxis([0 3]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
quiver(Lon(1:vectorSkip:end,1:vectorSkip:end), Lat(1:vectorSkip:end,1:vectorSkip:end), ...
       ug_s(1:vectorSkip:end,1:vectorSkip:end), vg_s(1:vectorSkip:end,1:vectorSkip:end), ...
       1.5, 'k');
title('Eddy-support score: local SLA + |zeta/f| + OWrot');
xlabel('Longitude'); ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, double(candidateSupport));
set(gca,'YDir','normal');
axis tight; grid on; box on;
colorbar;
caxis([0 1]);
hold on;
contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
quiver(Lon(1:vectorSkip:end,1:vectorSkip:end), Lat(1:vectorSkip:end,1:vectorSkip:end), ...
       ug_s(1:vectorSkip:end,1:vectorSkip:end), vg_s(1:vectorSkip:end,1:vectorSkip:end), ...
       1.5, 'k');
title('Candidate support mask, not automatic eddy detection');
xlabel('Longitude'); ylabel('Latitude');

sgtitle('Guided eddy interpretation');

exportgraphics(gcf, fullfile(outdir, 'Fig02_candidate_support_mask.png'), 'Resolution', 180);

%% ===================== Figure 3: Guided questions =====================

figure('Color','w','Position',[120 120 1450 760]);

imagesc(lonSub, latSub, etaLocal);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colormap(gca, redblue_colormap(256));
colorbar;
caxis([-0.4 0.4]);
hold on;

contour(Lon, Lat, eta, -1:0.1:1, 'k', 'LineWidth', 0.4);
contour(Lon, Lat, OWrot, [0 0], 'm', 'LineWidth', 1.0);

quiver(Lon(1:vectorSkip:end,1:vectorSkip:end), Lat(1:vectorSkip:end,1:vectorSkip:end), ...
       ug_s(1:vectorSkip:end,1:vectorSkip:end), vg_s(1:vectorSkip:end,1:vectorSkip:end), ...
       1.5, 'k');

title('Choose eddy candidates by eye: local SLA, vectors, and OW_{rot}=0');
xlabel('Longitude');
ylabel('Latitude');

text(min(lonSub)+0.5, min(latSub)+1.0, { ...
    'Q1: 閉じた正偏差または負偏差はどこにあるか?', ...
    'Q2: その周囲のベクトルは回転しているか?', ...
    'Q3: \zeta/f はまとまった符号を持っているか?', ...
    'Q4: OW_{rot}>0 の領域と対応しているか?', ...
    '注意: OW_{rot}>0 だけでは渦とは断定しない。'}, ...
    'BackgroundColor','w', ...
    'EdgeColor',[0.7 0.7 0.7], ...
    'FontSize',11);

exportgraphics(gcf, fullfile(outdir, 'Fig03_guided_questions.png'), 'Resolution', 180);

%% ===================== Interpretation guide printed to console =====================

fprintf('\nInterpretation guide:\n');
fprintf('  Eddy-like structure should satisfy several conditions:\n');
fprintf('    1. closed or nearly closed SLA/ADT anomaly\n');
fprintf('    2. rotating geostrophic vectors\n');
fprintf('    3. coherent relative vorticity patch\n');
fprintf('    4. OWrot > 0, meaning rotation-dominant\n');
fprintf('\n');
fprintf('  Northern Hemisphere sign guide:\n');
fprintf('    positive SLA/ADT anomaly -> anticyclonic tendency -> clockwise rotation\n');
fprintf('    negative SLA/ADT anomaly -> cyclonic tendency -> counterclockwise rotation\n');
fprintf('\n');
fprintf('  OWrot is only a diagnostic, not an eddy detector.\n');

%% ===================== Save =====================

outmat = fullfile(outdir, 'SLA_eddy_interpretation.mat');

save(outmat, ...
    'lonSub', 'latSub', 'Lon', 'Lat', ...
    'eta', 'etaLocal', 'eta_bg', 'etaLabel', ...
    'ug_s', 'vg_s', 'speed', ...
    'zetaOverF', 'OWrot', 'KE', 'KElabel', ...
    'mask_localSLA', 'mask_zeta', 'mask_OWrot', ...
    'candidateSupport', 'eddyScore', ...
    'localSLA_thr', 'zetaF_thr', 'OWrot_thr', 'bgWindow', ...
    '-v7.3');

fprintf('\nSaved: %s\n', outmat);
disp('Done.');

%% ===================== Local functions =====================

function B = smooth2_nan(A, win)

    if win <= 1
        B = A;
        return
    end

    if mod(win,2) == 0
        win = win + 1;
    end

    K = ones(win, win);
    M = isfinite(A);

    A0 = A;
    A0(~M) = 0;

    num = conv2(A0, K, 'same');
    den = conv2(double(M), K, 'same');

    B = num ./ den;
    B(den == 0) = NaN;
end

function cmap = redblue_colormap(n)

    if nargin < 1
        n = 256;
    end

    n1 = floor(n/2);
    n2 = n - n1;

    blue = [linspace(0,1,n1)' linspace(0,1,n1)' ones(n1,1)];
    red  = [ones(n2,1) linspace(1,0,n2)' linspace(1,0,n2)'];

    cmap = [blue; red];
end

useSmoothedEta = true;
etaSmoothWindow = 5;

etaForContour = eta;

if useSmoothedEta
    etaForContour = smooth2_nan(eta, etaSmoothWindow);
end

contourInterval = 0.05;
contourLevelMin = -1.5;
contourLevelMax =  1.5;

etaMin = min(etaForContour(:), [], 'omitnan');
etaMax = max(etaForContour(:), [], 'omitnan');

levMin = max(contourLevelMin, floor(etaMin / contourInterval) * contourInterval);
levMax = min(contourLevelMax, ceil(etaMax / contourInterval) * contourInterval);

levels = levMin:contourInterval:levMax;

C = contourc(lonSub, latSub, etaForContour, levels);

C(1,k) = contour level
C(2,k) = number of points
C(:, k+1 : k+npt) = contour coordinates

k = 1;
while k < size(C,2)

    contourLevel = C(1,k);
    npt = C(2,k);

    x = C(1, k+1:k+npt);
    y = C(2, k+1:k+npt);

    k = k + npt + 1;
end

closeDist_deg = hypot(x(1) - x(end), y(1) - y(end));

if closeDist_deg > 1e-8
    continue;
end

in = inpolygon(Lon, Lat, x, y);
nInside = sum(in(:));

if nInside < minInsidePoints
    continue;
end

etaInside = etaForContour(in);

etaMaxInside = max(etaInside, [], 'omitnan');
etaMinInside = min(etaInside, [], 'omitnan');

ampPos = etaMaxInside - contourLevel;
ampNeg = contourLevel - etaMinInside;

if amplitude < minAmplitude
    continue;
end

% Area-weighted centroid
wArea = cosd(Lat(in));
lonCentroid = sum(Lon(in) .* wArea, 'omitnan') / sum(wArea, 'omitnan');
latCentroid = sum(Lat(in) .* wArea, 'omitnan') / sum(wArea, 'omitnan');

area_km2 = approximate_area_km2(in, lonSub, latSub, Lon, Lat);
radius_km = sqrt(area_km2 / pi);

if radius_km < minRadius_km || radius_km > maxRadius_km
    continue;
end

zetaFInside = zetaOverF(in);
OWrotInside = OWrot(in);

meanZetaF = mean(zetaFInside, 'omitnan');
owrotFraction = sum(OWrotInside > 0 & isfinite(OWrotInside)) ./ ...
                sum(isfinite(OWrotInside));

duplicateCenterDist_km = 50;

if dkm <= duplicateCenterDist_km
    group(end+1) = j;
end

duplicateSelection = "largest_area";

%% run05d_closed_contour_eddy_detection_student.m
% ============================================================
% Ocean Remote Sensing Practical
% Step 4: Simple closed-contour eddy detection from SLA/ADT
%
% Input:
%   run05b_diagnostics_output/OW_EKE_vorticity.mat
%
% Purpose:
%   - Extract closed SLA/ADT contours
%   - Treat each closed contour as an eddy-region candidate
%   - Calculate:
%       1. contour polygon
%       2. area-weighted centroid
%       3. SLA/ADT maximum or minimum inside the contour
%       4. amplitude relative to contour level
%       5. approximate radius
%       6. mean zeta/f
%       7. fraction of OWrot > 0
%
% Important:
%   This is a simplified educational version.
%   It is NOT the full AVISO / py-eddy-tracker algorithm.
%
% Concept:
%   Anticyclonic candidate:
%       closed contour surrounding a local SLA/ADT maximum
%
%   Cyclonic candidate:
%       closed contour surrounding a local SLA/ADT minimum
%
%   Eddy centroid:
%       area-weighted geometric centroid of grid cells inside contour
%
%   Eddy center:
%       SLA/ADT maximum or minimum inside contour
%
% ============================================================

clear; close all; clc;

%% ============================================================
% 0. Settings
% ============================================================

inmat = fullfile('run05b_diagnostics_output', 'OW_EKE_vorticity.mat');

outdir = 'run05d_closed_contour_eddy_detection_output';
if ~exist(outdir, 'dir')
    mkdir(outdir);
end

% Use smoothed or unsmoothed SLA/ADT for contour detection?
% For noisy fields, smoothing is useful.
useSmoothedEta = true;

% Smoothing window for SLA/ADT before contour extraction.
% Must be odd. 1 means no smoothing.
etaSmoothWindow = 5;

% Contour interval [m].
% Smaller interval gives more contours, but also more duplicate candidates.
contourInterval = 0.05;

% Minimum and maximum contour levels [m].
% These will be clipped by actual data range.
contourLevelMin = -1.5;
contourLevelMax =  1.5;

% Minimum number of grid cells inside a closed contour.
minInsidePoints = 20;

% Minimum eddy amplitude [m].
% amplitude = |SLA extreme - contour level|
minAmplitude = 0.05;

% Approximate eddy radius limits [km].
minRadius_km = 20;
maxRadius_km = 300;

% OWrot support condition.
% This is not mandatory for extraction, but used for diagnostic filtering.
minOWrotFraction = 0.10;

% zeta/f support condition.
% This is not mandatory for extraction, but used for diagnostic filtering.
minAbsMeanZetaF = 0.01;

% Duplicate merging:
% If two candidates have extreme centers closer than this distance,
% they are treated as the same eddy candidate.
duplicateCenterDist_km = 50;

% For each duplicate group, choose:
%   "largest_area" or "largest_amplitude"
duplicateSelection = "largest_area";

% Plot settings
vectorSkip = 4;

fprintf('=== run05d: Closed-contour eddy detection ===\n');
fprintf('Input file: %s\n', inmat);

%% ============================================================
% 1. Load diagnostics from run05b
% ============================================================

if ~isfile(inmat)
    error('Input MAT file not found: %s', inmat);
end

S = load(inmat);

lonSub = S.lonSub;
latSub = S.latSub;
Lon = S.Lon;
Lat = S.Lat;

eta = S.eta;
etaLabel = S.etaLabel;

ug_s = S.ug_s;
vg_s = S.vg_s;
speed = S.speed;

zetaOverF = S.zetaOverF;
OWrot = S.OWrot;
KE = S.KE;

fprintf('Loaded: %s\n', inmat);
fprintf('etaLabel = %s\n', etaLabel);
fprintf('eta size = [%d %d] = [lat lon]\n', size(eta,1), size(eta,2));

[nlat, nlon] = size(eta);

%% ============================================================
% 2. Prepare SLA/ADT field for contour extraction
% ============================================================

etaForContour = eta;

if useSmoothedEta
    etaForContour = smooth2_nan(eta, etaSmoothWindow);
end

% Keep land/invalid mask
etaForContour(~isfinite(eta)) = NaN;

etaMin = min(etaForContour(:), [], 'omitnan');
etaMax = max(etaForContour(:), [], 'omitnan');

levMin = max(contourLevelMin, floor(etaMin / contourInterval) * contourInterval);
levMax = min(contourLevelMax, ceil(etaMax / contourInterval) * contourInterval);

levels = levMin:contourInterval:levMax;

fprintf('\nContour settings:\n');
fprintf('  eta range for contour = %.3f to %.3f m\n', etaMin, etaMax);
fprintf('  contour levels = %.3f to %.3f m, interval %.3f m\n', ...
    levMin, levMax, contourInterval);
fprintf('  number of levels = %d\n', numel(levels));

%% ============================================================
% 3. Extract all closed contours
% ============================================================

% contourc requires x vector, y vector, and Z matrix of size [length(y), length(x)].
C = contourc(lonSub, latSub, etaForContour, levels);

allCandidates = struct([]);

icand = 0;
k = 1;

fprintf('\nExtracting closed contours...\n');

while k < size(C,2)

    contourLevel = C(1,k);
    npt = C(2,k);

    x = C(1, k+1:k+npt);
    y = C(2, k+1:k+npt);

    k = k + npt + 1;

    % Need enough points to make a polygon
    if npt < 4
        continue;
    end

    % Check if contour is closed
    closeDist_deg = hypot(x(1) - x(end), y(1) - y(end));

    % Tolerance is small because contourc usually repeats the first point.
    if closeDist_deg > 1e-8
        continue;
    end

    % Skip contours touching domain boundary.
    % These are often not true closed eddy contours within the domain.
    touchesBoundary = any(abs(x - min(lonSub)) < 1e-8) || ...
                      any(abs(x - max(lonSub)) < 1e-8) || ...
                      any(abs(y - min(latSub)) < 1e-8) || ...
                      any(abs(y - max(latSub)) < 1e-8);

    if touchesBoundary
        continue;
    end

    % Points inside polygon
    in = inpolygon(Lon, Lat, x, y);

    nInside = sum(in(:));

    if nInside < minInsidePoints
        continue;
    end

    etaInside = etaForContour(in);

    if all(~isfinite(etaInside))
        continue;
    end

    % Extreme values inside contour
    etaMaxInside = max(etaInside, [], 'omitnan');
    etaMinInside = min(etaInside, [], 'omitnan');

    ampPos = etaMaxInside - contourLevel;
    ampNeg = contourLevel - etaMinInside;

    % Decide polarity from which side has larger amplitude.
    % polarity = +1 : positive SLA/ADT eddy candidate
    % polarity = -1 : negative SLA/ADT eddy candidate
    if ampPos >= ampNeg
        polarity = +1;
        amplitude = ampPos;

        tmp = etaForContour;
        tmp(~in) = NaN;
        [etaExtreme, indExtreme] = max(tmp(:), [], 'omitnan');

    else
        polarity = -1;
        amplitude = ampNeg;

        tmp = etaForContour;
        tmp(~in) = NaN;
        [etaExtreme, indExtreme] = min(tmp(:), [], 'omitnan');
    end

    if amplitude < minAmplitude
        continue;
    end

    [iyExtreme, ixExtreme] = ind2sub(size(etaForContour), indExtreme);

    lonExtreme = Lon(iyExtreme, ixExtreme);
    latExtreme = Lat(iyExtreme, ixExtreme);

    % Area-weighted centroid.
    % On lon-lat grid, approximate grid area weighting by cos(latitude).
    wArea = cosd(Lat(in));
    lonCentroid = sum(Lon(in) .* wArea, 'omitnan') / sum(wArea, 'omitnan');
    latCentroid = sum(Lat(in) .* wArea, 'omitnan') / sum(wArea, 'omitnan');

    % Approximate area and radius.
    % Estimate grid-cell area from local grid spacing.
    area_km2 = approximate_area_km2(in, lonSub, latSub, Lon, Lat);
    radius_km = sqrt(area_km2 / pi);

    if radius_km < minRadius_km || radius_km > maxRadius_km
        continue;
    end

    % Diagnostics inside contour
    zetaFInside = zetaOverF(in);
    OWrotInside = OWrot(in);
    speedInside = speed(in);
    KEInside = KE(in);

    meanZetaF = mean(zetaFInside, 'omitnan');
    medianZetaF = median(zetaFInside, 'omitnan');
    meanSpeed = mean(speedInside, 'omitnan');
    maxSpeed = max(speedInside, [], 'omitnan');
    meanKE = mean(KEInside, 'omitnan');

    owrotFraction = sum(OWrotInside > 0 & isfinite(OWrotInside)) ./ ...
                    sum(isfinite(OWrotInside));

    if isempty(owrotFraction) || isnan(owrotFraction)
        owrotFraction = 0;
    end

    % Dynamical support flag.
    % This is diagnostic only.
    supportFlag = (owrotFraction >= minOWrotFraction) && ...
                  (abs(meanZetaF) >= minAbsMeanZetaF);

    % Expected zeta/f sign:
    % In both hemispheres:
    %   anticyclonic tendency: zeta/f < 0
    %   cyclonic tendency    : zeta/f > 0
    %
    % But polarity from SLA sign is not always enough near strong fronts.
    if polarity > 0
        polarityText = "positive-SLA/ADT candidate";
    else
        polarityText = "negative-SLA/ADT candidate";
    end

    icand = icand + 1;

    allCandidates(icand).id_raw = icand;
    allCandidates(icand).level = contourLevel;
    allCandidates(icand).x = x;
    allCandidates(icand).y = y;
    allCandidates(icand).nContourPoints = npt;
    allCandidates(icand).nInside = nInside;

    allCandidates(icand).polarity = polarity;
    allCandidates(icand).polarityText = polarityText;

    allCandidates(icand).amplitude = amplitude;
    allCandidates(icand).etaExtreme = etaExtreme;
    allCandidates(icand).lonExtreme = lonExtreme;
    allCandidates(icand).latExtreme = latExtreme;
    allCandidates(icand).ixExtreme = ixExtreme;
    allCandidates(icand).iyExtreme = iyExtreme;

    allCandidates(icand).lonCentroid = lonCentroid;
    allCandidates(icand).latCentroid = latCentroid;

    allCandidates(icand).area_km2 = area_km2;
    allCandidates(icand).radius_km = radius_km;

    allCandidates(icand).meanZetaF = meanZetaF;
    allCandidates(icand).medianZetaF = medianZetaF;
    allCandidates(icand).owrotFraction = owrotFraction;
    allCandidates(icand).meanSpeed = meanSpeed;
    allCandidates(icand).maxSpeed = maxSpeed;
    allCandidates(icand).meanKE = meanKE;
    allCandidates(icand).supportFlag = supportFlag;
end

fprintf('  number of raw closed-contour candidates = %d\n', numel(allCandidates));

if isempty(allCandidates)
    warning('No closed-contour eddy candidates found. Try lowering thresholds.');
end

%% ============================================================
% 4. Merge duplicate nested contours
% ============================================================

% Closed contours are often nested around the same SLA maximum/minimum.
% Here we merge candidates whose extreme centers are close to each other.

% Add selection-related fields to all raw candidates first.
% This makes all structures have exactly the same fields.
for q = 1:numel(allCandidates)
    allCandidates(q).id = NaN;
    allCandidates(q).rawGroup = [];
    allCandidates(q).nMerged = 1;
end

% IMPORTANT:
% Initialize selectedCandidates as an empty structure array
% with exactly the same fields as allCandidates.
selectedCandidates = allCandidates([]);

if ~isempty(allCandidates)
    nRaw = numel(allCandidates);
    used = false(nRaw,1);

    isel = 0;

    for i = 1:nRaw

        if used(i)
            continue;
        end

        lon0 = allCandidates(i).lonExtreme;
        lat0 = allCandidates(i).latExtreme;
        pol0 = allCandidates(i).polarity;

        group = i;

        for j = i+1:nRaw

            if used(j)
                continue;
            end

            if allCandidates(j).polarity ~= pol0
                continue;
            end

            dkm = haversine_km(lon0, lat0, ...
                               allCandidates(j).lonExtreme, ...
                               allCandidates(j).latExtreme);

            if dkm <= duplicateCenterDist_km
                group(end+1) = j; %#ok<SAGROW>
            end
        end

        used(group) = true;

        % Choose representative contour from the group
        if duplicateSelection == "largest_amplitude"

            amps = arrayfun(@(s) s.amplitude, allCandidates(group));
            [~, ibest] = max(amps);

        else
            areas = arrayfun(@(s) s.area_km2, allCandidates(group));
            [~, ibest] = max(areas);
        end

        bestIndex = group(ibest);

        isel = isel + 1;

        % IMPORTANT:
        % Add the extra fields BEFORE assigning to selectedCandidates.
        % Otherwise MATLAB complains about dissimilar structures after
        % the first selected candidate.
        tmpCandidate = allCandidates(bestIndex);
        tmpCandidate.id = isel;
        tmpCandidate.rawGroup = group;
        tmpCandidate.nMerged = numel(group);

        selectedCandidates(isel) = tmpCandidate; %#ok<SAGROW>
    end
end

fprintf('  number of selected candidates after merging = %d\n', numel(selectedCandidates));

%% ============================================================
% 5. Create candidate table
% ============================================================

if isempty(selectedCandidates)

    T = table();

else

    nSel = numel(selectedCandidates);

    id = zeros(nSel,1);
    polarity = zeros(nSel,1);
    level = zeros(nSel,1);
    amplitude = zeros(nSel,1);
    lonCentroid = zeros(nSel,1);
    latCentroid = zeros(nSel,1);
    lonExtreme = zeros(nSel,1);
    latExtreme = zeros(nSel,1);
    etaExtreme = zeros(nSel,1);
    area_km2 = zeros(nSel,1);
    radius_km = zeros(nSel,1);
    meanZetaF = zeros(nSel,1);
    owrotFraction = zeros(nSel,1);
    meanSpeed = zeros(nSel,1);
    maxSpeed = zeros(nSel,1);
    supportFlag = false(nSel,1);
    nMerged = zeros(nSel,1);

    for i = 1:nSel
        id(i) = selectedCandidates(i).id;
        polarity(i) = selectedCandidates(i).polarity;
        level(i) = selectedCandidates(i).level;
        amplitude(i) = selectedCandidates(i).amplitude;
        lonCentroid(i) = selectedCandidates(i).lonCentroid;
        latCentroid(i) = selectedCandidates(i).latCentroid;
        lonExtreme(i) = selectedCandidates(i).lonExtreme;
        latExtreme(i) = selectedCandidates(i).latExtreme;
        etaExtreme(i) = selectedCandidates(i).etaExtreme;
        area_km2(i) = selectedCandidates(i).area_km2;
        radius_km(i) = selectedCandidates(i).radius_km;
        meanZetaF(i) = selectedCandidates(i).meanZetaF;
        owrotFraction(i) = selectedCandidates(i).owrotFraction;
        meanSpeed(i) = selectedCandidates(i).meanSpeed;
        maxSpeed(i) = selectedCandidates(i).maxSpeed;
        supportFlag(i) = selectedCandidates(i).supportFlag;
        nMerged(i) = selectedCandidates(i).nMerged;
    end

    T = table(id, polarity, level, amplitude, ...
        lonCentroid, latCentroid, lonExtreme, latExtreme, etaExtreme, ...
        area_km2, radius_km, meanZetaF, owrotFraction, ...
        meanSpeed, maxSpeed, supportFlag, nMerged);

    % Sort by amplitude descending
    T = sortrows(T, 'amplitude', 'descend');
end

csvFile = fullfile(outdir, 'closed_contour_eddy_candidates.csv');
writetable(T, csvFile);

fprintf('\nSaved candidate table:\n  %s\n', csvFile);

%% ============================================================
% 6. Figure: all raw closed contours
% ============================================================

figure('Color','w','Position',[80 80 1500 850]);

imagesc(lonSub, latSub, etaForContour);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colorbar;
caxis([-1 1]);
hold on;

contour(Lon, Lat, etaForContour, levels, 'Color', [0.4 0.4 0.4], 'LineWidth', 0.3);

for i = 1:numel(allCandidates)

    if allCandidates(i).polarity > 0
        col = [1.0 0.2 0.1];   % red-ish
    else
        col = [0.1 0.2 1.0];   % blue-ish
    end

    plot(allCandidates(i).x, allCandidates(i).y, '-', ...
        'Color', col, 'LineWidth', 0.8);
end

title('All raw closed SLA/ADT contour candidates');
xlabel('Longitude');
ylabel('Latitude');

exportgraphics(gcf, fullfile(outdir, 'Fig01_all_raw_closed_contours.png'), 'Resolution', 180);

%% ============================================================
% 7. Figure: selected eddy candidates
% ============================================================

figure('Color','w','Position',[80 80 1500 900]);

imagesc(lonSub, latSub, etaForContour);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colorbar;
caxis([-1 1]);
hold on;

contour(Lon, Lat, etaForContour, levels, 'Color', [0.55 0.55 0.55], 'LineWidth', 0.25);

% Velocity vectors
quiver(Lon(1:vectorSkip:end,1:vectorSkip:end), ...
       Lat(1:vectorSkip:end,1:vectorSkip:end), ...
       ug_s(1:vectorSkip:end,1:vectorSkip:end), ...
       vg_s(1:vectorSkip:end,1:vectorSkip:end), ...
       1.5, 'k');

for i = 1:numel(selectedCandidates)

    if selectedCandidates(i).polarity > 0
        col = [1.0 0.0 0.0];
        marker = '^';
    else
        col = [0.0 0.1 1.0];
        marker = 'v';
    end

    plot(selectedCandidates(i).x, selectedCandidates(i).y, '-', ...
        'Color', col, 'LineWidth', 1.8);

    % Centroid: circle
    plot(selectedCandidates(i).lonCentroid, selectedCandidates(i).latCentroid, ...
        'o', 'MarkerSize', 8, ...
        'MarkerFaceColor', 'w', ...
        'MarkerEdgeColor', col, ...
        'LineWidth', 1.5);

    % SLA/ADT extreme: triangle
    plot(selectedCandidates(i).lonExtreme, selectedCandidates(i).latExtreme, ...
        marker, 'MarkerSize', 8, ...
        'MarkerFaceColor', col, ...
        'MarkerEdgeColor', 'k', ...
        'LineWidth', 1.0);

    text(selectedCandidates(i).lonCentroid, selectedCandidates(i).latCentroid, ...
        sprintf(' %d', i), ...
        'Color', 'k', ...
        'FontSize', 9, ...
        'FontWeight', 'bold', ...
        'BackgroundColor', 'w', ...
        'Margin', 1);
end

title({'Selected closed-contour eddy candidates', ...
       'Contour = eddy boundary candidate; circle = centroid; triangle = SLA/ADT extreme'});
xlabel('Longitude');
ylabel('Latitude');

exportgraphics(gcf, fullfile(outdir, 'Fig02_selected_eddy_candidates.png'), 'Resolution', 180);

%% ============================================================
% 8. Figure: selected candidates over OWrot and zeta/f
% ============================================================

figure('Color','w','Position',[80 80 1500 750]);
tiledlayout(1,2,'TileSpacing','compact','Padding','compact');

nexttile;
imagesc(lonSub, latSub, zetaOverF);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colormap(gca, redblue_colormap(256));
colorbar;
caxis([-0.5 0.5]);
hold on;
contour(Lon, Lat, etaForContour, levels, 'Color', [0.35 0.35 0.35], 'LineWidth', 0.25);
plot_selected_candidates(selectedCandidates);
title('\zeta/f with selected eddy contours');
xlabel('Longitude');
ylabel('Latitude');

nexttile;
imagesc(lonSub, latSub, OWrot * 1e10);
set(gca,'YDir','normal');
axis tight; grid on; box on;
colormap(gca, redblue_colormap(256));
colorbar;
caxis([-5 5]);
hold on;
contour(Lon, Lat, OWrot, [0 0], 'm', 'LineWidth', 0.8);
contour(Lon, Lat, etaForContour, levels, 'Color', [0.35 0.35 0.35], 'LineWidth', 0.25);
plot_selected_candidates(selectedCandidates);
title('OW_{rot} with selected eddy contours');
xlabel('Longitude');
ylabel('Latitude');

sgtitle('Dynamical consistency check for selected closed-contour candidates');

exportgraphics(gcf, fullfile(outdir, 'Fig03_candidates_on_zetaF_OWrot.png'), 'Resolution', 180);

%% ============================================================
% 9. Print interpretation guide
% ============================================================

fprintf('\nInterpretation guide:\n');
fprintf('  Red contours  : positive SLA/ADT closed-contour candidates\n');
fprintf('  Blue contours : negative SLA/ADT closed-contour candidates\n');
fprintf('  Circle        : area-weighted centroid of the closed contour\n');
fprintf('  Triangle      : SLA/ADT maximum or minimum inside the contour\n');
fprintf('\n');
fprintf('Important:\n');
fprintf('  This is not a full AVISO eddy tracker.\n');
fprintf('  It is a simplified closed-contour detector for education.\n');
fprintf('  Check SLA shape, vectors, zeta/f, and OWrot before calling it an eddy.\n');

if ~isempty(T)
    fprintf('\nTop 10 candidates by amplitude:\n');
    disp(T(1:min(10,height(T)), :));
end

%% ============================================================
% 10. Save MAT file
% ============================================================

outmat = fullfile(outdir, 'closed_contour_eddy_detection.mat');

save(outmat, ...
    'lonSub', 'latSub', 'Lon', 'Lat', ...
    'eta', 'etaForContour', 'etaLabel', ...
    'ug_s', 'vg_s', 'speed', ...
    'zetaOverF', 'OWrot', 'KE', ...
    'allCandidates', 'selectedCandidates', 'T', ...
    'levels', 'contourInterval', ...
    'minInsidePoints', 'minAmplitude', ...
    'minRadius_km', 'maxRadius_km', ...
    'minOWrotFraction', 'minAbsMeanZetaF', ...
    'duplicateCenterDist_km', 'duplicateSelection', ...
    'useSmoothedEta', 'etaSmoothWindow', ...
    '-v7.3');

fprintf('\nSaved MAT file:\n  %s\n', outmat);
disp('Done.');

%% ============================================================
% Local functions
% ============================================================

function B = smooth2_nan(A, win)

    if win <= 1
        B = A;
        return
    end

    if mod(win,2) == 0
        win = win + 1;
    end

    K = ones(win, win);
    M = isfinite(A);

    A0 = A;
    A0(~M) = 0;

    num = conv2(A0, K, 'same');
    den = conv2(double(M), K, 'same');

    B = num ./ den;
    B(den == 0) = NaN;
end

function area_km2 = approximate_area_km2(in, lon, lat, Lon, Lat)

    Re = 6371000;  % m

    lon = double(lon(:));
    lat = double(lat(:));

    if numel(lon) < 2 || numel(lat) < 2
        area_km2 = NaN;
        return
    end

    dlon_rad = deg2rad(abs(median(diff(lon), 'omitnan')));
    dlat_rad = deg2rad(abs(median(diff(lat), 'omitnan')));

    % grid cell area varies with latitude
    cellArea_m2 = (Re^2) .* cosd(Lat) .* dlon_rad .* dlat_rad;
    cellArea_m2(~in) = NaN;

    area_km2 = sum(cellArea_m2(:), 'omitnan') / 1e6;
end

function dkm = haversine_km(lon1, lat1, lon2, lat2)

    R = 6371.0; % km

    lon1 = deg2rad(lon1);
    lat1 = deg2rad(lat1);
    lon2 = deg2rad(lon2);
    lat2 = deg2rad(lat2);

    dlon = lon2 - lon1;
    dlat = lat2 - lat1;

    a = sin(dlat/2).^2 + cos(lat1).*cos(lat2).*sin(dlon/2).^2;
    c = 2 * atan2(sqrt(a), sqrt(1-a));

    dkm = R * c;
end

function cmap = redblue_colormap(n)

    if nargin < 1
        n = 256;
    end

    n1 = floor(n/2);
    n2 = n - n1;

    blue = [linspace(0,1,n1)' linspace(0,1,n1)' ones(n1,1)];
    red  = [ones(n2,1) linspace(1,0,n2)' linspace(1,0,n2)'];

    cmap = [blue; red];
end

function plot_selected_candidates(selectedCandidates)

    for ii = 1:numel(selectedCandidates)

        if selectedCandidates(ii).polarity > 0
            col = [1.0 0.0 0.0];
            marker = '^';
        else
            col = [0.0 0.1 1.0];
            marker = 'v';
        end

        plot(selectedCandidates(ii).x, selectedCandidates(ii).y, '-', ...
            'Color', col, 'LineWidth', 1.5);

        plot(selectedCandidates(ii).lonCentroid, selectedCandidates(ii).latCentroid, ...
            'o', 'MarkerSize', 7, ...
            'MarkerFaceColor', 'w', ...
            'MarkerEdgeColor', col, ...
            'LineWidth', 1.2);

        plot(selectedCandidates(ii).lonExtreme, selectedCandidates(ii).latExtreme, ...
            marker, 'MarkerSize', 7, ...
            'MarkerFaceColor', col, ...
            'MarkerEdgeColor', 'k', ...
            'LineWidth', 0.8);
    end
end