function PLOT_ICAL4_CORR_CTS( l4 )
%PLOT_ICAL4_CORR_CTS Summary of this function goes here
%
% Example:
%    l4_corr_cts = READICA_L4_CORR_CTS; % Read a days worth of moment data
%    PLOT_ICAL4_CORR_CTS(l4_corr_cts)
%
% See also: IMPORT_L4_CORR_CTS, READICA_L4_CORR_CTS

if ~exist('energy.mat','file')
    disp('You need to run: ICAINITIALIZE before you run this script')
    return
end

Energy = load('energy.mat');
E = squeeze(Energy.E(:,2,:));

subplot(2,1,1)
[tsum, time_instances, sw] = sum_along_time(l4);
[tsum, time_instances, sw] = insert_nan_at_time_gap(tsum, time_instances, sw);
uniqsw = unique(sw);
for s=1:length(uniqsw)
    sw1=uniqsw(s);
    my_plot(tsum(:,sw==sw1), time_instances(sw==sw1), log10(squeeze(E(sw1,:))), 'time', 'E (eV)');
    hold on
end
hold off
set(gca,'XLim', time_instances([1 end]))
title(datestr(time_instances(1), 'yyyy-mm-dd'))
datetick x
axis tight
set(gca,'YLim', log10([max(min(squeeze(E(uniqsw,:)),[],2)) max(max(squeeze(E(uniqsw,:)),[],2))]))
cx=caxis;
cx(1) = 0;
caxis(cx);

subplot(2,3,4)
elsum = sum_along(l4.counts, l4.elev);
my_plot(elsum, 0:15, 0:95, 'Elevation bin', 'Energy bin')

subplot(2,3,5)
msum = sum_along(l4.counts, l4.mass);
my_plot(msum, 0:31, 0:95, 'Mass bin', 'Energy bin')


subplot(2,3,6)
azsum = sum_along(l4.counts, l4.azim);
my_plot(azsum, 0:15, 0:95, 'Azimuth bin', 'Energy bin')

    function [summed, time_instances, sw] = sum_along_time(l4)
        time_instances=unique(l4.time);
        summed=zeros(96,length(time_instances));
        sw=zeros(length(time_instances),1);
        for i=1:length(time_instances)
           summed(:,i) = sum(l4.counts(l4.time==time_instances(i),:), 'omitnan');
           sw(i) = l4.swversion(find(l4.time==time_instances(i),1));
        end
        summed(summed<0) = 0;
    end

    function summed = sum_along(m, v)
        max_v=max(v);
        summed=zeros(96,max_v+1);
        for i=0:max_v
           summed(:,i+1) = sum(m(v==i,:), 'omitnan');
        end
        summed(summed<0) = 0;
    end

    function my_plot(m, xtic, ytic, xlbl, ylbl)
        sz = size(m);
        pcolor(xtic, ytic, log10(m))
        shading flat
        ylabel(ylbl)
        xlabel(xlbl)
        colorbar
    end
    %%
    function [new_mx, new_t, new_sw] = insert_nan_at_time_gap(mx, t, sw)
        % find gaps larger than standard 192s and insert a new timestamp, t+192s,
        % and a column of nan's in both sunions & heavyions
        seconds=192;
        gaps = find(diff(t) > seconds/86400);
        i_add = 0;
        for gap = 1:length(gaps)
            ix = gaps(gap) + i_add;

            t = [t(1:ix); t(ix) + (t(ix)-t(ix-1))/86400; t(ix+1:end)];
            mx = [mx(:,1:ix) nan(96, 1) mx(:,ix+1:end)];
            sw = [sw(1:ix); sw(ix); sw(ix+1:end)];
            i_add = i_add + 1;
        end
        new_mx = mx;
        new_t = t;
        new_sw = sw;
    end
end

