sims = ["swop_20230304_27012_09_SD_AM_xC", "swop_20230304_27012_09_SD_UM_xC"];

direc_root = getenv('GEMINI_SIM_ROOT');
xg = aurogem.grid.read(fullfile(direc_root, sims(1)));
Bmag = median(xg.Bmag(:));
x = xg.x2(3:end-2);
y = xg.x3(3:end-2);
z = xg.x1(3:end-2);
dz = xg.dx1h;
[Xh, Yh] = ndgrid(x, y);
[Zv, Yv] = ndgrid(z, y);
[~, dX] = gradient(Xh);
[dY, ~] = gradient(Yh);

%%
for s = 'AB'
    disp(s)
    sim = sims(s - 64);
    direc = fullfile(direc_root, sim);
    cfg = gemini3d.read.config(direc);
    dat_h5 = fullfile(direc, gemini3d.datelab(cfg.times(end)) + ".h5");
    dat.(s).v2 = h5read(dat_h5, '/v2avgall');
    dat.(s).v3 = h5read(dat_h5, '/v3avgall');
    dat.(s).j2 = h5read(dat_h5, '/J2all');
    dat.(s).j3 = h5read(dat_h5, '/J3all');
    dat.(s).sigP = h5read(dat_h5, '/sigPall');
    dat.(s).sigH = abs(h5read(dat_h5, '/sigHall'));

    v2 = squeeze(dat.(s).v2(end, :, :));
    v3 = squeeze(dat.(s).v3(end, :, :));
    dat.(s).E2 = -permute(repmat(v3, [1, 1, length(z)]), [3, 1, 2]) * Bmag;
    dat.(s).E3 = permute(repmat(v2, [1, 1, length(z)]), [3, 1, 2]) * Bmag;
    
    Emag = sqrt(dat.(s).E2.^2 + dat.(s).E3.^2);
    dat.(s).jP = (dat.(s).j2 .* dat.(s).E2 + dat.(s).j3 .* dat.(s).E3) ./ Emag;
    dat.(s).jH = (dat.(s).j2 .* dat.(s).E3 - dat.(s).j3 .* dat.(s).E2) ./ Emag;

    dat.(s).SigP = squeeze(sum(dat.(s).sigP .* dz, 1));
    dat.(s).SigH = squeeze(sum(dat.(s).sigH .* dz, 1));

    [~, dE22] = gradient(squeeze(dat.(s).E2(1, :, :)));
    [dE33, ~] = gradient(squeeze(dat.(s).E3(1, :, :)));
    [dSIGH3, dSIGH2] = gradient(dat.(s).SigH);
    dat.(s).jA = dat.(s).SigP .* (dE22 ./ dX + dE33 ./ dY);
    dat.(s).jC = (dSIGH2 .* squeeze(dat.(s).E3(1, :, :)) ./ dX - dSIGH3 .* squeeze(dat.(s).E2(1, :, :)) ./ dY);
end

%%
lim.x = [-1, 1] * 120;
lim.y = [-1, 1] * 60;
lim.z = [80, 165];
lim.jP = [-1, 1] * 45;
lim.jH = [-1, 1] * 95;
lim.E = [-1, 1] * 48;
lim.S = [0, 85];
lim.jC = [-1, 1] * 80;

close all
figure('PaperUnits', 'inches', 'PaperPosition', [0, 0, 6.5, 7] * 4)
tiledlayout(6, 2, 'TileSpacing', 'compact', 'Padding', 'tight')

set(0, 'defaultSurfaceEdgeColor', 'flat')
aurogem.tools.setall(0, 'FontName', 'Arial')
aurogem.tools.setall(0, 'FontSize', 10 * 4)
aurogem.tools.setall(0, 'Multiplier', 1)
colorcet = @aurogem.tools.colorcet;

nexttile
pcolor(Yv / 1e3, Zv / 1e3, squeeze(dat.A.jP(:, length(x)/2, :)) * 1e6)
xlim(lim.y); ylim(lim.z); clim(lim.jP)
colormap(gca, colorcet('D1A'))
title('Accelerated')
xlabel('Mag. north (km)'); ylabel('Mag. up (km)')

nexttile
pcolor(Yv / 1e3, Zv / 1e3, squeeze(dat.B.jP(:, length(x)/2, :)) * 1e6)
xlim(lim.y); ylim(lim.z); clim(lim.jP)
colormap(gca, colorcet('D1A'))
clb = colorbar; clb.Label.String = 'j_P (uA/m^2)';
title('Unaccelerated')
xlabel('Mag. north (km)'); ylabel('Mag. up (km)')

