% Matlab script to plot ICA NRM data (mass look-up table)
% Note that it only works for mass lookup table data
% For fully resolved mass channels the values
% will be between 0 and 31 and one must use the 
% calibration mass lookup tables (in the CALIB directory)
% to find the mass channels which corresponds to different
% species.
% 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 6 (default)
% or 2. During Earth flyby (20050301) an alternating post acceleration
% mode was used, with 2 and 6 as reference levels.

if exist('iontime')

figure(11)
clf
set(gcf,'renderer','painters');
H=[]; Hpp=[]; O=[]; He=[]; Opp=[];
ICAINITIALIZE
% The spectrograms for different masses and sectors have the same time so
% find all unique time instances:
time_instances=iontime(find(diff(iontime)~=0)); 

pl=2; el=0;
for jj=1:length(time_instances)
  aa=find(masschannel(1,:)==0 & iontime==time_instances(jj) & ...
          pacclevel==pl & elev(1,:)==el);
  if ~isempty(aa)
    H(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
  else
    H(:,jj)=ones(96,1)*NaN;
  end
  aa=find(masschannel(1,:)==4 & iontime==time_instances(jj) & ...
          pacclevel==pl & elev(1,:)==el);
  if ~isempty(aa)
    Hpp(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
  else
    Hpp(:,jj)=ones(96,1)*NaN;
  end
   aa=find(masschannel(1,:)==2 & iontime==time_instances(jj) & ...
           pacclevel==pl & elev(1,:)==el);
  if ~isempty(aa)
    O(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
  else
    O(:,jj)=ones(96,1)*NaN;
  end 
   aa=find(masschannel(1,:)==3 & iontime==time_instances(jj) & ...
           pacclevel==pl & elev(1,:)==el);
  if ~isempty(aa)
    He(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
  else
    He(:,jj)=ones(96,1)*NaN;
  end   
   aa=find(masschannel(1,:)==5 & iontime==time_instances(jj) & ...
           pacclevel==pl & elev(1,:)==el);
  if ~isempty(aa)
    Opp(:,jj)=sum(ionspectra(:,aa),2)./length(aa);
  else
    Opp(:,jj)=ones(96,1)*NaN;
  end   
end
subplot(3,1,1)
pcolor(time_instances,E,log10(H)), shading flat, colorbar
datetick('x'); set(gca,'Yscale','log');
subplot(3,1,2)
pcolor(time_instances,E,log10(Hpp)), shading flat, colorbar
datetick('x'); set(gca,'Yscale','log');
subplot(3,1,3)
pcolor(time_instances,E,log10(O)), shading flat, colorbar
datetick('x'); set(gca,'Yscale','log');


% Some more samples where we add all data to get an average spectrogram
el=8;
pl=2;
aa=find(masschannel(1,:)==0 & pacclevel== pl & elev(1,:)==el);
Htot=sum(ionspectra(:,aa),2)/length(aa); 
aa=find(masschannel(1,:)==4 & pacclevel== pl & elev(1,:)==el);
Hpptot=sum(ionspectra(:,aa),2)/length(aa);  
aa=find(masschannel(1,:)==2 & pacclevel== pl & elev(1,:)==el);
Otot=sum(ionspectra(:,aa),2)/length(aa);
aa=find(masschannel(1,:)==3 & pacclevel== pl & elev(1,:)==el);
Hetot=sum(ionspectra(:,aa),2)/length(aa);
aa=find(masschannel(1,:)==5 & pacclevel== pl & elev(1,:)==el);
Opptot=sum(ionspectra(:,aa),2)/length(aa);  
  figure(22)
  semilogx(E,Htot,E,Hpptot,E,Otot,E,Hetot,E,Opptot)

% Same for default post acceleration level  
figure(23)
pl=6;
aa=find(masschannel(1,:)==0 & pacclevel== pl & elev(1,:)==el);
Htot=sum(ionspectra(:,aa),2)/length(aa); 
aa=find(masschannel(1,:)==4 & pacclevel== pl & elev(1,:)==el);
Hpptot=sum(ionspectra(:,aa),2)/length(aa);  
aa=find(masschannel(1,:)==2 & pacclevel== pl & elev(1,:)==el);
Otot=sum(ionspectra(:,aa),2)/length(aa);
aa=find(masschannel(1,:)==3 & pacclevel== pl & elev(1,:)==el);
Hetot=sum(ionspectra(:,aa),2)/length(aa);
aa=find(masschannel(1,:)==5 & pacclevel== pl & elev(1,:)==el);
Opptot=sum(ionspectra(:,aa),2)/length(aa);  
  semilogx(E,Htot,E,Hpptot,E,Otot,E,Hetot,E,Opptot)

% Same for other elevation hemisphere  
el=0;
pl=2;
aa=find(masschannel(1,:)==0 & pacclevel== pl & elev(1,:)==el);
Htot=sum(ionspectra(:,aa),2)/length(aa); 
aa=find(masschannel(1,:)==4 & pacclevel== pl & elev(1,:)==el);
Hpptot=sum(ionspectra(:,aa),2)/length(aa);  
aa=find(masschannel(1,:)==2 & pacclevel== pl & elev(1,:)==el);
Otot=sum(ionspectra(:,aa),2)/length(aa);
aa=find(masschannel(1,:)==3 & pacclevel== pl & elev(1,:)==el);
Hetot=sum(ionspectra(:,aa),2)/length(aa);
aa=find(masschannel(1,:)==5 & pacclevel== pl & elev(1,:)==el);
Opptot=sum(ionspectra(:,aa),2)/length(aa);  
  figure(12)
  semilogx(E,Htot,E,Hpptot,E,Otot,E,Hetot,E,Opptot)


figure(13)
pl=6; el=0;
aa=find(masschannel(1,:)==0 & pacclevel== pl & elev(1,:)==el);
Htot=sum(ionspectra(:,aa),2)/length(aa); 
aa=find(masschannel(1,:)==4 & pacclevel== pl & elev(1,:)==el);
Hpptot=sum(ionspectra(:,aa),2)/length(aa);  
aa=find(masschannel(1,:)==2 & pacclevel== pl & elev(1,:)==el);
Otot=sum(ionspectra(:,aa),2)/length(aa);
aa=find(masschannel(1,:)==3 & pacclevel== pl & elev(1,:)==el);
Hetot=sum(ionspectra(:,aa),2)/length(aa);
aa=find(masschannel(1,:)==5 & pacclevel== pl & elev(1,:)==el);
Opptot=sum(ionspectra(:,aa),2)/length(aa);
  semilogx(E,Htot,E,Hpptot,E,Otot,E,Hetot,E,Opptot)

else
    disp('No data to plot. To read data run READICA2.M.');
end
