%
% should be called from doDist.m which sets up the names of files etc to
% be processed.
%
% This script reads raw and filtered attitude solution from magnetometer(roll, pitch, and heading)
% and compare to truth.
% It computes the attitude error and plots as CDFs.
%

clearvars -global

% POS Pac Column index
dTime = 1;
dLat  = 2;
dLon  = 3;
dHgt  = 4;

if(exist('TruthLabel') == 0)
  TruthLabel = ['POSPac - ',Section];
else
  TruthLabel = [TruthLabel,' - ',Section];
end

%
% Following are newer configuration parameters make sure that older
% configuration files still work by setting parameters to default values
% if they don't exist.
%
if(exist('saveFig') == 0)
  saveFig = 0;
end

if(exist('addPercentiles') == 0)
  addPercentiles = 0;
end

if(exist('RollPitchCDFAxisMax') == 0)
  % Config parameter name changed
  RollPitchCDFAxisMax = RollPitchAxisMax;
end

if(exist('HeadingCDFAxisMax') == 0)
  % Config parameter name changed
  HeadingCDFAxisMax = HeadingAxisMax;
end

if(exist('HeadAxisMax') == 0)
  % Different name in earlier "CallDoDist.m" scripts
  if(exist('HeadingAxisMax'))
    HeadAxisMax = HeadingAxisMax;
  else
    HeadAxisMax = 1;
  end
end

if(exist('RPAxisMax') == 0)
  RPAxisMax = 0.5;
end

if(exist('plotMap') == 0)
  plotMap = 0;
end

if(exist('zoomLink') == 0)
  zoomLink = 0;
end

if(exist('ProcessOutliers') == 0)
  ProcessOutliers = 0;
else
  % ProcessOutliers may be enabled, make sure a threshold is enabled
  if(exist('d2DOutlierThres') == 0)
    d2DOutlierThres = -1;
  end
  if(exist('dHgtOutlierThres') == 0)
    dHgtOutlierThres = -1;
  end
  if(exist('d2DOutlierPer') == 0)
    d2DOutlierPer = -1;
  end
  if(exist('dHgtOutlierPer') == 0)
    dHgtOutlierPer = -1;
  end
end

% Set the POSPac indices for optional data
dEastVel = 5;
dNorthVel = 6;
dUpVel = 7;
dEastSig  = 8;
dNorthSig = 9;
dUpSig    = 10;
dRoll = 11;
dPitch = 12;
dHeading = 13;

% If any optional data is not vailable set the corresponding index to NaN
i = find (isnan(POSPac(:,dEastVel)) == 0);
if (length(i )== 0)
  dEastVel = NaN;
end
  
i = find (isnan(POSPac(:,dNorthVel)) == 0);
if (length(i )== 0)
  dNorthVel = NaN;
end
  
i = find (isnan(POSPac(:,dUpVel)) == 0);
if (length(i )== 0)
  dUpVel = NaN;
end
  
i = find (isnan(POSPac(:,dEastSig)) == 0);
if (length(i )== 0)
  dEastSig = NaN;
end
  
i = find (isnan(POSPac(:,dNorthSig)) == 0);
if (length(i )== 0)
  dNorthSig = NaN;
end
  
i = find (isnan(POSPac(:,dUpSig)) == 0);
if (length(i )== 0)
  dUpSig = NaN;
end
  
i = find (isnan(POSPac(:,dRoll)) == 0);
if (length(i )== 0)
  dRoll = NaN;
end
  
i = find (isnan(POSPac(:,dPitch)) == 0);
if (length(i )== 0)
  dPitch = NaN;
end
  
i = find (isnan(POSPac(:,dHeading)) == 0);
if (length(i )== 0)
  dHeading = NaN;
end
   
if(length(InputDataName) <= 10)
  % Canned colour array if we have 10 or less data sets to process
 % colorArray(1,1:3) = [1.0   0.75  0.0 ];
 % colorArray(2,1:3) = [1.0   0.0  0.0 ];
 % colorArray(3,1:3) = [0.0   0.0  1.0 ];

  colorArray(1,1:3) = [0.0   0.0  1.0 ];
  colorArray(2,1:3) = [1.0   0.0  0.0 ];
  colorArray(3,1:3) = [1.0   0.75  0.0 ];

  colorArray(4,1:3) = [0.0   1.0  0.0 ];
  colorArray(5,1:3) = [0.0   1.0  1.0 ];
  colorArray(6,1:3) = [1.0   0.0  1.0 ];
  colorArray(7,1:3) = [0.0   0.0  0.0 ];
  colorArray(8,1:3) = [1.0   1.0  0.0 ];
  colorArray(9,1:3) = [1.0   0.0  0.6 ];
  colorArray(10,1:3) = [0.0   0.5  0.0 ];
else
  % Attempt to form a colour array with as much difference between the
  % colours as possible
  colorArray = distinguishable_colors(length(InputDataName));
end


% Add NaN's if there is a gap (assume 10Hz data). The NaN forces a
% pen-up/pen-down
PPOut = penupGaps(POSPac,dTime,0.11);

%
% To give the data some context plot the raw truth.
%
figure
plot(PPOut(:,dLon), PPOut(:,dLat));
xlabel('Longitude [degree]');
ylabel('Latitude [degree]');
axis square;
grid
title([TruthLabel,' Position']);
outputFilename = [Section , ['_',TruthLabel,'_Position.png']];
savePlot(outputFilename,saveFig);

if(plotMap)
  % Load in map data as a background
  % - Note without a google map API key there may be a limit 
  % to the number of calls to this per day
  %
  % This can take a few seconds to download the map data from the server
  h = figure;
  set(h,'PaperPositionMode','auto');
  plot(PPOut(:,dLon), PPOut(:,dLat),'r','linewidth',3);
  disp('Please wait loading google map ...');
  plot_google_map('maptype','roadmap')
  title([TruthLabel,' Position']);
  xlabel('Longitude [degree]');
  ylabel('Latitude [degree]');
  outputFilename = [Section , ['_',TruthLabel,'_Position_Map.png']];
  savePlot(outputFilename,saveFig);

  h = figure;
  set(h,'PaperPositionMode','auto');
  plot(PPOut(:,dLon), PPOut(:,dLat),'r','linewidth',3);
  disp('Please wait loading google map ...');
  plot_google_map('maptype','satellite')
  title([TruthLabel,' Position']);
  xlabel('Longitude [degree]');
  ylabel('Latitude [degree]');
  outputFilename = [Section , ['_',TruthLabel,'_Position_satellite.png']];
  savePlot(outputFilename,saveFig);