nexttile
pcolor(Yv / 1e3, Zv / 1e3, squeeze(dat.A.jH(:, length(x)/2, :)) * 1e6)
xlim(lim.y); ylim(lim.z); clim(lim.jH)
colormap(gca, colorcet('D1A'))
xlabel('Mag. north (km)'); ylabel('Mag. up (km)')

nexttile
pcolor(Yv / 1e3, Zv / 1e3, squeeze(dat.B.jH(:, length(x)/2, :)) * 1e6)
xlim(lim.y); ylim(lim.z); clim(lim.jH)
colormap(gca, colorcet('D1A'))
clb = colorbar; clb.Label.String = 'j_H (uA/m^2)';
xlabel('Mag. north (km)'); ylabel('Mag. up (km)')

nexttile
pcolor(Xh / 1e3, Yh / 1e3, squeeze(dat.A.E3(end, :, :)) * 1e3)
xlim(lim.x); ylim(lim.y); clim(lim.E)
colormap(gca, colorcet('D13'))
xlabel('Mag. east (km)'); ylabel('Mag. north (km)')

nexttile
pcolor(Xh / 1e3, Yh / 1e3, squeeze(dat.B.E3(end, :, :)) * 1e3)
xlim(lim.x); ylim(lim.y); clim(lim.E)
colormap(gca, colorcet('D13'))
clb = colorbar; clb.Label.String = 'E_{north} (mV/m)';
xlabel('Mag. east (km)'); ylabel('Mag. north (km)')

nexttile
pcolor(Xh / 1e3, Yh / 1e3, dat.A.SigH)
xlim(lim.x); ylim(lim.y); clim(lim.S)
colormap(gca, colorcet('L18'))
xlabel('Mag. east (km)'); ylabel('Mag. north (km)')

nexttile
pcolor(Xh / 1e3, Yh / 1e3, dat.B.SigH)
xlim(lim.x); ylim(lim.y); clim(lim.S)
colormap(gca, colorcet('L18'))
clb = colorbar; clb.Label.String = '\Sigma_H (S)';
xlabel('Mag. east (km)'); ylabel('Mag. north (km)')

nexttile
pcolor(Xh / 1e3, Yh / 1e3, dat.A.jC * 1e6)
xlim(lim.x); ylim(lim.y); clim(lim.jC)
colormap(gca, colorcet('D1A'))
xlabel('Mag. east (km)'); ylabel('Mag. north (km)')

nexttile
pcolor(Xh / 1e3, Yh / 1e3, dat.B.jC * 1e6)
xlim(lim.x); ylim(lim.y); clim(lim.jC)
colormap(gca, colorcet('D1A'))
clb = colorbar; clb.Label.String = '(\bfb\rm\times\bfE\rm)\cdot\nabla\Sigma_H (uA/m^2)';
xlabel('Mag. east (km)'); ylabel('Mag. north (km)')

nexttile
pcolor(Xh / 1e3, Yh / 1e3, dat.A.jA * 1e6)
xlim(lim.x); ylim(lim.y); clim(lim.jC)
colormap(gca, colorcet('D1A'))
xlabel('Mag. east (km)'); ylabel('Mag. north (km)')

nexttile
pcolor(Xh / 1e3, Yh / 1e3, dat.B.jA * 1e6)
xlim(lim.x); ylim(lim.y); clim(lim.jC)
colormap(gca, colorcet('D1A'))
clb = colorbar; clb.Label.String = '(\Sigma_P\nabla\cdot\bfE\rm (uA/m^2)';
xlabel('Mag. east (km)'); ylabel('Mag. north (km)')

filename = 'plots/00_joule_heating_viii.png';
print(gcf, filename, '-dpng', '-r96');
close all

%%
ll = 0;
tt = -10;
ww = 1110;
hh = 439;
im = imread(filename);
for i = 1:12
    px = ll + ww * mod(i - 1, 2);
    py = tt + hh * floor((i - 1) / 2);
    im = insertText(im, [px, py], char(96 + i), ...
        'AnchorPoint', 'LeftTop', 'BoxOpacity', 0, 'FontSize', 12 * 4, 'FontColor', 'k');
end
imwrite(im, filename, 'png')