function PLOT_ICAL3(l3)
%PLOT_ICAL Plot Density and velocity of RPC ICA RAW data.
%
% Example:
%    l2 = READICA; % Read a days worth of moment data
%    PLOT_ICA(l2)
%
% See also: IMPORT_L3, READICA_L3, ICAINITIALIZE

% Matlab script to plot ICA calibrated data
%
%
% Working with the data can be always done in the same way as in this
% sample, finding data through a logical expression or the find command. 
% To find a certain sector one simply
% adds a azim(1,:)==desired sector or (azim(1,:)<= desired sector &
% azim(2,:)>=desired sector).

%
% Normal mode data with onboard binning of mass channels according to lookup tables
% can be identified using the mass table number, which is
% non-zero for mass lookup table data.
% An important note regarding Normal mode data using mass lookup tables: For early data,
% before 2014-10-30, the energy table of ICA did not correspond to
% expectations, there was an offset in the high voltages. This made both
% elevation values and mass lookup tables wrong. The elevation tables has
% been corrected for this, but the onboard lookup tables where what they
% were. Therefore the columns of the tables do not fully agree with the
% intended ion species, this is mainly a problem at low to intermediate
% energies. The data in the archive has been expanded to 32 mass channels,
% but the onboard binning may cause some artificial broadening of the mass
% channel distribution.
%

if ~exist('energy.mat','file')
    disp('You need to run: ICAINITIALIZE before you run this script')
    return
end
% mass_lookup = load('mass_lookup.mat'); % not used here
energy = load('energy.mat');
E = energy.E;

figure(33)
clf
set(gcf,'renderer','painters');

% The spectrograms for different masses and sectors have the same time so
% find all unique time instances:
[time_instances, ia, ~] = unique(l3.iontime);
%ia is index variable to pick out only one value per time instance

pl=6; % Use only high post acceleration (6 for a value between 0 and 7)


%% An example to create a mass matrix plot ( energy vs mass channel number)
% Most useful for burst mode data when mtable = 0

mass_matrix=zeros(96,32); % sum over all times with same swversion, azimuths, assume full resolution
ix=(l3.swversion == l3.swversion(1) & l3.mtable == 0 & l3.pacclevel == pl);
if ~isempty(ix)
    for mch=0:31
        aa=l3.masschannel(ix)==mch;
        mass_matrix(:,mch+1)=sum(l3.flux(aa,:),1);
    end
    figure(33)
    pcolor(0:31,squeeze(E(l3.swversion(1),2,:)),log10(mass_matrix));
    shading flat;
    colorbar;
    set(gca,'yscale','log');

    % Similarly add all mass channels and azimuths and make an energy spectrogram
    spectrogram=zeros(96,length(time_instances));
    for jj=1:length(time_instances)
        ix=l3.iontime==time_instances(jj); %
        spectrogram(:,jj)=sum(l3.flux(ix,:),1);
    end
    figure(44)
    pcolor(time_instances,squeeze(E(l3.swversion(1),2,:)),log10(spectrogram));
    shading flat;
    colorbar;
    set(gca,'yscale','log');

    datetick('x')
else
    disp('No mass matrix data to show')
end

end