end

%
% If we have attitude from the truth plot it
%
if(isnan(dRoll) == 0)
  figure
  plot(PPOut(:,dTime), PPOut(:,dRoll));
  hold on;
  plot(PPOut(:,dTime), PPOut(:,dPitch),'r');
  xlabel('GPS Time [seconds]');
  ylabel('Roll & Pitch [deg]');
  grid
  title([TruthLabel,' Roll & Pitch']);
  legend('Roll','Pitch','Location','Best');
  outputFilename = [Section , ['_',TruthLabel,'_RollPitch.png']];
  savePlot(outputFilename,saveFig);

  figure
  plot(PPOut(:,dTime), PPOut(:,dHeading));
  xlabel('GPS Time [seconds]');
  ylabel('Heading [deg]');
  grid
  title([TruthLabel,' Heading']);
  outputFilename = [Section , ['_',TruthLabel,'_Heading.png']];
  savePlot(outputFilename,saveFig);
end

clear PPOut;
clear CrossTrackData;

% Log command window output
diary OFF;
DiaryFileName = sprintf('%s_stats.txt', Section);
if exist(DiaryFileName,'file')
    delete (DiaryFileName);
end
diary ON;
diary (DiaryFileName);

%
% Loop around reading data and computing the 2D and Cross-track error
% for each data set.
%
for k=1:length(InputDataName)
    clear ENU;
    clear InputData;

    Ftdata = 0;
    DataStr = InputDataName{k}{2};
    
    disp('***************************************************************');
    disp('***************************************************************');
    disp(['                ',DataStr]);
    disp('***************************************************************');
    disp('***************************************************************');
      
    % Load the input data
    FullPathAndName = [InputDataPath '\' InputDataName{k}{3} '\' InputDataName{k}{1}];

    % First try to read a matlab version of the file, if that fails read
    % the original file and then write a matlab format. By caching to
    % .mat file format we can make the read over 100 times faster,
    % particularly useful if we are reading the same data file in
    % multiple times! Original name is of the format ABCDEF.XXX the mat
    % format ABCDEF.XXX.mat will be generated. We use the name format
    % .dut.mat to avoid issues if the reference is changed to be a DUT,
    % the reference uses a cached MAT format but has different data
    % format, we don't want to pick up an old reference file here
    if(exist([FullPathAndName,'.dut.mat'],'file') == 2)
      disp(['!!!!Warning MAT format ',FullPathAndName,' file - erase if newer TXT data available!!!!']);
      load([FullPathAndName,'.dut.mat']);
    end

    if(exist('InputData') == 0)
      % A mat version of the file doesn't exist or didn't contain the
      % correct variable, read in the data and save as a MAT file for next time.

      if (strcmpi(InputDataName{k}{4},'TXT'))
        InputData = load(FullPathAndName);        
        save([FullPathAndName,'.dut.mat'],'InputData');
      elseif (strcmpi(InputDataName{k}{4},'MAT'))
        % Matlab format from viewdat (--mat)
        InputData = load_viewdat(FullPathAndName);
      else
        error('Unknown file format');    
      end
    end

    % If we have defined start/stop bounds compute the yield and
    % datarate
    clear ActualCount;
    clear ExpectedCount;
    ActualCount(1) = NaN;
    for m=1:length(StartTime)
      % We can have multiple segments of data, loop around for each
      % segment and build up the yield.
      if (StartTime(m) >= 0 && EndTime(m) >= 0)
        TimeIndex = InputDataName{k}{5};
        if (m == 1)
            InputData(:,TimeIndex) = InputData(:,TimeIndex)/1000;
        end
        i = find( (InputData(:,TimeIndex) >= StartTime(m)) & (InputData(:,TimeIndex) <= EndTime(m)));
        ActualCount(m) = length(i);
        if ActualCount(m) > 0
           % Compute the expected epoch count
           dataRate(k) = min(diff(InputData(i,TimeIndex)));
           
           % Buffalo data had some strange time tags (e.g. 0.91
           % seconds). Find the closest time tag to expected data rates
           rates = [0.1 0.2 0.5 1 2 5 10];

           MinDelta = 100;
           for i = 1:length(rates)
             delta = abs((dataRate(k) -  rates(i)));
             if(delta < MinDelta)
               MinDelta = delta;
             end
           end

           % Don't always get zero due to numerical rounding
           if(MinDelta > 1e-6)
             disp(['Unexpected data rate = ',num2str(dataRate(k)),'s - re-estimating']);
             MinDelta = 100;
             MinIndex =   0;
             for i = 1:length(rates)
               delta = abs(dataRate(k) - rates(i));
               if(delta < MinDelta)
                 MinDelta = delta;
                 MinIndex = i;
               end
             end
             if(MinDelta ~= 100)
               dataRate(k) = rates(MinIndex);
               disp(['Estimated data rate = ',num2str(dataRate(k))]);
             end
           end

           ExpectedCount(m) = (EndTime(m) - StartTime(m)) / dataRate(k) + 1;
        end       
      end
    end

    if(sum(ActualCount) > 0)
      Yield(k) = 100.0*(sum(ActualCount) / sum(ExpectedCount));
      disp([DataStr, ' Solution availability : ', num2str(Yield(k)), '%']);
    else
      Yield(k) = NaN;
      continue;
    end
    
    % Do we need to round the time tag due to 1ms offsets?
    % Some of the consumer devices occasionally have time tags off by
    % 1ms
    if(InputDataName{k}{7} == 1) 
        % Time tag rounding (often have 1ms offset)
        InputData(:,InputDataName{k}{5}) = round(InputData(:,InputDataName{k}{5}).*10)./10;
    end    
    
    % Load common data points with matching time tags
    [C, A, B] = intersect(InputData(:,InputDataName{k}{5}), POSPac(:,dTime));


    % Realtime data == InputData(A,:);
    % POSPac   data == POSPac(B,:);

    %--------------------------------------------------------------------
    % Merge the DUT (device under test) and truth data time aligned into
    % a common 2D array.
    %--------------------------------------------------------------------

    key.time       =  1;
    key.truth_lat  =  2;
    key.truth_lon  =  3;
    key.truth_hgt  =  4;
    key.truth_eVel =  5;
    key.truth_nVel =  6;
    key.truth_uVel =  7;
    key.DUT_Roll       =  8;
    key.DUT_Pitch      =  9;
    key.DUT_Heading    = 10;
    key.truth_Roll     = 11;
    key.truth_Pitch    = 12;
    key.truth_Heading  = 13;
    key.DUT_sigRoll    = 14;
    key.DUT_sigPitch   = 15;
    key.DUT_sigHeading = 16;
    key.DUT_FtRoll    = 17;
    key.DUT_FtPitch   = 18;
    key.DUT_FtHeading = 19;
    key.DUT_FtSigRoll    = 20;
    key.DUT_FtSigPitch   = 21;
    key.DUT_FtSigHeading = 22;

    clear Merge;

    % Position, position stats
    Merge(:,key.time)    = InputData(A,InputDataName{k}{5});
    Merge(:,key.truth_lat) = POSPac(B,dLat);
    Merge(:,key.truth_lon) = POSPac(B,dLon);
    Merge(:,key.truth_hgt) = POSPac(B,dHgt);
    
    % Velocity
    if (isnan(dEastVel) == 0)
      Merge(:,key.truth_eVel) = POSPac(B,dEastVel);
    else
      Merge(:,key.truth_eVel) = NaN(1,length(B));
    end
    
    if (isnan(dNorthVel) == 0)
      Merge(:,key.truth_nVel) = POSPac(B,dNorthVel);
    else
      Merge(:,key.truth_nVel) = NaN(1,length(B));
    end
    
    if (isnan(dUpVel) == 0)
      Merge(:,key.truth_uVel) = POSPac(B,dUpVel);
    else
      Merge(:,key.truth_uVel) = NaN(1,length(B));
    end

    % Attitude
    % Raw mag attitude
    if(isnan(InputDataName{k}{8}) == 0)
      Merge(:,key.DUT_Roll) = InputData(A,InputDataName{k}{8});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_Roll) = Merge(:,key.DUT_Roll)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_Roll) = NaN(1,length(A));
    end
    
    if(isnan(InputDataName{k}{9}) == 0)
      Merge(:,key.DUT_Pitch) = InputData(A,InputDataName{k}{9});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_Pitch) = Merge(:,key.DUT_Pitch)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_Pitch) = NaN(1,length(A));
    end

    if(isnan(InputDataName{k}{10}) == 0)
      Merge(:,key.DUT_Heading) = InputData(A,InputDataName{k}{10});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_Heading) = Merge(:,key.DUT_Heading)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_Heading) = NaN(1,length(A));
    end
    
    % Filtered mag attitude
    if(isnan(InputDataName{k}{17}) == 0)
      Merge(:,key.DUT_FtRoll) = InputData(A,InputDataName{k}{17});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_FtRoll) = Merge(:,key.DUT_FtRoll)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_FtRoll) = NaN(1,length(A));
    end
    
    if(isnan(InputDataName{k}{18}) == 0)
      Merge(:,key.DUT_FtPitch) = InputData(A,InputDataName{k}{18});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_FtPitch) = Merge(:,key.DUT_FtPitch)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_FtPitch) = NaN(1,length(A));
    end

    if(isnan(InputDataName{k}{19}) == 0)
      Merge(:,key.DUT_FtHeading) = InputData(A,InputDataName{k}{19});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_FtHeading) = Merge(:,key.DUT_FtHeading)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_FtHeading) = NaN(1,length(A));
    end
    
    % True attitude
    if(isnan(dRoll) == 0)
      Merge(:,key.truth_Roll) = POSPac(B,dRoll);
    else
      Merge(:,key.truth_Roll) = NaN(1,length(B));     
    end
    
    if(isnan(dPitch) == 0)
      Merge(:,key.truth_Pitch)   = POSPac(B,dPitch);
    else
      Merge(:,key.truth_Pitch) = NaN(1,length(B));     
    end
    
    if(isnan(dHeading) == 0)
      Merge(:,key.truth_Heading) = POSPac(B,dHeading);
    else
      Merge(:,key.truth_Heading) = NaN(1,length(B));     
    end

    % Raw attitude sigma
    if(isnan(InputDataName{k}{14}) == 0)
      Merge(:,key.DUT_sigRoll) = InputData(A,InputDataName{k}{14});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_sigRoll) = Merge(:,key.DUT_sigRoll)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_sigRoll) = NaN(1,length(A));
    end
    
    if(isnan(InputDataName{k}{15}) == 0)
      Merge(:,key.DUT_sigPitch) = InputData(A,InputDataName{k}{15});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_sigPitch) = Merge(:,key.DUT_sigPitch)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_sigPitch) = NaN(1,length(A));
    end

    if(isnan(InputDataName{k}{16}) == 0)
      Merge(:,key.DUT_sigHeading) = InputData(A,InputDataName{k}{16});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_sigHeading) = Merge(:,key.DUT_sigHeading)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_sigHeading) = NaN(1,length(A));
    end
    
    % Filtered attitude sigma
    if(isnan(InputDataName{k}{20}) == 0)
      Merge(:,key.DUT_FtSigRoll) = InputData(A,InputDataName{k}{20});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_FtSigRoll) = Merge(:,key.DUT_FtSigRoll)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_FtSigRoll) = NaN(1,length(A));
    end
    
    if(isnan(InputDataName{k}{21}) == 0)
      Merge(:,key.DUT_FtSigPitch) = InputData(A,InputDataName{k}{21});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_FtSigPitch) = Merge(:,key.DUT_FtSigPitch)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_FtSigPitch) = NaN(1,length(A));
    end

    if(isnan(InputDataName{k}{22}) == 0)
      Merge(:,key.DUT_FtSigHeading) = InputData(A,InputDataName{k}{22});

      if(strcmpi(InputDataName{k}{6},'RAD')) 
        Merge(:,key.DUT_FtSigHeading) = Merge(:,key.DUT_FtSigHeading)  .* 180 / pi;
      end
    else
      Merge(:,key.DUT_FtSigHeading) = NaN(1,length(A));
    end

    %--------------------------------------------------------------------
    % DUT and Truth data are now merged in a common array, use key.* to
    % access the various elements
    %--------------------------------------------------------------------
    
    if(exist('dataRate'))
      localDataRate = dataRate(k);
    else
      % We don't have a data rate because we don't have start/stop
      % times. Based on the smallest time gap estimate the data rate.
      localDataRate = min(diff(Merge(:,key.time)));
    end

    % Insert NaN's into the data if there is a gap to force a pen up in
    % the plots.
    %
    % Caution: When forming any stats need to remove the NaNs, use
    % "GoodTimeTag" to get an index without NaNs.
    MergeOut = penupGaps(Merge,key.time, localDataRate + 0.001);

    % A data set from a uBlox8 had a lot of single issed epochs. "penupGaps"
    % failed and had several "0" time stamps returned but only at the
    % beginning of the data set. The following will fail at week roll,
    % but we don't test then so for now find 0 time tags and change to
    % NaN so they will be filtered out.
    i = find(MergeOut(:,key.time) == 0);
    MergeOut(i,:) = NaN;

    Merge = MergeOut;
    GoodTimeTag = find(isnan(Merge(:,key.time)) == 0);

    clear MergeOut;

    %
    % If we have truth and estimated attitude...
    if(  (isnan(InputDataName{k}{8}) == 0) ...
      && (isnan(InputDataName{k}{9}) == 0) ...
      && (isnan(InputDataName{k}{10}) == 0) ...
      && (isnan(dRoll) == 0) )
        
        % Check for corrections...
        if (isnan(InputDataName{k}{11}) == 0)
            RollCorr = InputDataName{k}{11};
        else
            RollCorr = 0;
        end
        if (isnan(InputDataName{k}{12}) == 0)
            PitchCorr = InputDataName{k}{12};
        else
            PitchCorr = 0;
        end
        if (isnan(InputDataName{k}{13}) == 0)
            HeadingCorr = InputDataName{k}{13};
        else
            HeadingCorr = 0;
        end
        
        % If roll(X), pitch(Y) or heading(Z) rotations are non-zero...
        if (RollCorr ~= 0 || PitchCorr ~= 0 || HeadingCorr ~= 0)
            % Correct the roll,pitch and heading estimates
            [RollHat, PitchHat, HeadingHat] = correctRollPitchHeading( Merge(GoodTimeTag,key.DUT_Roll),     ... % Roll
                                                                       Merge(GoodTimeTag,key.DUT_Pitch),    ... % Pitch
                                                                       Merge(GoodTimeTag,key.DUT_Heading),  ... % Heading
                                                                       RollCorr,  ...                 % Xrotation Correction
                                                                       PitchCorr,  ...                % Yrotation Correction
                                                                       HeadingCorr);                % Yrotation Correction
        else
            RollHat = Merge(GoodTimeTag,key.DUT_Roll);
            PitchHat = Merge(GoodTimeTag,key.DUT_Pitch);
            HeadingHat = Merge(GoodTimeTag,key.DUT_Heading);
        end
        
        
        % now compute the roll/pitch/heading error
        %  i = find(isnan(RollHat) == 0);
        [RollErr, PitchErr, HeadingErr] = compareRollPitchHeading( Merge(GoodTimeTag,key.truth_Roll), ...
                                                                   Merge(GoodTimeTag,key.truth_Pitch), ...
                                                                   Merge(GoodTimeTag,key.truth_Heading), ...
                                                                   RollHat(:)', ...
                                                                   PitchHat(:)', ...
                                                                   HeadingHat(:)');
        
        i = find(isnan(RollErr) == 0);
        HeadingErr = HeadingErr';
        RollErr = RollErr';
        PitchErr = PitchErr';
        
        % Raw roll error w.r.t. POSPac
        figure
        plot(Merge(GoodTimeTag,key.time), RollErr);
        hold on;
        if(isnan(InputDataName{k}{14}) == 0)
            % +/- 1-sigma estimates
            plot(Merge(GoodTimeTag,key.time), RollErr + Merge(GoodTimeTag,key.DUT_sigRoll),'m');
            plot(Merge(GoodTimeTag,key.time), RollErr - Merge(GoodTimeTag,key.DUT_sigRoll),'m');
            title(['Roll Error - ',DataStr,' with 1-sigma estimate bounds - ',Section]);
        else
            title(['Roll Error - ',DataStr,' - ',Section]);
        end
        xlabel('Time [GPS Seconds]');
        ylabel('Error [Deg]');
        Ax = axis;
        Ax(3) = -RPAxisMax;
        Ax(4) =  RPAxisMax;
        axis(Ax);
        grid
        outputFileName = [Section '_' DataStr '_RollErr.png'];
        savePlot(outputFileName,saveFig);
        
        disp(['Roll 1-sigma   = ',num2str(std(RollErr(i)),'%.3f'),'deg']);
        disp(['Roll mean      = ',num2str(mean(RollErr(i)),'%.3f'),'deg']);
        RollErrSort{k} = sort(abs(RollErr(i))); % Sort for the CDF
        
        % Raw pitch error w.r.t. POSPac
        figure
        plot(Merge(GoodTimeTag,key.time), PitchErr);
        hold on;
        if(isnan(InputDataName{k}{15}) == 0)
            % +/- 1-sigma estimates
            plot(Merge(GoodTimeTag,key.time), PitchErr + Merge(GoodTimeTag,key.DUT_sigPitch),'m');
            plot(Merge(GoodTimeTag,key.time), PitchErr - Merge(GoodTimeTag,key.DUT_sigPitch),'m');
            title(['Pitch Error - ',DataStr,' with 1-sigma estimate bounds - ',Section]);
        else
            title(['Pitch Error - ',DataStr,' - ',Section]);
        end
        xlabel('Time [GPS Seconds]');
        ylabel('Error [Deg]');
        grid
        Ax = axis;
        Ax(3) = -RPAxisMax;
        Ax(4) =  RPAxisMax;
        axis(Ax);
        outputFileName = [Section '_' DataStr '_PitchErr.png'];
        savePlot(outputFileName,saveFig);
        
        disp(['Pitch 1-sigma  = ',num2str(std(PitchErr(i)),'%.3f'),'deg']);
        disp(['Pitch mean     = ',num2str(mean(PitchErr(i)),'%.3f'),'deg']);
        PitchErrSort{k} = sort(abs(PitchErr(i))); % Sort for the CDF
        
        % Raw heading error w.r.t. POSPac
        j = find(HeadingErr > 180);
        HeadingErr(j) = HeadingErr(j) - 360;
        figure
        plot(Merge(GoodTimeTag,key.time), HeadingErr);
        hold on;
        if(isnan(InputDataName{k}{16}) == 0)
            % +/- 1-sigma estimates
            plot(Merge(GoodTimeTag,key.time), HeadingErr + Merge(GoodTimeTag,key.DUT_sigHeading),'m');
            plot(Merge(GoodTimeTag,key.time), HeadingErr - Merge(GoodTimeTag,key.DUT_sigHeading),'m')
            title(['Heading Error - ',DataStr,' with 1-sigma estimate bounds - ',Section]);
        else
            title(['Heading Error - ',DataStr,' - ',Section]);
        end
        xlabel('Time [GPS Seconds]');
        ylabel('Error [Deg]');
        grid
        Ax = axis;
        Ax(3) = -HeadAxisMax;
        Ax(4) =  HeadAxisMax;
        axis(Ax);
        outputFileName = [Section '_' DataStr '_HeadingErr.png'];
        savePlot(outputFileName,saveFig);
        
        disp(['Heading 1-sigma= ',num2str(std(HeadingErr(i)),'%.3f'),'deg']);
        disp(['Heading mean   = ',num2str(mean(HeadingErr(i)),'%.3f'),'deg']);
        HeadingErrSort{k} = sort(abs(HeadingErr(i))); % Sort for the CDF
        
        % Raw Roll
        RollData = sort(abs(RollErr(i)));
        RollStats(k,1) = RollData(floor(0.50*length(RollData)));
        RollStats(k,2) = RollData(floor(0.68*length(RollData)));
        RollStats(k,3) = RollData(floor(0.95*length(RollData)));
        RollStats(k,4) = RollData(floor(0.99*length(RollData)));
        RollStats(k,5) = RollData(end);
        % Raw Pitch
        PitchData = sort(abs(PitchErr(i)));
        PitchStats(k,1) = PitchData(floor(0.50*length(PitchData)));
        PitchStats(k,2) = PitchData(floor(0.68*length(PitchData)));
        PitchStats(k,3) = PitchData(floor(0.95*length(PitchData)));
        PitchStats(k,4) = PitchData(floor(0.99*length(PitchData)));
        PitchStats(k,5) = PitchData(end);
        % Raw Heading
        HeadingData = sort(abs(HeadingErr(i)));
        HeadingStats(k,1) = HeadingData(floor(0.50*length(HeadingData)));
        HeadingStats(k,2) = HeadingData(floor(0.68*length(HeadingData)));
        HeadingStats(k,3) = HeadingData(floor(0.95*length(HeadingData)));
        HeadingStats(k,4) = HeadingData(floor(0.99*length(HeadingData)));
        HeadingStats(k,5) = HeadingData(end);
        
    else
        RollStats(k,1:5) = NaN;
        PitchStats(k,1:5) = NaN;
        HeadingStats(k,1:5) = NaN;
    end
   
    % If we have truth and filtered attitude...
    if(  (isnan(InputDataName{k}{17}) == 0) ...
      && (isnan(InputDataName{k}{18}) == 0) ...
      && (isnan(InputDataName{k}{19}) == 0) ...
      && (isnan(dRoll) == 0) )
        Ftdata = 1;
        
        % Check for corrections...
        if (isnan(InputDataName{k}{11}) == 0)
            RollCorr = InputDataName{k}{11};
        else
            RollCorr = 0;
        end
        if (isnan(InputDataName{k}{12}) == 0)
            PitchCorr = InputDataName{k}{12};
        else
            PitchCorr = 0;
        end
        if (isnan(InputDataName{k}{13}) == 0)
            HeadingCorr = InputDataName{k}{13};
        else
            HeadingCorr = 0;
        end
        
        % If roll(X), pitch(Y) or heading(Z) rotations are non-zero...
        if (RollCorr ~= 0 || PitchCorr ~= 0 || HeadingCorr ~= 0)
            % Correct the roll,pitch and heading estimates
            [RollHatFt, PitchHatFt,  HeadingHatFt] = correctRollPitchHeading( Merge(GoodTimeTag,key.DUT_FtRoll),...   % Roll
                                                                              Merge(GoodTimeTag,key.DUT_FtPitch),...  % Pitch
                                                                              Merge(GoodTimeTag,key.DUT_FtHeading),...% Heading
                                                                              RollCorr,  ...                    % Xrotation Correction
                                                                              PitchCorr,  ...                   % Yrotation Correction
                                                                              HeadingCorr);                     % Zrotation Correction
        else
            RollHatFt = Merge(GoodTimeTag,key.DUT_FtRoll);
            PitchHatFt = Merge(GoodTimeTag,key.DUT_FtPitch);
            HeadingHatFt = Merge(GoodTimeTag,key.DUT_FtHeading);
        end
        
        
        % now compute the filtered roll/pitch/heading error
        [FtRollErr, FtPitchErr, FtHeadingErr] = compareRollPitchHeading( Merge(GoodTimeTag,key.truth_Roll), ...
                                                                         Merge(GoodTimeTag,key.truth_Pitch), ...
                                                                         Merge(GoodTimeTag,key.truth_Heading), ...
                                                                         RollHatFt(:)', ...
                                                                         PitchHatFt(:)', ...
                                                                         HeadingHatFt(:)');
        iFt = find(isnan(FtRollErr) == 0);
        FtRollErr = FtRollErr';
        FtPitchErr = FtPitchErr';
        FtHeadingErr = FtHeadingErr';
        
        % Filtered roll error w.r.t. POSPac
        figure
        plot(Merge(GoodTimeTag,key.time), FtRollErr);
        hold on;
        if(isnan(InputDataName{k}{20}) == 0)
            % +/- 1-sigma estimates
            plot(Merge(GoodTimeTag,key.time), FtRollErr + Merge(GoodTimeTag,key.DUT_FtSigRoll),'m');
            plot(Merge(GoodTimeTag,key.time), FtRollErr - Merge(GoodTimeTag,key.DUT_FtSigRoll),'m');
            title(['Filtered Roll Error - ',DataStr,' with 1-sigma estimate bounds - ',Section]);
        else
            title(['Filtered Roll Error - ',DataStr,' - ',Section]);
        end
        xlabel('Time [GPS Seconds]');
        ylabel('Error [Deg]');
        Ax = axis;
        Ax(3) = -RPAxisMax;
        Ax(4) =  RPAxisMax;
        axis(Ax);
        grid
        outputFileName = [Section '_' DataStr '_FilterRollErr.png'];
        savePlot(outputFileName,saveFig);
        
        disp(['Filtered Roll 1-sigma   = ',num2str(std(FtRollErr(iFt)),'%.3f'),'deg']);
        disp(['Filtered Roll mean      = ',num2str(mean(FtRollErr(iFt)),'%.3f'),'deg']);
        FtRollErrSort{k} = sort(abs(FtRollErr(iFt))); % Sort for the CDF
        
        % Filtered pitch error w.r.t. POSPac
        figure
        plot(Merge(GoodTimeTag,key.time), FtPitchErr);
        hold on;
        if(isnan(InputDataName{k}{15}) == 0)
            % +/- 1-sigma estimates
            plot(Merge(GoodTimeTag,key.time), FtPitchErr + Merge(GoodTimeTag,key.DUT_FtSigPitch),'m');
            plot(Merge(GoodTimeTag,key.time), FtPitchErr - Merge(GoodTimeTag,key.DUT_FtSigPitch),'m');
            title(['Filtered Pitch Error - ',DataStr,' with 1-sigma estimate bounds - ',Section]);
        else
            title(['Filtered Pitch Error - ',DataStr,' - ',Section]);
        end
        xlabel('Time [GPS Seconds]');
        ylabel('Error [Deg]');
        grid
        Ax = axis;
        Ax(3) = -RPAxisMax;
        Ax(4) =  RPAxisMax;
        axis(Ax);
        outputFileName = [Section '_' DataStr '_FilterPitchErr.png'];
        savePlot(outputFileName,saveFig);
        
        disp(['Filtered Pitch 1-sigma  = ',num2str(std(FtPitchErr(iFt)),'%.3f'),'deg']);
        disp(['Filtered Pitch mean     = ',num2str(mean(FtPitchErr(iFt)),'%.3f'),'deg']);
        FtPitchErrSort{k} = sort(abs(FtPitchErr(iFt))); % Sort for the CDF
        
        % Filtered heading error w.r.t. POSPac
        j = find(FtHeadingErr > 180);
        FtHeadingErr(j) = FtHeadingErr(j) - 360;
        figure
        plot(Merge(GoodTimeTag,key.time), FtHeadingErr);
        hold on;
        if(isnan(InputDataName{k}{16}) == 0)
            % +/- 1-sigma estimates
            plot(Merge(GoodTimeTag,key.time), FtHeadingErr + Merge(GoodTimeTag,key.DUT_FtSigHeading),'m');
            plot(Merge(GoodTimeTag,key.time), FtHeadingErr - Merge(GoodTimeTag,key.DUT_FtSigHeading),'m')
            title(['Filtered Heading Error - ',DataStr,' with 1-sigma estimate bounds - ',Section]);
        else
            title(['Filtered Heading Error - ',DataStr,' - ',Section]);
        end
        xlabel('Time [GPS Seconds]');
        ylabel('Error [Deg]');
        grid
        Ax = axis;
        Ax(3) = -HeadAxisMax;
        Ax(4) =  HeadAxisMax;
        axis(Ax);
        outputFileName = [Section '_' DataStr '_FilterHeadingErr.png'];
        savePlot(outputFileName,saveFig);
        
        disp(['Filtered Heading 1-sigma= ',num2str(std(FtHeadingErr(iFt)),'%.3f'),'deg']);
        disp(['Filtered Heading mean   = ',num2str(mean(FtHeadingErr(iFt)),'%.3f'),'deg']);
        FtHeadingErrSort{k} = sort(abs(FtHeadingErr(iFt))); % Sort for the CDF

        % Filtered Roll
        RollData = sort(abs(FtRollErr(iFt)));
        FtRollStats(k,1) = RollData(floor(0.50*length(RollData)));
        FtRollStats(k,2) = RollData(floor(0.68*length(RollData)));
        FtRollStats(k,3) = RollData(floor(0.95*length(RollData)));
        FtRollStats(k,4) = RollData(floor(0.99*length(RollData)));
        FtRollStats(k,5) = RollData(end);
        %Filtered Pitch
        PitchData = sort(abs(FtPitchErr(iFt)));
        FtPitchStats(k,1) = PitchData(floor(0.50*length(PitchData)));
        FtPitchStats(k,2) = PitchData(floor(0.68*length(PitchData)));
        FtPitchStats(k,3) = PitchData(floor(0.95*length(PitchData)));
        FtPitchStats(k,4) = PitchData(floor(0.99*length(PitchData)));
        FtPitchStats(k,5) = PitchData(end);
        % Filtered Heading
        HeadingData = sort(abs(FtHeadingErr(iFt)));
        FtHeadingStats(k,1) = HeadingData(floor(0.50*length(HeadingData)));
        FtHeadingStats(k,2) = HeadingData(floor(0.68*length(HeadingData)));
        FtHeadingStats(k,3) = HeadingData(floor(0.95*length(HeadingData)));
        FtHeadingStats(k,4) = HeadingData(floor(0.99*length(HeadingData)));
        FtHeadingStats(k,5) = HeadingData(end);
    else
        FtRollStats(k,1:5) = NaN;
        FtPitchStats(k,1:5) = NaN;
        FtHeadingStats(k,1:5) = NaN;
    end

    % Save the name of this data set for later
    NameArray{k} = DataStr;

end %for k=1:length(InputDataName)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We have data for all receivers, now plot several CDFs to compare

% Plot Roll/Pitch/Heading CDFs
if( (isnan(dRoll) == 0) & (exist('RollErrSort')) )
  createCDF(RollErrSort, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            RollPitchCDFAxisMax, ...
            0.0, ...
            'abs(Roll Error) [deg]', ...
            ['Cumulative Distribution : abs(Roll) CDF (', Section,')'], ...
            [Section, '_RollCDF.'], ...
            saveFig, ...
            addPercentiles);
end

if( (isnan(dPitch) == 0) & (exist('PitchErrSort')) )
  createCDF(PitchErrSort, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            RollPitchCDFAxisMax, ...
            0.0, ...
            'abs(Pitch Error) [deg]', ...
            ['Cumulative Distribution : abs(Pitch) CDF (', Section,')'], ...
            [Section, '_PitchCDF.'], ...
            saveFig, ...
            addPercentiles);
end

if( (isnan(dHeading) == 0) & (exist('HeadingErrSort')) )
  createCDF(HeadingErrSort, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            HeadingCDFAxisMax, ...
            0.0, ...
            'abs(Heading Error) [deg]', ...
            ['Cumulative Distribution : abs(Heading) CDF (', Section,')'], ...
            [Section, '_HeadingCDF.'], ...
            saveFig, ...
            addPercentiles);
end
% Filtered CDFs
if( (isnan(dRoll) == 0) & (exist('FtRollErrSort')) )
  createCDF(FtRollErrSort, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            RollPitchCDFAxisMax, ...
            0.0, ...
            'abs(Roll Error) [deg]', ...
            ['Cumulative Distribution : abs(Roll) CDF (', Section,')'], ...
            [Section, '_FilterRollCDF.'], ...
            saveFig, ...
            addPercentiles);
end

if( (isnan(dPitch) == 0) & (exist('FtPitchErrSort')) )
  createCDF(FtPitchErrSort, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            RollPitchCDFAxisMax, ...
            0.0, ...
            'abs(Pitch Error) [deg]', ...
            ['Cumulative Distribution : abs(Pitch) CDF (', Section,')'], ...
            [Section, '_FilterPitchCDF.'], ...
            saveFig, ...
            addPercentiles);
end

if( (isnan(dHeading) == 0) & (exist('FtHeadingErrSort')) )
  createCDF(FtHeadingErrSort, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            HeadingCDFAxisMax, ...
            0.0, ...
            'abs(Heading Error) [deg]', ...
            ['Cumulative Distribution : abs(Heading) CDF (', Section,')'], ...
            [Section, '_FilterHeadingCDF.'], ...
            saveFig, ...
            addPercentiles);
end

diary OFF;

%
% Collate the key stats as HTML (XXX_stats.html):-
%
name = [Section,'_mag_stats.html'];
fid = fopen(name,'w');
fprintf(fid,'<!DOCTYPE html>\n');
fprintf(fid,'<html>\n');
fprintf(fid,'<head>\n');
fprintf(fid,'<title>%s - Mag Analysis</title>\n',Section);

fprintf(fid,'<style>\n');
fprintf(fid,'table\n');
fprintf(fid,'{\n');
fprintf(fid,'border-collapse:collapse;\n');
fprintf(fid,'}\n');
fprintf(fid,'table, td, th\n');
fprintf(fid,'{\n');
fprintf(fid,'border:1px solid black;\n');
fprintf(fid,'}\n');
fprintf(fid,'</style>\n');

fprintf(fid,'</head>\n');
fprintf(fid,'<body>\n');

fprintf(fid,'<h2>Analysis Directory: %s</h2>',pwd);

if(plotMap)
    File = [Section,'_POSPac - ',Section,'_Position_Map.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
    File = [Section,'_POSPac - ',Section,'_Position_satellite.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
    fprintf(fid,'<br/>');
end

fprintf(fid,'<br/><hr/>');

AttitudeAvailable = length(find(isnan(RollStats) == 0));
if(AttitudeAvailable)
  MaxPlots = 3;
end

fprintf(fid,'<h3>Result Plots</h3>\n');
fprintf(fid,'<table>\n');
fprintf(fid,'<tr><th>Data Set</th>');

fprintf(fid,'<th>Roll Err</th><th>Pitch Err</th><th>Heading Err</th><th>Filter Roll Err</th><th>Filter Pitch Err</th><th>Filter Heading Err</th>');

for i=1:k
    if (~isnan(Yield(i)) && dataRate(i) > 0)
        StatAxis{i} = InputDataName{i}{2};
        
        fprintf(fid,'<tr>\n');
        
        fprintf(fid,['<td>',InputDataName{i}{2},'</td>\n']);
        
        StrPNG = [Section,'_',InputDataName{i}{2},'_RollErr.png'];
        StrFIG = [Section,'_',InputDataName{i}{2},'_RollErr.fig'];
        fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>&nbsp;<a href="',StrFIG,'">FIG</a></td>']);
        
        StrPNG = [Section,'_',InputDataName{i}{2},'_PitchErr.png'];
        StrFIG = [Section,'_',InputDataName{i}{2},'_PitchErr.fig'];
        fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>&nbsp;<a href="',StrFIG,'">FIG</a></td>']);
        
        StrPNG = [Section,'_',InputDataName{i}{2},'_HeadingErr.png'];
        StrFIG = [Section,'_',InputDataName{i}{2},'_HeadingErr.fig'];
        fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>&nbsp;<a href="',StrFIG,'">FIG</a></td>']);
        
        if (Ftdata)
            StrPNG = [Section,'_',InputDataName{i}{2},'_FilterRollErr.png'];
            StrFIG = [Section,'_',InputDataName{i}{2},'_FilterRollErr.fig'];
            fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>&nbsp;<a href="',StrFIG,'">FIG</a></td>']);
            
            StrPNG = [Section,'_',InputDataName{i}{2},'_FilterPitchErr.png'];
            StrFIG = [Section,'_',InputDataName{i}{2},'_FilterPitchErr.fig'];
            fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>&nbsp;<a href="',StrFIG,'">FIG</a></td>']);
            
            StrPNG = [Section,'_',InputDataName{i}{2},'_FilterHeadingErr.png'];
            StrFIG = [Section,'_',InputDataName{i}{2},'_FilterHeadingErr.fig'];
            fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>&nbsp;<a href="',StrFIG,'">FIG</a></td>']);
        else
            fprintf(fid,'<th>N/A</th><th>N/A</th><th>N/A</th>');
        end
        
        fprintf(fid,'</tr>\n');
    end
end

fprintf(fid,'</table>\n');

fprintf(fid,'<hr/>\n');

for s=1:MaxPlots
    if(s == 1)
        fprintf(fid,'<h3>Roll Stats:</h3>\n');
        stats = RollStats;
        ftstats = FtRollStats;
    elseif(s == 2)
        fprintf(fid,'<h3>Pitch Stats:</h3>\n');
        stats = PitchStats;
        ftstats = FtPitchStats;
    elseif(s == 3)
        fprintf(fid,'<h3>Heading Stats:</h3>\n');
        stats = HeadingStats;
        ftstats = FtHeadingStats;
    end
    
    fprintf(fid,'<table>\n');
    fprintf(fid,'<tr><th>Data Set</th>');
    if (StartTime(1) >= 0 && EndTime(1) >= 0)
        fprintf(fid,'<th>Yield [%%]*</th><th>Rate [Hz]</th>');
    end
    
    fprintf(fid,'<th>50%%</th><th>68%%</th><th>95%%</th><th>99%%</th><th>100%%</th>');
    fprintf(fid,'<th>Primary Soln**</th>');
    fprintf(fid,'<th>50%%</th><th>68%%</th><th>95%%</th><th>99%%</th><th>100%%</th>');
    fprintf(fid,'<th>Primary Soln**</th></tr>\n');
    
    
    RowCol = size(stats);
    clear StatAxis;
    
    for i=1:k
        if (~isnan(Yield(i)) && dataRate(i) > 0 && ~isnan(stats(i,1)))
            StatAxis{i} = InputDataName{i}{2};

            fprintf(fid,'<tr>\n');
            
            % Provide the name for each row of the table ... for some rows link
            % to a plot
            fprintf(fid,'<td>%s</td>\n',InputDataName{i}{2});
            
            if (StartTime(1) >= 0 & EndTime(1) >= 0)
                fprintf(fid,'<td>%.3f</td>\n',Yield(i));
                fprintf(fid,'<td>%.0f</td>\n',1.0/dataRate(i));
            end
            
            for j=1:RowCol(2)
                % Compare this element to the minimum for this statistic. If this
                % element is the "best" (smallest number) then colour the cell
                % green.
                if(stats(i,j) == min(stats(:,j)))
                    fprintf(fid,'<td bgcolor=#7cfc00>%.3f</td>\n',stats(i,j));
                elseif(stats(i,j) == max(stats(:,j)))
                    fprintf(fid,'<td bgcolor=#ff9999>%.3f</td>\n',stats(i,j));
                else
                    fprintf(fid,'<td>%.3f</td>\n',stats(i,j));
                end
            end
            fprintf(fid,'<td>Magnetometer</td>\n');
            
            if (~isnan(ftstats(i,1)))
                for j=1:RowCol(2)
                    % Compare this element to the minimum for this statistic. If this
                    % element is the "best" (smallest number) then colour the cell
                    % green.
                    if(ftstats(i,j) == min(ftstats(:,j)))
                        fprintf(fid,'<td bgcolor=#7cfc00>%.3f</td>\n',ftstats(i,j));
                    elseif(ftstats(i,j) == max(ftstats(:,j)))
                        fprintf(fid,'<td bgcolor=#ff9999>%.3f</td>\n',ftstats(i,j));
                    else
                        fprintf(fid,'<td>%.3f</td>\n',ftstats(i,j));
                    end
                end
            else
                for j=1:RowCol(2)
                    fprintf(fid,'<td>N/A</td>\n');
                end
            end
            fprintf(fid,'<td>Filtered</td>\n');
            
            fprintf(fid,'</tr>\n');
        end
    end
    
    fprintf(fid,'</table>\n');
    
    if (StartTime(1) >= 0 & EndTime(1) >= 0)
        fprintf(fid,['<h5>* - Yield indicates the percent of positions versus expected in '  ...
            'each data set, it is not an indicator of the percentage of positions '  ...
            'for a given position type. e.g. assume the data set is an RTK data set, '  ...
            'it may have 100%% yield but unless the environment is extremely benign '  ...
            'far fewer than 100%% of the data will be RTK, the rest is likely to be '  ...
            'Stinger. To determine the per position type yield review the <i>Primary ' ...
            'Soln</i> for the most frequent position type or for full details ', ...
            'the T012KML summary output. The statistics are ' ...
            'for the entire data set irrespective of position type</h5>\n']);
    end
    fprintf(fid,['<h5>** - Primary Soln is the position type that occurs the most in the ' ...
        'file. For example assume the ' ...
        'expected epochs is 1000, the actual is 800 with ' ...
        '400 RTK fixed, 300 DGNSS and 100 SBAS. The yield would be ' ...
        '80%% (800 out of 1000), the primary solution would be RTK ' ...
        'fixed.</h5>\n']);
    
    [row, col] = size(stats);
    if( (row>1) && (row <= 14) )
        % rotateXLabels() seems to break down if the number of labels is
        % greater than ~14. Don't create the bar charts if we have too much
        % data.
        figure
        ftstatsIndex = find(~isnan(ftstats(:,1)));
        
        bar([stats(:,1:4); ftstats(ftstatsIndex,1:4)]);
        for ss = 1:length(ftstatsIndex)
            StatAxis{row + ss} = ['Filtered ',StatAxis{ftstatsIndex(ss)}];
        end

        A = axis;
        if(s==1)
            title([Section,' Roll Stats']);
            name = [Section,'_Roll_summary.png'];
            A(4) = RPAxisMax;
            ylabel('Error [deg]');
        elseif(s==2)
            title([Section,' Pitch Stats']);
            name = [Section,'_Pitch_summary.png'];
            A(4) = RPAxisMax;
            ylabel('Error [deg]');
        elseif(s==3)
            title([Section,' Heading Stats']);
            name = [Section,'_Heading_summary.png'];
            A(4) = HeadAxisMax;
            ylabel('Error [deg]');
        end
        axis(A);
        set( gca(),'XTickLabel',StatAxis);
        rotateXLabels( gca(), 45);
        legend('50%','68%','95%','99%','location','best');
        grid;
        print('-dpng',name);
        fprintf(fid,['<a href="',name,'"><img src="',name,'" height="600"/></a>']);
    end
    
    fprintf(fid,'<hr/>\n');
end

fprintf(fid,'</body>\n');
fprintf(fid,'</html>\n');

fclose(fid);


