function PLOT_ICAL4_CORR( l4 )
%PLOT_ICAL4_CORR Summary of this function goes here
%
% Example:
%    l4_corr = READICA_L4_CORR; % Read a days worth of moment data
%    PLOT_ICAL5(l5)
%
% See also: IMPORT_L4_CORR, READICA_L4_CORR

figure(1111);
set(gcf, 'Name','H+','NumberTitle','off');
pos=get(gcf, 'Position');
subplot(2,1,1)
[tsum, uniqt] = sum_along_time(l4.H, l4.time);
[tsum, uniqt] = insert_nan_at_time_gap(tsum, uniqt);
my_plot(tsum, uniqt, 0:95, [], 'Energy bin');
title(['H^{+} ' datestr(uniqt(1), 'yyyy-mm-dd')])
datetick x
axis tight

subplot(2,2,3)
my_plot(sum_along(l4.H, l4.elev), 0:15, 0:95, 'Elevation bin', 'Energy bin')
subplot(2,2,4)
my_plot(sum_along(l4.H, l4.azim), 0:15, 0:95, 'Azimuth bin', 'Energy bin')

figure(2222);
set(gcf, 'Name','He++','NumberTitle','off');
set(gcf, 'Position', pos+[20 -20 0 0]);
subplot(2,1,1)
[tsum, uniqt] = sum_along_time(l4.He2, l4.time);
[tsum, uniqt] = insert_nan_at_time_gap(tsum, uniqt);
my_plot(tsum, uniqt, 0:95, [], 'Energy bin');
title(['He^{2+} ' datestr(uniqt(1), 'yyyy-mm-dd')])
datetick x
axis tight

subplot(2,2,3)
my_plot(sum_along(l4.He2, l4.elev), 0:15, 0:95, 'Elevation bin', 'Energy bin')
subplot(2,2,4)
my_plot(sum_along(l4.He2, l4.azim), 0:15, 0:95, 'Azimuth bin', 'Energy bin')

figure(3333);
set(gcf, 'Name','He+','NumberTitle','off');
set(gcf, 'Position', pos+[40 -40 0 0]);
subplot(2,1,1)
[tsum, uniqt] = sum_along_time(l4.He, l4.time);
[tsum, uniqt] = insert_nan_at_time_gap(tsum, uniqt);
my_plot(tsum, uniqt, 0:95, [], 'Energy bin');
title(['He^{+} ' datestr(uniqt(1), 'yyyy-mm-dd')])
datetick x
axis tight

subplot(2,2,3)
my_plot(sum_along(l4.He, l4.elev), 0:15, 0:95, 'Elevation bin', 'Energy bin')
subplot(2,2,4)
my_plot(sum_along(l4.He, l4.azim), 0:15, 0:95, 'Azimuth bin', 'Energy bin')

figure(4444);
set(gcf, 'Name','Heavy','NumberTitle','off');
set(gcf, 'Position', pos+[60 -60 0 0]);
subplot(2,1,1)
[tsum, uniqt] = sum_along_time(l4.Heavy, l4.time);
[tsum, uniqt] = insert_nan_at_time_gap(tsum, uniqt);
my_plot(tsum, uniqt, 0:95, [], 'Energy bin');
title(['Heavy ' datestr(uniqt(1), 'yyyy-mm-dd')])
datetick x
axis tight

subplot(2,2,3)
my_plot(sum_along(l4.Heavy, l4.elev), 0:15, 0:95, 'Elevation bin', 'Energy bin')
subplot(2,2,4)
my_plot(sum_along(l4.Heavy, l4.azim), 0:15, 0:95, 'Azimuth bin', 'Energy bin')

    function [summed, uniqt] = sum_along_time(m, t)
        uniqt=unique(t);
        summed=zeros(96,length(uniqt));
        for i=1:length(uniqt)
           summed(:,i) = sum(m(t==uniqt(i),:), 'omitnan');
        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) 
        pcolor(xtic, ytic, log10(m))
        shading flat
        ylabel(ylbl)
        if ~isempty(xlbl); xlabel(xlbl); end
        colorbar
    end

    %%
    function [new_mx, new_t] = insert_nan_at_time_gap(mx, t)
        % 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)];
            i_add = i_add + 1;
        end
        new_mx = mx;
        new_t = t;
    end
end

