% Matlab script to plot ICA data
%
% For fully resolved mass channels the values
% will be between 0 and 31. The sample mass spectra shown in
% Nilsson et al., Birth of a comet magnetosphere: a spring of
% water ions, Science, 2015 is then the currently best guide to mass
% separation of data. Please contact the team if you have a need for
% mass separation. I future higher level product will be delivered to the
% archive providing ready made separation into the major ion species
% H+, He2+ and heavier ions of cometary origin ions (mainly water group ions).
%
% Working with the data is always done in the same way as in this
% sample, using 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).
% In this sample all elevation data is supposed to be binned into two
% bins only, starting on 0 or 8 (which was the case for the Earth flyby
% data). One can choose which bin with the
% variable el. To find the size of bins for a particular data set
% plot(elev(1,:)') and plot(azim(1,:)') to see the start elevation and
% azimuth indices. For highest resolution data all 16 sectors and
% elevations appears as start indices, whereas for lowest resolution it
% is only 0 and 8 (two bins).
% Similarly post-acceleration can be choosen as one of three intervals,
% defined in the code below.
% During Earth flyby (20050301) an alternating post acceleration
% mode was used, with 2 and 6 as reference levels.
% Tested with Matlab R2011a

% Working with mass lookup tables:
% At times the number of transmitted mass channel intervals is too low for
% just binning of adjacent mass channels to be meaningful. Then data is
% binned according to onboard lookup tables. The mass channels then
% actually reflects which column in the lookup table the binned data
% corresponds to. The lookup table in the CALIB directory contains information
% about which mass channels were binned together.
%
% Lookup table data 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 inermediate
% energies. The nominal order of physical mass (m/q amu) [1 32 16 4 2 8]
% For newer data, software version 9, the lookup tables do not
% correspond to assumed physical masses in the same way as the original
% table. They rather try to maximise the useful information.
%




if exist('iontime','var')
    
    figure(11)
    clf
    set(gcf,'renderer','painters');
    if exist('ICAINITIALIZE', 'file')
        ICAINITIALIZE
    elseif exist('icainitialize', 'file')
        icainitialize
    end
    % The spectrograms for different masses and sectors have the same time so
    % find all unique time instances:
    time_instances=iontime(find(diff(iontime)~=0));
    H = NaN * ones(96,length(time_instances)); % Nominally protons
    alpha = NaN * ones(96,length(time_instances)); % Nominally He++
    Heavy = NaN * ones(96,length(time_instances)); % Heavy ions, i.e. water group and heavier
    He = NaN * ones(96,length(time_instances)); % He+
    
    pl=6;
    el=0;
    
    %% Find software version
    if exist('GET_SW_VERSION', 'file')
        sw = GET_SW_VERSION(promsection(1), iontime(1));
    elseif exist('get_sw_version', 'file')
        sw = get_sw_version(promsection(1), iontime(1));
    end
    %% Find correct table
    pacc_levels = [1, 2, 2, 2, 2, 3, 3, 3];
    pacc = pacc_levels(pl+1);
    %% Get correct mass_lookup
    mass_lookup = squeeze(LookUp(pacc, sw, :, :));
    
    %% Useful when mtable is not 0, i.e for lookup table data
    % Note that we give the nominal ion species for the orignal tables
    % Care must be taken with this, both because of the problem with the table
    % and because there exist cross talk which must be properly treated
    for jj=1:length(time_instances)
        % filter on common causes
        ix = iontime==time_instances(jj) & pacclevel==pl & elev(1,:)==el & mtable~=0;
        if sum(ix) == 0
            continue        % optimize
        end
        aa=find(masschannel(1,:)==0 & ix);
        if ~isempty(aa)
            H(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
        end
        aa=find(masschannel(1,:)==4 & ix);
        if ~isempty(aa)
            alpha(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
        end
        aa=find((masschannel(1,:)==1 | masschannel(1,:)==2 | masschannel(1,:)==5) & ix);
        if ~isempty(aa)
            Heavy(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
        end
        aa=find(masschannel(1,:)==3 & ix);
        if ~isempty(aa)
            He(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
        end
    end
    EE=squeeze(E(sw,2,:))';
    subplot(3,1,1)
    pcolor(time_instances,EE,log10(H)); shading flat; colorbar;
    datetick('x'); set(gca,'Yscale','log');
    subplot(3,1,2)
    pcolor(time_instances,EE,log10(alpha)); shading flat; colorbar;
    datetick('x'); set(gca,'Yscale','log');
    subplot(3,1,3)
    pcolor(time_instances,EE,log10(Heavy)); shading flat; colorbar;
    datetick('x'); set(gca,'Yscale','log');
    
    
    %% Some more samples where we add all data to get an average spectrogram
    pl=6;
    el=8;
    ix = pacclevel==pl & elev(1,:)==el;
    aa=find(masschannel(1,:)==0 & ix);
    Htot=sum(ionspectra(:,aa),2)/length(aa);
    aa=find(masschannel(1,:)==4 & ix);
    alphatot=sum(ionspectra(:,aa),2)/length(aa);
    aa=find(masschannel(1,:)==2 & ix);
    Heavytot=sum(ionspectra(:,aa),2)/length(aa);
    aa=find(masschannel(1,:)==3 & ix);
    Hetot=sum(ionspectra(:,aa),2)/length(aa);
    aa=find(masschannel(1,:)==5 & ix);
    Opptot=sum(ionspectra(:,aa),2)/length(aa);
    figure(22)
    semilogx(...
        EE,Htot,...
        EE,alphatot,...
        EE,Heavytot,...
        EE,Hetot,...
        EE,Opptot)
    
    
    
    %% AN example to create a mass matrix plot ( energy vs mass channel number)
    % Useful for burst mode data when mtable = 0
    
    mass_matrix=zeros(96,32); % sum over all times, azimuths, assume full resolution
    ix=mtable==0;
    if ~isempty(ix)
        for mch=0:31
            aa=find(masschannel(1,:)==mch & ix);
            mass_matrix(:,mch+1)=sum(ionspectra(:,aa),2);
        end
        figure(33)
        pcolor(0:31,EE,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=iontime==time_instances(jj); %
            spectrogram(:,jj)=sum(ionspectra(:,ix),2);
        end
        figure(44)
        pcolor(time_instances,EE,log10(spectrogram)), shading flat, colorbar, set(gca,'yscale','log')
        datetick('x')
    else
        disp('No mass matrix data to show')
    end
    
else
    disp('No data to plot. To read data run READICA2.M.');
end
