%
% This script plots the *.xtr output from Anubis, which is an open source 
% GNSS observable analysis program. For details on Anubis see this link:
% https://www.epos-ip.org/sites/default/files/2017-EUREF-LAC-AnubisTutorial.pdf
%
% See command_line.sh in this directory for a simple script to process
% a T04 file.
% 
% The code was inspired by these two scripts:
% https://github.com/spanikp/anubis-matlab
% http://www.westernexplorers.us/GNSSplotters/
%
% Copyright Trimble Inc. 2019
%

%filename = 'BD935_1_ROVER-20191110-All.xtr';
%filename = 'Z3G_KDDI_1_1123.xtr';
filename = 'Z3G_KDDI_2_1123.xtr';


% When 1 plots C/No instead of MPX
doSNR = 0;

if(doSNR)
  field = 'S';
else
  % Multipath
  field = 'M';
end

fid = fopen(filename,'r');
raw = textscan(fid,'%s','Delimiter','\n','whitespace','');

MaxPRN = 36; % ToDo - see if we can get Anubis to process BeiDou PRN > 36
data = raw{1,1};

% Get a list of all the Multipath (or SNR) observable sets
MPXTypes = {};
for i=1:length(data)
  line = data(i);
  str = line{1};
  if((length(str) >= 7) & (str(1) == '=') & (str(5) == field) & str2num(str(6)))
    MPXTypes{end+1} = str(2:7);
  end
end

MetricMean = [];
% Loop around for each observable type and plot it
for types=1:length(MPXTypes)
  ElevData = {};
  AziData  = {};
  index = 1;
  ElAzMP = [];
  pastHeader = 0;
  
  disp(MPXTypes{types});
  CurrSys = MPXTypes{types}(1:3);

  for i=1:length(data)
    line = data(i);
    str = line{1};

    if(pastHeader == 0)
      if( (length(str) >= 27) & (strcmp(str(1:27),'#====== Elevation & Azimuth')) )
        pastHeader = 1;
      end
    elseif((length(str) >= 7) & (strcmp(str(1:7),[' ',CurrSys,'ELE'])) )
      % Get the elevation mapping leaving this as a line string
      ElevData{end+1} = str;
    elseif((length(str) >= 7) & (strcmp(str(1:7),[' ',CurrSys,'AZI'])) )
      % Get the azimuth mapping leaving this as a line string
      AziData{end+1} = str;
    elseif((length(str) >= 7) & (strcmp(str(1:7),['=',MPXTypes{types}])) )
      % Extract the mean values
      [elem, status] = str2num(str(28:35));
      if(status)
        MetricMean(types) = elem;
      else
        MetricMean(types) = NaN;
      end
    elseif((length(str) >= 7) & (strcmp(str(1:7),[' ',MPXTypes{types}])) )
      % We got a line of "metric" data. "j" loops for each satellite. If
      % we have data for the satellite match it with elevation and azimuth
      % data (el/az is early in the file).
      for j=36:4:(36+4*MaxPRN -1)
        [elem, status] = str2num(str(j:(j+3)));
        if(status)
          % We have MP data - now get the elev / azi
          strEl = ElevData{index};
          strAz = AziData{index};
          [elemEl, statusEl] = str2num(strEl(j:(j+3)));
          [elemAz, statusAz] = str2num(strAz(j:(j+3)));
          % Check we got numbers, remove 0 el / 0 az one data had some
          % corrupt data at this location for GLONASS???!!
          if(statusEl & statusAz & ~(elemEl == 0 & elemAz == 0))
            ElAzMP(end+1,1:4) = [j/4-8, elem, elemEl, elemAz];
          end
        end
      end
      index = index + 1;
    end
  end

  % We've parsed all the data. Make sure we have enough data to process further
  if(length(ElAzMP) > 50)
    s = scatteredInterpolant(ElAzMP(:,4), ElAzMP(:,3), ElAzMP(:,2),'linear','none');
    if(length(s.Points) > 50 )
      [azig, eleg] = meshgrid(0:360, 0:90);
      mpg = s(azig,eleg);

      % The interpolation can cover most of the "Northern hole", check the elevation/azimuth
      % of the data and invalidate interpolated data where there's no actual data in this 
      % region. Perform this in five degree steps to avoid challenges with sparse data.
      %
      % NOTE: This will fail for Southern hemisphere data!!!
      for el=5:5:90
        i = find( (ElAzMP(:,3) >= (el-5)) & (ElAzMP(:,3) < el) );
        if(length(i) > 10)
          az = min(ElAzMP(i,4));
          for j=0:4
            mpg(el-j,1:az) = NaN;
            mpg(el-j,(360-az):end) = NaN;
          end
        end
      end


      %figure
      %polar(ElAzMP(:,4) * pi / 180, 90 - ElAzMP(:,3),'.');

      % Try to color the data based on the value in a simple polar plot
      figure
      cn = ceil(max(ElAzMP(:,2)));
      cm = colormap(parula(cn));
      polarscatter(ElAzMP(:,4) * pi / 180, 90 - ElAzMP(:,3),[],cm(fix(ElAzMP(:,2))+1),'filled');
      % Get the 0 azimuth at the top, we also need to change the
      % direction to clockwise.
      set(gca,'ThetaZeroLocation','top')
      set(gca,'ThetaDir','clockwise')
      if(doSNR)
        title([MPXTypes{types},' - Mean = ',num2str(MetricMean(types)),'dBHz'])
      else
        title([MPXTypes{types},' - Mean = ',num2str(MetricMean(types)),'cm'])
      end
      %figure
      %surf(mpg)
      %view(0,90)
      %colorbar

      figure
      polarplot3d(flipud(mpg),'PlotType','surfn','RadialRange',[0 90],'PolarGrid',{6,12},...
                              'GridStyle',':','AxisLocation','surf','TickSpacing',15);

      MPMax = 20*ceil(max(ElAzMP(:,2)) / 20);

      view(90,-90)
      c = colorbar;
      colLimits = [0,MPMax];
      colLimits(1) = colLimits(1) + 5;
      c.Limits = colLimits;
      c.Ticks = 0:20:MPMax;

      %c.Position = [c.Position(1)*1.02, c.Position(2)*1.4, 0.8*c.Position(3), c.Position(4)*0.9];
      %c.TickDirection = 'in';
      %c.LineWidth = 1.1;
      %c.FontSize = 10;
      %caxis([0,MPMax])
      %ylabel(c,sprintf('%s RMS MP%s value (cm)',GNScell{i},MPcode),'fontsize',10,'fontname','arial')
          
      axis equal
      axis tight
      set(gca,'XTickLabel',[]);
      set(gca,'YTickLabel',[]);
      hold on
      text(60,0,-100,'30','FontSize',10,'HorizontalAlignment','center','background','w','fontname','arial','FontWeight','bold')
      text(30,0,-100,'60','FontSize',10,'HorizontalAlignment','center','background','w','fontname','arial','FontWeight','bold')
      c = colorbar;
      if(doSNR)
        title([MPXTypes{types},' - Mean = ',num2str(MetricMean(types)),'dBHz'])
        ylabel(c,[MPXTypes{types},' [dBHz]']);
      else
        title([MPXTypes{types},' - Mean = ',num2str(MetricMean(types)),'cm'])
        ylabel(c,[MPXTypes{types},' [cm]']);
      end
    end
  end
end      
