comparisons = [...
["swop_20230210_35487_09_SD_AM_AC"; "swop_20230210_35487_09_PF_AM_AC"], ... # I
["swop_20230212_37331_09_SD_AM_xC"; "swop_20230212_37331_09_PF_AM_xC"], ... # II
["swop_20230304_27012_09_SD_AM_xC"; "swop_20230304_27012_09_NB_AM_xC"], ... # III
["swop_20230314_24547_09_SD_AM_AC"; "swop_20230314_24547_09_PF_AM_AC"], ... # IV
["swop_20230319_30210_09_SD_AM_xB"; "swop_20230319_30210_09_PF_AM_xB"], ... # V
...
["swop_20230210_35487_09_SD_AM_AC"; "swop_20230210_35487_09_SD_UM_AC"], ... # VI
["swop_20230212_37331_09_SD_AM_xC"; "swop_20230212_37331_09_SD_UM_xC"], ... # VII
["swop_20230304_27012_09_SD_AM_xC"; "swop_20230304_27012_09_SD_UM_xC"], ... # VIII
["swop_20230304_36829_09_SD_AM_xB"; "swop_20230304_36829_09_SD_UM_xB"], ... # IX
...
["swop_20230210_35487_09_SD_AM_AC"; "swop_20230210_35487_09_SD_AM_xA"], ... # X
["swop_20230314_24547_09_SD_AM_AC"; "swop_20230314_24547_09_SD_AM_xA"], ... # XI
];

resf = 4;
linw = 2;
fnts = 10 * resf;
pprw = 7 * resf;
pprh = 1.35 * resf;
fntn = 'Arial';
ar = [3, 1, 1];
letter_pos = [[0, 457]; [310, 457]] * resf;
direc_root = getenv('GEMINI_SIM_ROOT');

for comp = comparisons(:, 8)
for s = 'AB'
sim.(s) = char(comp(s-64));
direc = fullfile(direc_root, sim.(s));
cfg = gemini3d.read.config(direc);
if s == 'A'
    xg = aurogem.grid.read(direc);
    x = xg.x2(3:end-2) / 1e3;
    y = xg.x3(3:end-2) / 1e3;
    z = xg.x1(3:end-2) / 1e3;
    [X, Y] = ndgrid(x, y);
    dz = xg.dx1h / 1e3;
end
dat_file = fullfile(direc, gemini3d.datelab(cfg.times(end)) + ".h5");
dat.v2 = h5read(dat_file, '/v2avgall');
dat.v3 = h5read(dat_file, '/v3avgall');
dat.J2 = h5read(dat_file, '/J2all');
dat.J3 = h5read(dat_file, '/J3all');

[base, ~] = aurogem.tools.flux_tube_pars(direc);

v2 = squeeze(dat.v2(end, :, :));
v3 = squeeze(dat.v3(end, :, :));
E2 = -permute(repmat(v3, [1, 1, length(z)]), [3, 1, 2]) .* xg.Bmag;
E3 = permute(repmat(v2, [1, 1, length(z)]), [3, 1, 2]) .* xg.Bmag;
j2 = dat.J2; % A/m^2
j3 = dat.J3;

joule_3d = j2.*E2 + j3.*E3; % W/m^3
joule.(s) = squeeze(sum(joule_3d.*dz)) * 1e6; % mW/m^2

lim.u.(s) = [quantile(joule.(s)(:), 0.01), quantile(joule.(s)(:), 0.99)];
fprintf('%s quantiles: %.2f, %.2f\n', s, lim.u.(s)(1), lim.u.(s)(2))
lim.u.(s) = round(lim.u.(s) * 5) / 5 + [-1, 1] * 0.05;

ctrs.(s) = round(quantile(joule.(s)(:), [0.5, 0.95]), 1);

lim.x = base.limx * 0.98;
lim.y = base.limy * 0.98;
end

fprintf('%s v %s\n', sim.A, sim.B)

% hold on
% plot(joule.A(144, :))
% plot(joule.B(144, :))
% input('?')

close all
reset(0)
set(0, 'defaultSurfaceEdgeColor', 'flat')
set(0, 'defaultLineLineWidth', linw)
aurogem.tools.setall(0, 'FontName', fntn)
aurogem.tools.setall(0, 'FontSize', fnts)
aurogem.tools.setall(0, 'Multiplier', 1)
colorcet = @aurogem.tools.colorcet;

figure
set(gcf, 'PaperUnits', 'inches', 'PaperPosition', [0, 0, pprw, pprh], ...
    'Position', [10, 30, 900, 900*pprh/pprw])
tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'tight');

nexttile
hold on
pcolor(X, Y, joule.A)
contour(X, Y, joule.A, ctrs.A, 'ShowText', true, 'EdgeColor', 'k')
xlim(lim.x); ylim(lim.y); clim(lim.u.A)
xlabel('Mag. east (km)')
ylabel('Mag. north (km)')
colormap(gca, colorcet('L19'))
clb = colorbar(gca);
clb.Label.String = '\int_{||} \bfj\cdotE\rm (mW/m^2)';
pbaspect(ar)
set([gca, clb], 'FontSize', fnts)

nexttile
hold on
pcolor(X, Y, joule.B)
contour(X, Y, joule.B, ctrs.B, 'ShowText', true, 'EdgeColor', 'k')
colormap(gca, colorcet('L19'))
xlim(lim.x); ylim(lim.y); clim(lim.u.B)
xlabel('Mag. east (km)')
clb = colorbar(gca);
clb.Label.String = '\int_{||} \bfj\cdotE\rm (mW/m^2)';
pbaspect(ar)
set([gca, clb], 'FontSize', fnts)

filename = sprintf('plots/00_joule_heating_%s_v_%s.png', ...
    sim.A(6:end), sim.B(24:end));
print(gcf, filename, '-dpng', '-r96');
close all

im = imread(filename);
filename_compare = sprintf('plots/00_gemini_compare_%s_v_%s.png', ...
    sim.A(6:end), sim.B(24:end));
im_compare = imread(filename_compare);
im_new = uint8(ones(size(im, 1)+size(im_compare, 1), size(im_compare, 2), 3)) * 255;
im_new(1:size(im_compare, 1), :, :) = im_compare;
im_new(size(im_compare, 1)+1:end, 1:size(im, 2), :) = im;

for i = 1:length(letter_pos)
    im_new = insertText(im_new, letter_pos(i, :), char(96 + i + 4), ...
        'AnchorPoint', 'LeftTop', 'BoxOpacity', 0, 'FontSize', 14 * resf);
end

imwrite(im_new, filename, 'png')
end