%
% should be called from doDist.m which sets up the names of files etc to
% be processed.
%
% This script reads in a variety of position data and compares to truth.
% It computes the 2D and cross-track error and plots as CDFs.
%

clearvars -global

% T02 position type strings (as a function of "fix type" parameter)
PosStrings = [
    {'Autonomous'}
    {'DGPS'}
    {'SBAS'}
    {'RTK float'}
    {'RTK fixed'}
    {'Wide-area RTK float'}
    {'Wide-area RTK fixed'}
    {'Omnistar'} % Need to check fix_info to determine which OmniSTAR service
    {'CDGPS'}
    {'Autonomous - Kalman'}
    {'DGPS - Kalman'}
    {'SBAS - Kalman'}
    {'CDGPS - Kalman'}
    {'SBAS+ - Kalman'}
    {'CDGPS+ - Kalman'}
    {'RTX (GoldenEye)'}
    {'TLite XPS Hybrid Coarse'}
    {'TLite XPS Hybrid Fine'}
    {'INS Autonomous'}
    {'INS DGNSS'}
    {'INS RTK'}
    {'INS SBAS'}
    {'INS RTX-Carrier'}
    {'INS RTX-Code'}
    {'INS OmniSTAR-Precise'}
    {'INS Dead-Reckoning'}
    {'RTX-Stinger (Code solution)'}
    {'RTX-Fast RAIN'}
    {'xFillX'}
    {'RTX-RangePoint'}
    {'RTX-ViewPoint'}
    {'RTX-FieldPoint'}
    {'RTX-Fast RAM'}
    {'Titan Autonomous'}
    {'Titan DGNSS'}
    {'Titan SBAS'}
    {'SLAS'}
    {'SLAS - Kalman'}
];

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

if(exist('TruthLabel','var') == 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','var') == 0)
  saveFig = 0;
end

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

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

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

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

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

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

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

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

% Closes all the figure after data for each receiver has plotted
if(exist('closeOnComplete','var') == 0)
  closeOnComplete = 0;
end

% remove duplicate records (a known problem for BlueBottle)
if(exist('dup_rec_fix','var') == 0)
  dup_rec_fix = 0;
end

if(exist('ProcessAttitude','var') == 0)
  ProcessAttitude = 1;
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 available set the corresponding index to NaN
if (isempty( find(isnan(POSPac(:,dEastVel)) == 0, 1) ))
  dEastVel = NaN;
end
  
if (isempty( find(isnan(POSPac(:,dNorthVel)) == 0, 1) ))
  dNorthVel = NaN;
end
  
if (isempty( find(isnan(POSPac(:,dUpVel)) == 0, 1) ))
  dUpVel = NaN;
end
  
if (isempty( find(isnan(POSPac(:,dEastSig)) == 0, 1) ))
  dEastSig = NaN;
end
  
if (isempty( find(isnan(POSPac(:,dNorthSig)) == 0, 1) ))
  dNorthSig = NaN;
end
  
if (isempty( find(isnan(POSPac(:,dUpSig)) == 0, 1) ))
  dUpSig = NaN;
end
  
if (isempty( find(isnan(POSPac(:,dRoll)) == 0, 1) ))
  dRoll = NaN;
end
  
if (isempty( find(isnan(POSPac(:,dPitch)) == 0, 1) ))
  dPitch = NaN;
end
  
if (isempty( find(isnan(POSPac(:,dHeading)) == 0, 1) ))
  dHeading = NaN;
end
   
% Store the synced Time/Lat/Long for each data set, this is used in the
% "zoom" plot where based on the zoom of the E/N/U graph we plot a
% sub-set of the data.
clear TimeLatLong;
clear TimeLatLongAll;
% Need to make global so the call back can access the data
global TimeLatLong;
global TimeLatLongAll;

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);

%
% Plot POS Pac estimated errors
%
if(isnan(dEastSig) == 0)
  figure
  plot(PPOut(:,dTime), PPOut(:,dEastSig));
  hold on;
  plot(PPOut(:,dTime), PPOut(:,dNorthSig),'r');
  plot(PPOut(:,dTime), PPOut(:,dUpSig),'g');
  grid
  xlabel('GPS Time [Sec]');
  ylabel('Estimated 1-\sigma Error [m]');
  legend('East','North','Height','Location','Best');
  title([TruthLabel,' Estimated Error']);
  outputFileName = [Section , '_POSPac_Error.png'];
  savePlot(strrep(outputFileName,'\',''),saveFig);
end

%
% 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(strrep(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');
  set(h,'defaultLegendAutoUpdate','off')
  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(strrep(outputFileName,'\',''),saveFig);

  h = figure;
  set(h,'PaperPositionMode','auto');
  set(h,'defaultLegendAutoUpdate','off')
  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(strrep(outputFileName,'\',''),saveFig);
end


figure
plot(PPOut(:,dTime), PPOut(:,dHgt));
xlabel('GPS Time [seconds]');
ylabel('Height [m]');
grid
title([TruthLabel,' Height']);
outputFileName = [Section , ['_',TruthLabel,'_Hgt.png']];
savePlot(strrep(outputFileName,'\',''),saveFig);

%
% 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(strrep(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(strrep(outputFileName,'\',''),saveFig);
end

%
% If we have truth velocity plot it
%
if(isnan(dEastVel) == 0)
  figure
  plot(PPOut(:,dTime), sqrt(PPOut(:,dEastVel).^2 + PPOut(:,dNorthVel).^2));
  xlabel('GPS Time [seconds]');
  ylabel('Horizontal Velocity[m/s]');
  grid
  title([TruthLabel,' Horizontal Velocity']);
  outputFileName = [Section , ['_',TruthLabel,'_HorzVel.png']];
  savePlot(strrep(outputFileName,'\',''),saveFig);


  figure
  % Same plot as above but convert m/s to miles per hour
  plot(PPOut(:,dTime), 2.23694 * sqrt(PPOut(:,dEastVel).^2 + PPOut(:,dNorthVel).^2));
  xlabel('GPS Time [seconds]');
  ylabel('Horizontal Velocity[mph]');
  grid
  title([TruthLabel,' Horizontal Velocity']);
  outputFileName = [Section , ['_',TruthLabel,'_HorzVel_mph.png']];
  savePlot(strrep(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);



% Intialize the stats arrays and cells
dataRate = NaN(1, length(InputDataName));
Yield = NaN(1, length(InputDataName));
EastVel = cell(1, length(InputDataName));
NorthVel = cell(1, length(InputDataName));
HorzVelStats = NaN(length(InputDataName), 5);
HorzVel = cell(1, length(InputDataName));
UpVelStats = NaN(length(InputDataName), 5);
UpVel = cell(1, length(InputDataName));
EastVelStats = NaN(length(InputDataName), 5);
NorthVelStats = NaN(length(InputDataName), 5);
RollErrSort = cell(1, length(InputDataName));
PitchErrSort = cell(1, length(InputDataName));
HeadingErrSave = cell(1, length(InputDataName));
HeadingErrSort = cell(1, length(InputDataName));
RollStats = NaN(length(InputDataName), 5);
PitchStats = NaN(length(InputDataName), 5);
HeadingStats = NaN(length(InputDataName), 5);
PrimarySoln = NaN(1, length(InputDataName));
PrimarySolnPercent = NaN(1, length(InputDataName));
mean_std_rms = NaN(length(InputDataName), 9, 3);
CrossTrackSummary = cell(1, length(InputDataName));
AlongTrackSummary = cell(1, length(InputDataName));
CrossStats = NaN(length(InputDataName), 5);
AlongStats = NaN(length(InputDataName), 5);
VertSummary = cell(1, length(InputDataName));
VertStats = NaN(length(InputDataName), 5);
TwoDStats = NaN(length(InputDataName), 5);
EastStats = NaN(length(InputDataName), 5);
NorthStats = NaN(length(InputDataName), 5);
TwoDSummary = cell(1, length(InputDataName));
Data2D = cell(1, length(InputDataName));
CrossTrackData = cell(1, length(InputDataName));
AlongTrackData = cell(1, length(InputDataName));
DUTData = cell(1, length(InputDataName));
HeightSummary = cell(1, length(InputDataName)); 
HeightData = cell(1, length(InputDataName));
NameArray = cell(1, length(InputDataName));
StatAxis = cell(1, length(InputDataName));

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

    % Handle an optional E/N/U datum shift.
    NorthShift  = 0.0;
    EastShift   = 0.0;
    HeightShift = 0.0;
        
    DataStr = InputDataName{k}{2};
    
    disp('***************************************************************');
    disp('***************************************************************');
    disp(['                ',DataStr]);
    disp('***************************************************************');
    disp('***************************************************************');
      
    if (length(InputDataName{k}) >= 31)
      if(isnan(InputDataName{k}{29}) == 0)
        NorthShift  = InputDataName{k}{29};
      end
      if(isnan(InputDataName{k}{30}) == 0)
        EastShift   = InputDataName{k}{30};
      end
      if(isnan(InputDataName{k}{31}) == 0)
        HeightShift = InputDataName{k}{31};
      end
    end
 
    % 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','var') == 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);
      elseif (strcmpi(InputDataName{k}{4},'VNAV'))
        %1  time            Record time of validity (seconds)           8
        %2  latitude        GPS latitude (degrees)                      8
        %3  longitude       GPS longitude (degrees)                     8
        %4  altitude        GPS ellipsoidal altitude (meters)           8
        %5  velocity[0]     North velocity (meters/second)              8
        %6  velocity[1]     East  velocity (meters/second)              8
        %7  velocity[2]     Down  velocity (meters/second)              8
        %8  attitude[0]     Roll    (degrees)                           8
        %9  attitude[1]     Pitch   (degrees)                           8
        %10 attitude[2]     Heading (degrees)                           8
        %11 wander_angle    Wander angle (degrees)                      8
        %12 acceleration[0] X Body frame acceleration (meters/second2)  8
        %13 acceleration[1] Y Body frame acceleration (meters/second2)  8
        %14 acceleration[2] Z Body frame acceleration (meters/second2)  8
        %15 angular_rate[0] X Body frame angular rate (degrees/second)  8
        %16 angular_rate[1] Y Body frame angular rate (degrees/second)  8
        %17 angular_rate[3] Z Body frame angular rate (degrees/second)  8
        tmp= load_data(FullPathAndName,16);        
        % Round the time tags
        FullData(:,1) = round(tmp(:,1) * 10) / 10;
        % perform the actual decimation to 10Hz by only including data
        % when the time delta changes
        D = find(diff(FullData(:,1)) ~= 0);
        InputData = FullData(D,:);
        clear FullData;

        disp('***************************************************************');
        disp('***************************************************************');
        disp('***** VNAV data is interpolated to the correct time stamp *****');
        disp('***************************************************************');
        disp('***************************************************************');

        % Convert the heading to degrees for the unwrap
        if(strcmpi(InputDataName{k}{9},'RAD')) 
          tmp(:,10) = tmp(:,10) .* 180/pi;
        end

        for i=2:17
          if(i == 10)
            % Make sure we unwrap the heading to avoid spikes
            InputData(:,i) = mod(interp1(tmp(:,1),unwrap(tmp(:,i).*pi/180).*180/pi,InputData(:,1)),360);
          else
            InputData(:,i) = interp1(tmp(:,1),tmp(:,i),InputData(:,1));
          end
        end
          
        if(strcmpi(InputDataName{k}{9},'RAD')) 
          % Convert back to Radians to be consistent with the "RAD"
          % description, later we'll adjust all angles.
          InputData(:,10) = InputData(:,10) .* pi/180;
        end

        save([FullPathAndName,'.dut.mat'],'InputData');
      elseif (strcmpi(InputDataName{k}{4},'XLX'))
        % XLX file format, ignore the header (implemented by Sean Zhi)  
        % SR - can't get this to process XLX files from TrimbleStudio?
        InputData = dlmread(FullPathAndName,'\t',1,0);  
        save([FullPathAndName,'.dut.mat'],'InputData');
      elseif (strcmpi(InputDataName{k}{4},'INSPOS'))
        % WinAstra.inspos, reader from Johannes
        tmp = LoadInsposFile(FullPathAndName);

        i = find(    (isnan(tmp.d(:,4) ) == 0) ...
                  &  (isnan(tmp.d(:,5) ) == 0) ...
                  &  (isnan(tmp.d(:,6) ) == 0) ...
                  &  (isnan(tmp.d(:,7) ) == 0) ...
                  &  (isnan(tmp.d(:,8) ) == 0) ...
                  &  (isnan(tmp.d(:,9) ) == 0) ...
                  &  (isnan(tmp.d(:,10)) == 0) ...
                  &  (isnan(tmp.d(:,11)) == 0) ...
                  &  (isnan(tmp.d(:,12)) == 0) ...
                  &  (isnan(tmp.d(:,26)) == 0) ...
                  &  (isnan(tmp.d(:,27)) == 0) ...
                  &  (isnan(tmp.d(:,28)) == 0) );

        % Time         
        InputData(:,1) = tmp.d(i,2);

        % XYZ -> LLH
        InputData(:,2:4) = xyz2geo(tmp.d(i,4), tmp.d(i,5), tmp.d(i,6));
        InputData(:,2) = InputData(:,2) .* 180 / pi;
        InputData(:,3) = InputData(:,3) .* 180 / pi;
        
        % XYZ Velocity -> ENU Velocity
        InputData(:,5:7) = xyzvel2enu(InputData(:,2), InputData(:,3), ...
                                      tmp.d(i,7), tmp.d(i,8), tmp.d(i,9));
        
        % Roll/Pitch/Heading
        InputData(:,8:10)  = tmp.d(i,10:12);
%         InputData(:,11:13) = tmp.d(i,26:28);
%         Eric - Attitude Sigma
        InputData(:,11:13) = tmp.d(i,38:40);

        % For the XYZ covariance matrix (Upper triangle)
        XYZCov(:,1,1:3) = tmp.d(i,14:16);
        XYZCov(:,2,2:3) = tmp.d(i,17:18);
        XYZCov(:,3,3)   = tmp.d(i,19);
        % Fill in the rest of the matrix
        XYZCov(:,2,1) = XYZCov(:,1,2);
        XYZCov(:,3,1) = XYZCov(:,1,3);
        XYZCov(:,3,2) = XYZCov(:,2,3);

        [ENUCov, semiMinor, semiMajor, orientation] = xyzcov2enucov( InputData(:,2), InputData(:,3), XYZCov);
%         InputData(:,14) = ENUCov(:,1,1);
%         InputData(:,15) = ENUCov(:,2,2);
%         InputData(:,16) = ENUCov(:,3,3);
%         Eric - Pos Sigma
        InputData(:,14) = sqrt(tmp.d(i,20));
        InputData(:,15) = sqrt(tmp.d(i,23));
        InputData(:,16) = sqrt(tmp.d(i,25));


        InputData(:,17) = tmp.d(i,44);
        
%         Eric - Vel Sigma
        InputData(:,18) = sqrt(tmp.d(i,32));
        InputData(:,19) = sqrt(tmp.d(i,35));
        InputData(:,20) = sqrt(tmp.d(i,37));
        
        clear tmp;
        clear i;
          
        save([FullPathAndName,'.dut.mat'],'InputData');
      else
        error('Unknown file format');    
      end
    end
    
    % If the fix type is available..
    if(isnan(InputDataName{k}{19}) == 0)
      % If only considering certain fix types...
      if exist('considerPosType','var')
        % Find unwanted fix types
        j = ~ismember( InputData(:,InputDataName{k}{19}), considerPosType );
        if (length(j))
          % and get rid of them
          InputData(j,:) = [];
        end
      end
    end
    
    % If the alignment status is available...
    if(   (length(InputDataName{k}) >= 37)   ...
       && (isnan(InputDataName{k}{37}) == 0) )
      % Find all solutions with an alignment status > 1
      j = find ( InputData(:,InputDataName{k}{37}) > 1);
      % If the alignment status is > 1 than the INS engine had not 
      % determined precise alignment. We don't want to consider 
      % non-aligned INS solutions in the analysis so delete them.
      if (length(j))
        InputData(j,:) = [];
      end
    end

    % If we have defined start/stop bounds compute the yield and
    % datarate
    clear ActualCount;
    clear ExpectedCount;
    ActualCount = zeros(1,length(StartTime));
    ExpectedCount = zeros(1,length(StartTime));
    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};
        
        i = find( (InputData(:,TimeIndex) >= StartTime(m)) & (InputData(:,TimeIndex) <= EndTime(m)));
        
        % remove duplicate records (a known BlueBottle problem)
        if dup_rec_fix
            TMP02=diff(InputData(i,TimeIndex));
            TMP03=find(TMP02 <= 0);
            if isempty(TMP03) == 0
                InputData(TMP03+1,:)=[];
                i = find( (InputData(:,TimeIndex) >= StartTime(m)) & (InputData(:,TimeIndex) <= EndTime(m)));
            end
        end

        ActualCount(m) = length(i);
        if ActualCount(m) > 0
           % Compute the expected epoch count
           dataRate(k) = min(diff(InputData(i,TimeIndex)));

           if(dataRate(k) < 0.05)
             % It is very slow to process really high rate data so downsample
             % based on the rate
             if(dataRate(k) <= 0.005)
               DownSample = 10;
             elseif(dataRate(k) <= 0.01)
               DownSample = 5;
             elseif(dataRate(k) <= 0.02)
               DownSample = 3;
             else
               DownSample = 1;
             end
           else
             DownSample = 1;
           end

           % Buffalo data had some strange time tags (e.g. 0.91
           % seconds). Find the closest time tag to expected data rates
           rates = [0.005 0.01 0.02 0.05 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}{12} == 1) 
        % Time tag rounding (often have 1ms offset)
        InputData(:,InputDataName{k}{5}) = round(InputData(:,InputDataName{k}{5}).*10)./10;
    end    
    
    % Some consumer devices provide position in mean sea level (geoid
    % height). Adjust by the ellipsoid to geoid difference in Sunnyvale
    % !!!! If this is used elsewhere this needs to be modified !!!!!
    if(strcmpi(InputDataName{k}{10},'MSL')) 
        % MSL offset for Sunnyvale
        InputData(:,InputDataName{k}{8}) = InputData(:,InputDataName{k}{8})-31.894;
    end    
    
    % Novatel data occasionallly returns 0 in Lat/Long via the GGA.
    % Remove the data, not clear whether this is a receiver issue or our
    % scripts
    i = find(   (InputData(:,InputDataName{k}{6}) == 0) ...
              & (InputData(:,InputDataName{k}{7}) == 0) );
    
    InputData(i,InputDataName{k}{5}) = NaN;
    InputData(i,InputDataName{k}{6}) = NaN;
    InputData(i,InputDataName{k}{7}) = NaN;

    % Load common data points with matching time tags
    if( issorted(InputData(:,InputDataName{k}{5})) && issorted(POSPac(:,dTime)) )
        % Use a faster intersect_sorted function if both the time vectors are sorted
        [A, B] = intersect_sorted(InputData(:,InputDataName{k}{5}), POSPac(:,dTime));
    else
        % If somehow they are not sorted, use the default intersect
        [C, A, B] = intersect(InputData(:,InputDataName{k}{5}), POSPac(:,dTime));
    end

    % 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.DUT_lat    =  2;
    key.DUT_lon    =  3;
    key.DUT_hgt    =  4;
    key.truth_lat  =  5;
    key.truth_lon  =  6;
    key.truth_hgt  =  7;
    key.DUT_eVel   =  8;
    key.DUT_nVel   =  9;
    key.DUT_uVel   = 10;
    key.truth_eVel = 11;
    key.truth_nVel = 12;
    key.truth_uVel = 13;
    key.DUT_Roll      = 14;
    key.DUT_Pitch     = 15;
    key.DUT_Heading   = 16;
    key.truth_Roll    = 17;
    key.truth_Pitch   = 18;
    key.truth_Heading = 19;
    key.DUT_sigRoll    = 20;
    key.DUT_sigPitch   = 21;
    key.DUT_sigHeading = 22;
    key.DUT_sigEast  = 23;
    key.DUT_sigNorth = 24;
    key.DUT_sigUp = 25; % Eric - add    
    key.DUT_posType  = 26;

    key.DUT_HDOP = 27;
    key.DUT_VDOP = 28;
    key.DUT_eVelSig = 29;
    key.DUT_nVelSig = 30;
    key.DUT_uVelSig = 31;

    clear Merge;
    Merge = NaN(length(A), length(fieldnames(key)));
    
    % Position, position stats and type
    Merge(:,key.time)    = InputData(A,InputDataName{k}{5});
    Merge(:,key.DUT_lat) = InputData(A,InputDataName{k}{6});
    Merge(:,key.DUT_lon) = InputData(A,InputDataName{k}{7});
    if(strcmpi(InputDataName{k}{9},'RAD'))
        Merge(:,key.DUT_lat) = Merge(:,key.DUT_lat) .* 180/pi;
        Merge(:,key.DUT_lon) = Merge(:,key.DUT_lon) .* 180/pi;
    end
    Merge(:,key.DUT_hgt) = InputData(A,InputDataName{k}{8});
    
    if(isnan(InputDataName{k}{16}) == 0)
        Merge(:,key.DUT_sigEast) = InputData(A,InputDataName{k}{16});
    end
    
    if(isnan(InputDataName{k}{17}) == 0)
        Merge(:,key.DUT_sigNorth) = InputData(A,InputDataName{k}{17});
    end
    
    % Eric - Add sigUp
    if(isnan(InputDataName{k}{18}) == 0)
        Merge(:,key.DUT_sigUp) = InputData(A,InputDataName{k}{18});
    end
    
    if(isnan(InputDataName{k}{19}) == 0)
        Merge(:,key.DUT_posType) = InputData(A,InputDataName{k}{19});
    end
    
    if(   (length(InputDataName{k}) >= 33)   ...
            && (isnan(InputDataName{k}{32}) == 0) ...
            && (isnan(InputDataName{k}{33}) == 0) )
        Merge(:,key.DUT_HDOP) = InputData(A,InputDataName{k}{32});
        Merge(:,key.DUT_VDOP) = InputData(A,InputDataName{k}{33});
    end
    
    if(    (length(InputDataName{k}) >= 37)   ...
            &&  isnan(InputDataName{k}{37}) == 0)
        Merge(:,key.DUT_eVelSig) = InputData(A,InputDataName{k}{37});
    end
    
    
    if(   (length(InputDataName{k}) >= 38)   ...
            && isnan(InputDataName{k}{38}) == 0)
        Merge(:,key.DUT_nVelSig) = InputData(A,InputDataName{k}{38});
    end
    
    if(   (length(InputDataName{k}) >= 39)   ...
            && isnan(InputDataName{k}{39}) == 0)
        Merge(:,key.DUT_uVelSig) = InputData(A,InputDataName{k}{39});
    end
    
    Merge(:,key.truth_lat) = POSPac(B,dLat);
    Merge(:,key.truth_lon) = POSPac(B,dLon);
    Merge(:,key.truth_hgt) = POSPac(B,dHgt);
    
    % Velocity
    if(isnan(InputDataName{k}{14}) == 0)
        Merge(:,key.DUT_eVel) = InputData(A,InputDataName{k}{14});
    end
    
    if(isnan(InputDataName{k}{13}) == 0)
        Merge(:,key.DUT_nVel) = InputData(A,InputDataName{k}{13});
    end
    
    if(isnan(InputDataName{k}{15}) == 0)
        Merge(:,key.DUT_uVel) = InputData(A,InputDataName{k}{15});
    end
    
    if (isnan(dEastVel) == 0)
        Merge(:,key.truth_eVel) = POSPac(B,dEastVel);
    end
    
    if (isnan(dNorthVel) == 0)
        Merge(:,key.truth_nVel) = POSPac(B,dNorthVel);
    end
    
    if (isnan(dUpVel) == 0)
        Merge(:,key.truth_uVel) = POSPac(B,dUpVel);
    end
    
    % Attitude
    if(isnan(InputDataName{k}{20}) == 0)
        Merge(:,key.DUT_Roll) = InputData(A,InputDataName{k}{20});
        
        if(strcmpi(InputDataName{k}{9},'RAD'))
            Merge(:,key.DUT_Roll) = Merge(:,key.DUT_Roll)  .* 180 / pi;
        end
    end
    
    if(isnan(InputDataName{k}{21}) == 0)
        Merge(:,key.DUT_Pitch) = InputData(A,InputDataName{k}{21});
        
        if(strcmpi(InputDataName{k}{9},'RAD'))
            Merge(:,key.DUT_Pitch) = Merge(:,key.DUT_Pitch)  .* 180 / pi;
        end
    end
    
    if(isnan(InputDataName{k}{22}) == 0)
        Merge(:,key.DUT_Heading) = InputData(A,InputDataName{k}{22});
        
        if(strcmpi(InputDataName{k}{9},'RAD'))
            Merge(:,key.DUT_Heading) = Merge(:,key.DUT_Heading)  .* 180 / pi;
        end
    end
    
    if(isnan(dRoll) == 0)
        Merge(:,key.truth_Roll) = POSPac(B,dRoll);
    end
    
    if(isnan(dPitch) == 0)
        Merge(:,key.truth_Pitch)   = POSPac(B,dPitch);
    end
    
    if(isnan(dHeading) == 0)
        Merge(:,key.truth_Heading) = POSPac(B,dHeading);
    end
    
    % Attitude sigma
    if(isnan(InputDataName{k}{26}) == 0)
        Merge(:,key.DUT_sigRoll) = InputData(A,InputDataName{k}{26});
        
        if(strcmpi(InputDataName{k}{9},'RAD'))
            Merge(:,key.DUT_sigRoll) = Merge(:,key.DUT_sigRoll)  .* 180 / pi;
        end
    end
    
    if(isnan(InputDataName{k}{27}) == 0)
        Merge(:,key.DUT_sigPitch) = InputData(A,InputDataName{k}{27});
        
        if(strcmpi(InputDataName{k}{9},'RAD'))
            Merge(:,key.DUT_sigPitch) = Merge(:,key.DUT_sigPitch)  .* 180 / pi;
        end
    end
    
    if(isnan(InputDataName{k}{28}) == 0)
        Merge(:,key.DUT_sigHeading) = InputData(A,InputDataName{k}{28});
        
        if(strcmpi(InputDataName{k}{9},'RAD'))
            Merge(:,key.DUT_sigHeading) = Merge(:,key.DUT_sigHeading)  .* 180 / pi;
        end
    end


    %--------------------------------------------------------------------
    % DUT and Truth data are now merged in a common array, use key.* to
    % access the various elements
    %--------------------------------------------------------------------
    
    if(exist('dataRate','var'))
      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(   (length(InputDataName{k}) >= 33)   ...
       && (isnan(InputDataName{k}{32}) == 0) ...
       && (isnan(InputDataName{k}{33}) == 0) )
      figure;
      % In 2D mode we don't have VDOP, it gets saved as zero filter out
      % to prevent 0 VDOP in the graphs
      z = find(Merge(:,key.DUT_VDOP) < 0.5);
      VDOP = Merge(:,key.DUT_VDOP);
      VDOP(z) = NaN; % force a pen-up
      plot(Merge(:,key.time), Merge(:,key.DUT_HDOP),'b');
      hold on;
      plot(Merge(:,key.time), VDOP,'r');
      Ax = axis;
      Ax(3) = 0;
      Ax(4) = 5;
      axis(Ax);
      grid;
      title(['DOP - ',DataStr,' - ',Section]);
      xlabel('Time [GPS Seconds]');
      ylabel('DOP');
      legend(['HDOP \mu=',num2str(mean(Merge(GoodTimeTag,key.DUT_HDOP)),'%.1f')], ...
             ['VDOP \mu=',num2str(mean(Merge(GoodTimeTag,key.DUT_VDOP)),'%.1f')], ...
             'location','best');
      outputFileName = [Section '_' DataStr '_DOP.png'];
      savePlot(strrep(outputFileName,'\',''),saveFig);
    end
    
    % velocity plot
    if(  (isnan(InputDataName{k}{14}) == 0) ...
           && (isnan(dEastVel) == 0) )
        figure
        plot(Merge(:,key.time), Merge(:,key.truth_eVel) - Merge(:,key.DUT_eVel));
        hold on;
        plot(Merge(:,key.time), Merge(:,key.truth_nVel) - Merge(:,key.DUT_nVel),'r');
        plot(Merge(:,key.time), Merge(:,key.truth_uVel) - Merge(:,key.DUT_uVel),'g');
        Ax = axis;
        Ax(3) = -VelAxisMax;
        Ax(4) =  VelAxisMax;
        axis(Ax)
        grid
        xlabel('GPS Time [Sec]');
        ylabel('Velocity Error [m/s]');
        legend('East','North','Height','Location','Best');
        title(['Velocity Error - ',DataStr,' - ',Section]);
        outputFileName = [Section '_' DataStr '_Vel.png'];
        savePlot(strrep(outputFileName,'\',''),saveFig);
        
        if(    (length(InputDataName{k}) >= 37)   ...
           &&  isnan(InputDataName{k}{37}) == 0)
          figure
          plot(Merge(:,key.time), Merge(:,key.truth_eVel) - Merge(:,key.DUT_eVel));
          hold on
          plot(Merge(:,key.time), Merge(:,key.DUT_eVelSig),'r');
          plot(Merge(:,key.time), -Merge(:,key.DUT_eVelSig),'r');
          legend('East Vel Error','+East Sigma Est','-East Sigma Est','Location','Best');
          title(['East Velocity Error vs Sigma - ',DataStr,' - ',Section]);
          xlabel('GPS Time [Sec]');
          ylabel('Velocity Error [m/s]');
          outputFileName = [Section '_' DataStr '_EVelVsSig.png'];
          savePlot(strrep(outputFileName,'\',''),0);
        end


        if(    (length(InputDataName{k}) >= 38)   ...
           &&  isnan(InputDataName{k}{38}) == 0)
          figure
          plot(Merge(:,key.time), Merge(:,key.truth_nVel) - Merge(:,key.DUT_nVel));
          hold on
          plot(Merge(:,key.time), Merge(:,key.DUT_nVelSig),'r');
          plot(Merge(:,key.time), -Merge(:,key.DUT_nVelSig),'r');
          legend('North Vel Error','+North Sigma Est','-North Sigma Est','Location','Best');
          title(['North Velocity Error vs Sigma - ',DataStr,' - ',Section]);
          xlabel('GPS Time [Sec]');
          ylabel('Velocity Error [m/s]');
          outputFileName = [Section '_' DataStr '_NVelVsSig.png'];
          savePlot(strrep(outputFileName,'\',''),0);
        end


        if(    (length(InputDataName{k}) >= 39)   ...
           &&  isnan(InputDataName{k}{39}) == 0)
          figure
          plot(Merge(:,key.time), Merge(:,key.truth_uVel) - Merge(:,key.DUT_uVel));
          hold on
          plot(Merge(:,key.time), Merge(:,key.DUT_uVelSig),'r');
          plot(Merge(:,key.time), -Merge(:,key.DUT_uVelSig),'r');
          legend('Up Vel Error','+Up Sigma Est','-Up Sigma Est','Location','Best');
          title(['Up Velocity Error vs Sigma - ',DataStr,' - ',Section]);
          xlabel('GPS Time [Sec]');
          ylabel('Velocity Error [m/s]');
          outputFileName = [Section '_' DataStr '_UVelVsSig.png'];
          savePlot(strrep(outputFileName,'\',''),0);
        end

        i = find(isnan(Merge(:,key.DUT_eVel)) == 0);
        VEast = Merge(i,key.truth_eVel) - Merge(i,key.DUT_eVel);
        disp([DataStr, ' RMS(E_Vel) :', num2str(rms(VEast)), ' m/s']);
        disp([DataStr, ' max(abs(E_Vel)) :', num2str(max(abs(VEast))), ' m/s']);
        EastVel{k}  = sort(abs(VEast))';
        
        i = find(isnan(Merge(:,key.DUT_nVel)) == 0);
        VNorth = Merge(i,key.truth_nVel) - Merge(i,key.DUT_nVel);
        disp([DataStr, ' RMS(N_Vel) :', num2str(rms(VNorth)), ' m/s']);
        disp([DataStr, ' max(abs(N_Vel)) :', num2str(max(abs(VNorth))), ' m/s']);
        NorthVel{k} = sort(abs(VNorth))';
        
        HzVel = sort(sqrt(VEast.^2 + VNorth.^2))';
        HorzVelStats(k,1) = HzVel(floor(0.50*length(HzVel)));
        HorzVelStats(k,2) = HzVel(floor(0.68*length(HzVel)));
        HorzVelStats(k,3) = HzVel(floor(0.95*length(HzVel)));
        HorzVelStats(k,4) = HzVel(floor(0.99*length(HzVel)));
        HorzVelStats(k,5) = HzVel(end);
        HorzVel{k} = HzVel;
        
        i = find(isnan(Merge(:,key.DUT_uVel)) == 0);
        VUp = Merge(i,key.truth_uVel) - Merge(i,key.DUT_uVel);
        if(sum(isfinite(VUp)))
          disp([DataStr, ' RMS(U_Vel) :', num2str(rms(VUp)), ' m/s']);
          disp([DataStr, ' max(abs(U_Vel)) :', num2str(max(abs(VUp))), ' m/s']);
          UpV    = sort(abs(VUp))';
          UpVelStats(k,1) = UpV(floor(0.50*length(UpV)));
          UpVelStats(k,2) = UpV(floor(0.68*length(UpV)));
          UpVelStats(k,3) = UpV(floor(0.95*length(UpV)));
          UpVelStats(k,4) = UpV(floor(0.99*length(UpV)));
          UpVelStats(k,5) = UpV(end);
          UpVel{k}    = UpV;
        end

        % Eric - Add Velocity Stats East and North
        EastV    = sort(abs(VEast))';
        EastVelStats(k,1) = EastV(floor(0.50*length(EastV)));
        EastVelStats(k,2) = EastV(floor(0.68*length(EastV)));
        EastVelStats(k,3) = EastV(floor(0.95*length(EastV)));
        EastVelStats(k,4) = EastV(floor(0.99*length(EastV)));
        EastVelStats(k,5) = EastV(end);
        
        NorthV    = sort(abs(VNorth))';
        NorthVelStats(k,1) = NorthV(floor(0.50*length(NorthV)));
        NorthVelStats(k,2) = NorthV(floor(0.68*length(NorthV)));
        NorthVelStats(k,3) = NorthV(floor(0.95*length(NorthV)));
        NorthVelStats(k,4) = NorthV(floor(0.99*length(NorthV)));
        NorthVelStats(k,5) = NorthV(end);

    end;

%     % East velocity versus truth
%     if(  (isnan(InputDataName{k}{14}) == 0) ...
%        & (isnan(dEastVel) == 0) )
%       
%       figure
%       plot(Merge(:,key.time), Merge(:,key.truth_eVel) - Merge(:,key.DUT_eVel));
%       title(['East Velocity Error - ',DataStr,' - ',Section]);
%       xlabel('Time [GPS Seconds]');
%       ylabel('Velocity Error [m/s]');
%       Ax = axis;
%       Ax(3) = -1;
%       Ax(4) =  1;
%       axis(Ax)
%       grid
%       outputFileName = [Section '_' DataStr '_East_Vel.png'];
%       savePlot(strrep(outputFileName,'\',''),saveFig);
% 
%       i = find(isnan(Merge(:,key.DUT_eVel)) == 0);
%       VEast = Merge(i,key.truth_eVel) - Merge(i,key.DUT_eVel);
%       disp([DataStr, ' RMS(E_Vel) :', num2str(rms(VEast)), ' m/s']);
%       disp([DataStr, ' max(abs(E_Vel)) :', num2str(max(abs(VEast))), ' m/s']);
%       EastVel{k}  = sort(abs(VEast))';
%     end
% 
%     % North velocity versus truth
%     if(  (isnan(InputDataName{k}{13}) == 0) ...
%        & (isnan(dNorthVel) == 0) )
% 
%       figure
%       plot(Merge(:,key.time), Merge(:,key.truth_nVel) - Merge(:,key.DUT_nVel));
%       title(['North Velocity Error - ',DataStr,' - ',Section]);
%       xlabel('Time [GPS Seconds]');
%       ylabel('Velocity Error [m/s]');
%       Ax = axis;
%       Ax(3) = -1;
%       Ax(4) =  1;
%       axis(Ax)
%       grid
%       outputFileName = [Section '_' DataStr '_North_Vel.png'];
%       savePlot(strrep(outputFileName,'\',''),saveFig);
%       
%       i = find(isnan(Merge(:,key.DUT_nVel)) == 0);
%       VNorth = Merge(i,key.truth_nVel) - Merge(i,key.DUT_nVel);
%       disp([DataStr, ' RMS(N_Vel) :', num2str(rms(VNorth)), ' m/s']);
%       disp([DataStr, ' max(abs(N_Vel)) :', num2str(max(abs(VNorth))), ' m/s']);
%       NorthVel{k} = sort(abs(VNorth))';
%     end
% 
%     % Up velocity versus truth
%     if(  (isnan(InputDataName{k}{15}) == 0) ...
%        & (isnan(dUpVel) == 0) )
% 
%       figure
%       plot(Merge(:,key.time), Merge(:,key.truth_uVel) - Merge(:,key.DUT_uVel));
%       title(['Hgt Velocity Error - ',DataStr,' - ',Section]);
%       xlabel('Time [GPS Seconds]');
%       ylabel('Velocity Error [m/s]');
%       Ax = axis;
%       Ax(3) = -1;
%       Ax(4) =  1;
%       axis(Ax)
%       grid
%       outputFileName = [Section '_' DataStr '_Hgt_Vel.png'];
%       savePlot(strrep(outputFileName,'\',''),saveFig);
%       
%       i = find(isnan(Merge(:,key.DUT_uVel)) == 0);
%       VUp = Merge(i,key.truth_uVel) - Merge(i,key.DUT_uVel);
%       disp([DataStr, ' RMS(U_Vel) :', num2str(rms(VUp)), ' m/s']);
%       disp([DataStr, ' max(abs(U_Vel)) :', num2str(max(abs(VUp))), ' m/s']);
%       UpVel{k}    = sort(abs(VUp))';
%     end

    %
    % If we have truth and estimated attitude...
    if(  (isnan(InputDataName{k}{20}) == 0) ...
      && (isnan(InputDataName{k}{21}) == 0) ...
      && (isnan(InputDataName{k}{22}) == 0) ...
      && (isnan(dRoll) == 0) )

      % Check for corrections...
      if (isnan(InputDataName{k}{23}) == 0)
        RollCorr = InputDataName{k}{23};
      else
        RollCorr = 0;
      end
      if (isnan(InputDataName{k}{24}) == 0)
        PitchCorr = InputDataName{k}{24};
      else
        PitchCorr = 0;
      end
      if (isnan(InputDataName{k}{25}) == 0)
        HeadingCorr = InputDataName{k}{25};
      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(:,key.DUT_Roll),     ... % Roll
                                                                   Merge(:,key.DUT_Pitch),    ... % Pitch
                                                                   Merge(:,key.DUT_Heading),  ... % Heading
                                                                   RollCorr,  ...                 % Xrotation Correction
                                                                   PitchCorr,  ...                % Yrotation Correction
                                                                   HeadingCorr);                  % Zrotation Correction
      else
        RollHat = Merge(:,key.DUT_Roll);
        PitchHat = Merge(:,key.DUT_Pitch);
        HeadingHat = Merge(:,key.DUT_Heading);
      end


      % now compute the roll/pitch/heading error
    %  i = find(isnan(RollHat) == 0);    
      [RollErr, PitchErr, HeadingErr] = compareRollPitchHeading( Merge(:,key.truth_Roll), ...
                                                                 Merge(:,key.truth_Pitch), ...
                                                                 Merge(:,key.truth_Heading), ...
                                                                 RollHat(:)', ...
                                                                 PitchHat(:)', ...
                                                                 HeadingHat(:)');


      i = find(isnan(RollErr) == 0);    
      if( (isempty(i) == 0) & (ProcessAttitude == 1) )
          
          RollErr = RollErr';
          PitchErr = PitchErr';
          HeadingErr = HeadingErr';
          
          % Roll error w.r.t. POSPac
          figure
          plot(Merge(:,key.time), RollErr);
          hold on;
          if(isnan(InputDataName{k}{26}) == 0)
              % +/- 1-sigma estimates
              plot(Merge(:,key.time), RollErr + Merge(:,key.DUT_sigRoll),'m');
              plot(Merge(:,key.time), RollErr - Merge(:,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(strrep(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
          
          % Pitch error w.r.t. POSPac
          figure
          plot(Merge(:,key.time), PitchErr);
          hold on;
          if(isnan(InputDataName{k}{27}) == 0)
              % +/- 1-sigma estimates
              plot(Merge(:,key.time), PitchErr + Merge(:,key.DUT_sigPitch),'m');
              plot(Merge(:,key.time), PitchErr - Merge(:,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(strrep(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
          
          % Heading error w.r.t. POSPac
          j = find(HeadingErr > 180);
          HeadingErr(j) = HeadingErr(j) - 360;
          figure
          plot(Merge(:,key.time), HeadingErr);
          hold on;
          if(isnan(InputDataName{k}{28}) == 0)
              % +/- 1-sigma estimates
              plot(Merge(:,key.time), HeadingErr + Merge(:,key.DUT_sigHeading),'m');
              plot(Merge(:,key.time), HeadingErr - Merge(:,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(strrep(outputFileName,'\',''),saveFig);
          
          disp(['Heading 1-sigma= ',num2str(std(HeadingErr(i)),'%.3f'),'deg']);
          disp(['Heading mean   = ',num2str(mean(HeadingErr(i)),'%.3f'),'deg']);
          HeadingErrSave{k} = HeadingErr(i);
          HeadingErrSort{k} = sort(abs(HeadingErr(i))); % Sort for the CDF
          
          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);
          
          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);
          
          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);
          
          disp(['Median Roll, Pitch, Heading (deg) = ',num2str(median(RollErr(i)),'%.3f'),', ',num2str(median(PitchErr(i)),'%.3f'),', ',num2str(median(HeadingErr(i)),'%.3f')]);
      end;
    end
   
    % Calculate position errors
    ENU = LLH2ENU(Merge(:,key.DUT_lat) .* pi/180,   ...
                  Merge(:,key.DUT_lon) .* pi/180,   ...
                  Merge(:,key.DUT_hgt),             ...
                  Merge(:,key.truth_lat) .* pi/180, ...
                  Merge(:,key.truth_lon) .* pi/180, ...
                  Merge(:,key.truth_hgt));
    
    % Shift to the truth datum
    ENU(GoodTimeTag,1) = ENU(GoodTimeTag,1) - EastShift;
    ENU(GoodTimeTag,2) = ENU(GoodTimeTag,2) - NorthShift;
    ENU(GoodTimeTag,3) = ENU(GoodTimeTag,3) - HeightShift;
    
    % If heading is available...
    if(isnan(dHeading) == 0)
        
      % If roll and pitch are also available...
      if (isnan(dRoll) == 0 && isnan(dPitch) == 0)          
        % Compute the full vehicle to navigation frame rotation matrix
        %
        %   cp*sy   cr*cy + sr*sp*sy  -sr*cy + cr*sp*sy
        %   cp*cy  -cr*sy + sr*sp*cy   sr*sy + cr*sp*cy 
        %      sp             -sr*cp             -cr*cp
        %          
        roll = Merge(:,key.truth_Roll) * pi / 180;
        pitch = Merge(:,key.truth_Pitch) * pi / 180;
        yaw = Merge(:,key.truth_Heading) * pi / 180;
        cr = cos(roll);
        sr = sin(roll);
        cp = cos(pitch);
        sp = sin(pitch);
        cy = cos(yaw);
        sy = sin(yaw);
        V2N11 = cp.*sy ;
        V2N12 = cr.*cy + sr.*sp.*sy;
        V2N13 = -sr.*cy + cr.*sp.*sy;
        V2N21 = cp.*cy;
        V2N22 = -cr.*sy + sr.*sp.*cy;
        V2N23 = sr.*sy + cr.*sp.*cy;
        V2N31 = sp;
        V2N32 = -sr.*cp;
        V2N33 = -cr.*cp;
        % Get angular rate from truth attitude
        OmegaX = sin(diff(roll));
        OmegaY = sin(diff(pitch));
        OmegaZ = sin(diff(yaw));
      else
        % Compute the partial body frame to ENU rotation matrix
        %
        %   sy   cy   0
        %   cy  -sy   0
        %    0    0  -1
        %          
        yaw = Merge(:,key.truth_Heading) * pi / 180;
        cy = cos(yaw);
        sy = sin(yaw);
        z = zeros(size(yaw));
        o = ones(size(yaw));
        V2N11 = sy;
        V2N12 = cy;
        V2N13 = z;
        V2N21 = cy;
        V2N22 = -sy;
        V2N23 = z;
        V2N31 = z;
        V2N32 = z;
        V2N33 = -o;
        % Get angular rate from truth attitude
        OmegaX = zeros(length(yaw) - 1,1);
        OmegaY = zeros(length(yaw) - 1,1);
        OmegaZ = sin(diff(yaw));
      end

      % Get vehicle frame XYZ (along track, cross track and vectical)
      % corrections if available. These are corrections to the position 
      % of the sensor in the vehicle frame which is centered at the 
      % primary GNSS antenna
      % If corrections for all three axis are available... 
      if(   (length(InputDataName{k}) >= 36)   ...
         && (isnan(InputDataName{k}{34}) == 0) ...
         && (isnan(InputDataName{k}{35}) == 0) ...
         && (isnan(InputDataName{k}{36}) == 0) )
        Xv = InputDataName{k}{34};
        Yv = InputDataName{k}{35};
        Zv = InputDataName{k}{36};
      else
        Xv = 0;
        Yv = 0;
        Zv = 0;    
        % Special case of along track only corrections
        if ( isnan(InputDataName{k}{11}) == 0 )
          Xv = InputDataName{k}{11};
        end
      end
      
      % Apply the vehicle frame corrections to the ENU errors
      ENU(GoodTimeTag,1) = ENU(GoodTimeTag,1) - (Xv .* V2N11(GoodTimeTag) ...
                                               + Yv .* V2N12(GoodTimeTag) ...
                                               + Zv .* V2N13(GoodTimeTag));
      ENU(GoodTimeTag,2) = ENU(GoodTimeTag,2) - (Xv .* V2N21(GoodTimeTag) ...
                                               + Yv .* V2N22(GoodTimeTag) ...
                                               + Zv .* V2N23(GoodTimeTag));
      ENU(GoodTimeTag,3) = ENU(GoodTimeTag,3) - (Xv .* V2N31(GoodTimeTag) ...
                                               + Yv .* V2N32(GoodTimeTag) ...
                                               + Zv .* V2N33(GoodTimeTag));

      % If using internal antenna which is different from the truth
      % The velocity of the receiver need to be corrected for lever arm to get
      % the velocity at the external antenna position (where truth is generated)
      % Lever arm correction for velocity
      if (Xv ~= 0 || Yv ~= 0 || Zv ~= 0)
          deltaT = diff(Merge(:,key.time));
          vx_corr = (-OmegaZ .* Yv + OmegaY .* Zv)./deltaT;
          vy_corr = (OmegaZ .* Xv - OmegaX .* Zv)./deltaT;
          vZ_corr = (-OmegaY .* Xv + OmegaX .* Yv)./deltaT;
          % Estimate the velocity at external antenna
          Merge(2:end,key.DUT_eVel) = Merge(2:end,key.DUT_eVel) ...
                                    - (vx_corr .* V2N11(2:end) ...
                                    + vy_corr .* V2N12(2:end) ...
                                    + vZ_corr .* V2N13(2:end));
          Merge(2:end,key.DUT_nVel) = Merge(2:end,key.DUT_nVel) ...
                                    - (vx_corr .* V2N21(2:end) ...
                                    + vy_corr .* V2N22(2:end) ...
                                    + vZ_corr .* V2N23(2:end));
          Merge(2:end,key.DUT_uVel) = Merge(2:end,key.DUT_uVel) ...
                                    - (vx_corr .* V2N31(2:end) ...
                                    + vy_corr .* V2N32(2:end) ...
                                    + vZ_corr .* V2N33(2:end));
          
          % plot corrected velocity error
          if(  (isnan(InputDataName{k}{14}) == 0) ...
                  && (isnan(dEastVel) == 0) )
              figure
              plot(Merge(:,key.time), Merge(:,key.truth_eVel) - Merge(:,key.DUT_eVel));
              hold on;
              plot(Merge(:,key.time), Merge(:,key.truth_nVel) - Merge(:,key.DUT_nVel),'r');
              plot(Merge(:,key.time), Merge(:,key.truth_uVel) - Merge(:,key.DUT_uVel),'g');
              Ax = axis;
              Ax(3) = -VelAxisMax;
              Ax(4) =  VelAxisMax;
              axis(Ax)
              grid
              xlabel('GPS Time [Sec]');
              ylabel('Velocity Error [m/s]');
              legend('East','North','Height','Location','Best');
              title(['Velocity Error - ',DataStr,' - ',Section]);
              outputFileName = [Section '_' DataStr '_Vel.png'];
              savePlot(strrep(outputFileName,'\',''),saveFig);
              
              if(    (length(InputDataName{k}) >= 37)   ...
                      &&  isnan(InputDataName{k}{37}) == 0)
                  figure
                  plot(Merge(:,key.time), Merge(:,key.truth_eVel) - Merge(:,key.DUT_eVel));
                  hold on
                  plot(Merge(:,key.time), Merge(:,key.DUT_eVelSig),'r');
                  plot(Merge(:,key.time), -Merge(:,key.DUT_eVelSig),'r');
                  legend('East Vel Error','+East Sigma Est','-East Sigma Est','Location','Best');
                  title(['East Velocity Error vs Sigma - ',DataStr,' - ',Section]);
                  xlabel('GPS Time [Sec]');
                  ylabel('Velocity Error [m/s]');
                  outputFileName = [Section '_' DataStr '_EVelVsSig.png'];
                  savePlot(strrep(outputFileName,'\',''),0);
              end
              
              
              if(    (length(InputDataName{k}) >= 38)   ...
                      &&  isnan(InputDataName{k}{38}) == 0)
                  figure
                  plot(Merge(:,key.time), Merge(:,key.truth_nVel) - Merge(:,key.DUT_nVel));
                  hold on
                  plot(Merge(:,key.time), Merge(:,key.DUT_nVelSig),'r');
                  plot(Merge(:,key.time), -Merge(:,key.DUT_nVelSig),'r');
                  legend('North Vel Error','+North Sigma Est','-North Sigma Est','Location','Best');
                  title(['North Velocity Error vs Sigma - ',DataStr,' - ',Section]);
                  xlabel('GPS Time [Sec]');
                  ylabel('Velocity Error [m/s]');
                  outputFileName = [Section '_' DataStr '_NVelVsSig.png'];
                  savePlot(strrep(outputFileName,'\',''),0);
              end
              
              
              if(    (length(InputDataName{k}) >= 39)   ...
                      &&  isnan(InputDataName{k}{39}) == 0)
                  figure
                  plot(Merge(:,key.time), Merge(:,key.truth_uVel) - Merge(:,key.DUT_uVel));
                  hold on
                  plot(Merge(:,key.time), Merge(:,key.DUT_uVelSig),'r');
                  plot(Merge(:,key.time), -Merge(:,key.DUT_uVelSig),'r');
                  legend('Up Vel Error','+Up Sigma Est','-Up Sigma Est','Location','Best');
                  title(['Up Velocity Error vs Sigma - ',DataStr,' - ',Section]);
                  xlabel('GPS Time [Sec]');
                  ylabel('Velocity Error [m/s]');
                  outputFileName = [Section '_' DataStr '_UVelVsSig.png'];
                  savePlot(strrep(outputFileName,'\',''),0);
              end
              
              i = find(isnan(Merge(:,key.DUT_eVel)) == 0);
              VEast = Merge(i,key.truth_eVel) - Merge(i,key.DUT_eVel);
              disp([DataStr, ' RMS(E_Vel) :', num2str(rms(VEast)), ' m/s']);
              disp([DataStr, ' max(abs(E_Vel)) :', num2str(max(abs(VEast))), ' m/s']);
              EastVel{k}  = sort(abs(VEast))';
              
              i = find(isnan(Merge(:,key.DUT_nVel)) == 0);
              VNorth = Merge(i,key.truth_nVel) - Merge(i,key.DUT_nVel);
              disp([DataStr, ' RMS(N_Vel) :', num2str(rms(VNorth)), ' m/s']);
              disp([DataStr, ' max(abs(N_Vel)) :', num2str(max(abs(VNorth))), ' m/s']);
              NorthVel{k} = sort(abs(VNorth))';
              
              HzVel = sort(sqrt(VEast.^2 + VNorth.^2))';
              HorzVelStats(k,1) = HzVel(floor(0.50*length(HzVel)));
              HorzVelStats(k,2) = HzVel(floor(0.68*length(HzVel)));
              HorzVelStats(k,3) = HzVel(floor(0.95*length(HzVel)));
              HorzVelStats(k,4) = HzVel(floor(0.99*length(HzVel)));
              HorzVelStats(k,5) = HzVel(end);
              HorzVel{k} = HzVel;
              
              i = find(isnan(Merge(:,key.DUT_uVel)) == 0);
              VUp = Merge(i,key.truth_uVel) - Merge(i,key.DUT_uVel);
              if(sum(isfinite(VUp)))
                  disp([DataStr, ' RMS(U_Vel) :', num2str(rms(VUp)), ' m/s']);
                  disp([DataStr, ' max(abs(U_Vel)) :', num2str(max(abs(VUp))), ' m/s']);
                  UpV    = sort(abs(VUp))';
                  UpVelStats(k,1) = UpV(floor(0.50*length(UpV)));
                  UpVelStats(k,2) = UpV(floor(0.68*length(UpV)));
                  UpVelStats(k,3) = UpV(floor(0.95*length(UpV)));
                  UpVelStats(k,4) = UpV(floor(0.99*length(UpV)));
                  UpVelStats(k,5) = UpV(end);
                  UpVel{k}    = UpV;
              end
              
              % Eric - Add Velocity Stats East and North
              EastV    = sort(abs(VEast))';
              EastVelStats(k,1) = EastV(floor(0.50*length(EastV)));
              EastVelStats(k,2) = EastV(floor(0.68*length(EastV)));
              EastVelStats(k,3) = EastV(floor(0.95*length(EastV)));
              EastVelStats(k,4) = EastV(floor(0.99*length(EastV)));
              EastVelStats(k,5) = EastV(end);
              
              NorthV    = sort(abs(VNorth))';
              NorthVelStats(k,1) = NorthV(floor(0.50*length(NorthV)));
              NorthVelStats(k,2) = NorthV(floor(0.68*length(NorthV)));
              NorthVelStats(k,3) = NorthV(floor(0.95*length(NorthV)));
              NorthVelStats(k,4) = NorthV(floor(0.99*length(NorthV)));
              NorthVelStats(k,5) = NorthV(end);
              
          end
      end
    end

    % Plot the horizontal error versus sigma
    if(   (isnan(InputDataName{k}{16}) == 0) ...
       && (isnan(InputDataName{k}{17}) == 0) )
      % East/North Sigma
      figure
      HorzErr = sqrt(ENU(:,1).^2 + ENU(:,2).^2);
      plot(Merge(:,key.time), HorzErr);
      hold on;
      
      EN_Sigma = sqrt(  Merge(:,key.DUT_sigEast).^2 ...
                      + Merge(:,key.DUT_sigNorth).^2);
      
      i = find(HorzErr <= EN_Sigma);
      disp(['Horizontal Error less than Sigma Est = ',num2str(100*length(i)/length(HorzErr(GoodTimeTag))),'%']);
      disp(['Mean Horizontal Error = ',num2str(mean(HorzErr(GoodTimeTag))),'m']);
      disp(['Mean Horizontal Error Est = ',num2str(mean(EN_Sigma(GoodTimeTag))),'m']);

      plot(Merge(:,key.time), EN_Sigma,'g');
      grid
      Ax = axis;
      Ax(3) = 0;
      Ax(4) = 2;
      axis(Ax);
      title([DataStr,' - ',Section]);
      legend('Error','Estimated Error');
      ylabel('Horizontal Error [m]');
      xlabel('Time [GPS Seconds]');
      outputFileName = [Section '_' DataStr '_ENErrVsEst.png'];
      savePlot(strrep(outputFileName,'\',''),saveFig);
    end
    
    % To be expanded - analysis based on position type
    if(isnan(InputDataName{k}{19}) == 0)
      % T02 data type
      PosType = Merge(:,key.DUT_posType);

      MaxPos = max(PosType);
      MinPos = min(PosType);
      
      PrimarySoln(k) = 0;
      PrimarySolnPercent(k) = 0;

      disp('****************** Pos Type ******************************');
      % Display percentages of each position type and cache the
      % percentage for the position solution that occurs the mode
      PosType_length=length(find(isnan(PosType)==0));
      TypeIndex = 1;
      for type=MinPos:MaxPos
        i = find(PosType == type);
        if(~isempty(i))
          Percent = 100*length(i)/PosType_length;
          disp(['Type:',num2str(type),' - ',PosStrings{type+1},': ',num2str(Percent,'%.3f'),'%%']);
          % Save as a fraction. Before we report in the web page we'll
          % adjust based on the yield and report the "no position" percent
          % as well.
          PosTypeArray{k}(TypeIndex,1) = type; % The position type
          PosTypeArray{k}(TypeIndex,2) = Percent/100;
          TypeIndex = TypeIndex+1;

          if(Percent > PrimarySolnPercent(k))
            PrimarySoln(k) = type;
            PrimarySolnPercent(k) = Percent;
          end
        end
      end
      disp('**********************************************************');
    end

    TimeLatLong{k} = [Merge(:,key.time),Merge(:,key.DUT_lon),Merge(:,key.DUT_lat)];
      
    MapFig = figure;
    set(MapFig,'PaperPositionMode','auto');
    set(MapFig,'defaultLegendAutoUpdate','off')
    h = plot(TimeLatLong{k}(:,2), TimeLatLong{k}(:,3),'r','linewidth',3);
    xlabel('Longitude [degree]');
    ylabel('Latitude [degree]');
    title([DataStr,' - ',Section]);
    if(plotMap)
      disp('Please wait loading google map ...');
      plot_google_map('maptype','roadmap');

      outputFileName = [Section '_' DataStr '_EN_MapView.png'];
      savePlot(strrep(outputFileName,'\',''),saveFig);

      outputFileName = [Section '_' DataStr '_EN_MapView.html'];
      plotOpenStreetMap(outputFileName,TimeLatLong{k}(:,3), TimeLatLong{k}(:,2));
    end

    if( (zoomLink ==1) && (length(InputDataName) <= 17) )
      % If zoomLink is set to 1 we're going to link two graphs. Zooming 
      % in/out of the E/N/U graph will cause the data to change in a Lat/Long 
      % graph we are about to draw. The Lat/Long graph has optionally a 
      % google maps background. This gives spatial situational data to 
      % errors. For example if you spot an outlier in the E/N/U graph 
      % zoom into it and only the data in the x-axis range (time) of 
      % that sub-plot will be shown on the Lat/Long graph. To do this we 
      % use the matlab "linkdata" feature together with the ActionPostCallback
      % connected to the zoom of the E/N/U graph.

      % The callback copies a subset of the full data set into
      % TimeLatLong, we start of with all the data.

      TimeLatLongAll{k} = TimeLatLong{k};

      % Associate the X/Y data series to variables - there has to be a
      % better way to do this. We want to associate a different set of
      % data with each plot. If you use TimeLatLong{k}(:,2) when k
      % changes the data changes so each plot ends up with the same
      % data. Only by explicitly setting the element number in {} do we
      % get the correct data associated with each plot!
      if(k==1)
        set(h,'XDataSource','TimeLatLong{1}(:,2)'); set(h,'YDataSource','TimeLatLong{1}(:,3)');
      elseif(k==2)
        set(h,'XDataSource','TimeLatLong{2}(:,2)'); set(h,'YDataSource','TimeLatLong{2}(:,3)');
      elseif(k==3)
        set(h,'XDataSource','TimeLatLong{3}(:,2)'); set(h,'YDataSource','TimeLatLong{3}(:,3)');
      elseif(k==4)
        set(h,'XDataSource','TimeLatLong{4}(:,2)'); set(h,'YDataSource','TimeLatLong{4}(:,3)');
      elseif(k==5)
        set(h,'XDataSource','TimeLatLong{5}(:,2)'); set(h,'YDataSource','TimeLatLong{5}(:,3)');
      elseif(k==6)
        set(h,'XDataSource','TimeLatLong{6}(:,2)'); set(h,'YDataSource','TimeLatLong{6}(:,3)');
      elseif(k==7)
        set(h,'XDataSource','TimeLatLong{7}(:,2)'); set(h,'YDataSource','TimeLatLong{7}(:,3)');
      elseif(k==8)
        set(h,'XDataSource','TimeLatLong{8}(:,2)'); set(h,'YDataSource','TimeLatLong{8}(:,3)');
      elseif(k==9)
        set(h,'XDataSource','TimeLatLong{9}(:,2)'); set(h,'YDataSource','TimeLatLong{9}(:,3)');
      elseif(k==10)
        set(h,'XDataSource','TimeLatLong{10}(:,2)'); set(h,'YDataSource','TimeLatLong{10}(:,3)');
      elseif(k==11)
        set(h,'XDataSource','TimeLatLong{11}(:,2)'); set(h,'YDataSource','TimeLatLong{11}(:,3)');
      elseif(k==12)
        set(h,'XDataSource','TimeLatLong{12}(:,2)'); set(h,'YDataSource','TimeLatLong{12}(:,3)');
      elseif(k==13)
        set(h,'XDataSource','TimeLatLong{13}(:,2)'); set(h,'YDataSource','TimeLatLong{13}(:,3)');
      elseif(k==14)
        set(h,'XDataSource','TimeLatLong{14}(:,2)'); set(h,'YDataSource','TimeLatLong{14}(:,3)');
      elseif(k==15)
        set(h,'XDataSource','TimeLatLong{15}(:,2)'); set(h,'YDataSource','TimeLatLong{15}(:,3)');
      elseif(k==16)
        set(h,'XDataSource','TimeLatLong{16}(:,2)'); set(h,'YDataSource','TimeLatLong{16}(:,3)');
      elseif(k==17)
        set(h,'XDataSource','TimeLatLong{17}(:,2)'); set(h,'YDataSource','TimeLatLong{17}(:,3)');
      end
      linkdata on;
    end

    % Plot ENU errors
    figure;
    subplot(3,1,1)
    plot(Merge(:,key.time),ENU(:,1)); % East
%    plot(Merge(:,key.time),VEast); % East

    ylabel('East Error [m]');
    grid
    title([DataStr,' - ',Section]);
    ENMax = HorzAxisMax;
    UMax = VertAxisMax;
    A = axis;
    A(3) = -ENMax;
    A(4) =  ENMax;
    axis(A);
    disp(['Mean East Error: ',num2str(mean(ENU(GoodTimeTag,1)),'%.3f')]);

    subplot(3,1,2)
    plot(Merge(:,key.time),ENU(:,2)); % North
%    plot(Merge(:,key.time),VNorth); % North

    ylabel('North Error [m]');
    grid
    A = axis;
    A(3) = -ENMax;
    A(4) =  ENMax;
    axis(A);
    subplot(3,1,3)
    disp(['Mean North Error: ',num2str(mean(ENU(GoodTimeTag,2)),'%.3f')]);

    plot(Merge(:,key.time),ENU(:,3)); % Up
%    plot(Merge(:,key.time),VUp); % Up

    ylabel('Height Error [m]');
    grid
    A = axis;
    A(3) = -2.0 * UMax;
    A(4) =  2.0 * UMax;
    axis(A);
    xlabel('Time [GPS Seconds]');

    outputFileName = [Section '_' DataStr '_ENU.png'];
    print ('-dpng','-r300',strrep(outputFileName,'\',''));
    if(saveFig)
      % savePlot.m not yet updated for subplots so save the file with
      % code in this script
      filenameFig = regexprep(strrep(outputFileName,'\',''),'.png','.fig');
      saveas(gcf,strrep(filenameFig,'\',''),'fig');
    end    
    
    %% Eric - AINS Generate ENU Errors Plot
    figure;
    subplot(3,1,1)
    plot(Merge(:,key.time),ENU(:,1)); % East
    hold all;
    plot(Merge(:,key.time), Merge(:,key.DUT_sigEast),'r');
    plot(Merge(:,key.time),-Merge(:,key.DUT_sigEast),'r');
%     legend('Error','Estimated Error');

    ylabel('East Error [m]');
    grid
    title(['Position Error - ',DataStr,' - ',Section]);
    ENMax = HorzAxisMax;
    UMax = VertAxisMax;
    A = axis;
    A(3) = -ENMax;
    A(4) =  ENMax;
    axis(A);
%     disp(['Mean East Error: ',num2str(mean(ENU(GoodTimeTag,1)),'%.3f')]);

    subplot(3,1,2)
    plot(Merge(:,key.time),ENU(:,2)); % North
    hold all;
    plot(Merge(:,key.time), Merge(:,key.DUT_sigNorth),'r');
    plot(Merge(:,key.time),-Merge(:,key.DUT_sigNorth),'r');

    ylabel('North Error [m]');
    grid
    A = axis;
    A(3) = -ENMax;
    A(4) =  ENMax;
    axis(A);
    
    subplot(3,1,3)
%     disp(['Mean North Error: ',num2str(mean(ENU(GoodTimeTag,2)),'%.3f')]);

    plot(Merge(:,key.time),ENU(:,3)); % Up
    hold all;
    plot(Merge(:,key.time), Merge(:,key.DUT_sigUp),'r');
    plot(Merge(:,key.time),-Merge(:,key.DUT_sigUp),'r');

    ylabel('Height Error [m]');
    grid
    A = axis;
    A(3) = -2.0 * UMax;
    A(4) =  2.0 * UMax;
    axis(A);
    xlabel('Time [GPS Seconds]');

    outputFileName = ['AINS_' Section '_' DataStr '_ENU.png'];
    print ('-dpng','-r300',strrep(outputFileName,'\',''));
    if(saveFig)
      % savePlot.m not yet updated for subplots so save the file with
      % code in this script
      filenameFig = regexprep(strrep(outputFileName,'\',''),'.png','.fig');
      saveas(gcf,strrep(filenameFig,'\',''),'fig');
    end
    
    %% Eric - AINS Generate Vel Errors Plot
    figure;
    subplot(3,1,1)
    plot(Merge(:,key.time), Merge(:,key.truth_eVel) - Merge(:,key.DUT_eVel),'g'); % East
    hold all;
    plot(Merge(:,key.time), Merge(:,key.DUT_eVelSig),'r');
    plot(Merge(:,key.time),-Merge(:,key.DUT_eVelSig),'r');
%     legend('Error','Estimated Error');

    ylabel('East Error [m/s]');
    grid
    title(['Velocity Error - ',DataStr,' - ',Section]);
    Ax = axis;
    if exist('VelHorxAxisMax','var')
      Ax(3) = -VelHorzAxisMax;
      Ax(4) =  VelHorzAxisMax;
    else
      Ax(3) = -VelAxisMax;
      Ax(4) =  VelAxisMax;
    end
    axis(Ax)
%     disp(['Mean East Error: ',num2str(mean(ENU(GoodTimeTag,1)),'%.3f')]);

    subplot(3,1,2)
    plot(Merge(:,key.time), Merge(:,key.truth_nVel) - Merge(:,key.DUT_nVel),'g'); % North
    hold all;
    plot(Merge(:,key.time), Merge(:,key.DUT_nVelSig),'r');
    plot(Merge(:,key.time),-Merge(:,key.DUT_nVelSig),'r');

    ylabel('North Error [m/s]');
    grid
    Ax = axis;
    if exist('VelHorzAxisMax','var')
      Ax(3) = -VelHorzAxisMax;
      Ax(4) =  VelHorzAxisMax;
    else
      Ax(3) = -VelAxisMax;
      Ax(4) =  VelAxisMax;
    end
    axis(Ax)
    
    subplot(3,1,3)
%     disp(['Mean North Error: ',num2str(mean(ENU(GoodTimeTag,2)),'%.3f')]);

    plot(Merge(:,key.time), Merge(:,key.truth_uVel) - Merge(:,key.DUT_uVel),'g'); % Up
    hold all;
    plot(Merge(:,key.time), Merge(:,key.DUT_uVelSig),'r');
    plot(Merge(:,key.time),-Merge(:,key.DUT_uVelSig),'r');

    ylabel('Height Error [m/s]');
    grid
    Ax = axis;
    if exist('VelVertAxisMax','var')
      Ax(3) = -VelVertAxisMax;
      Ax(4) =  VelVertAxisMax;
    else
      Ax(3) = -VelAxisMax;
      Ax(4) =  VelAxisMax;
    end
    axis(Ax)
    xlabel('Time [GPS Seconds]');

    outputFileName = ['AINS_' Section '_' DataStr '_Vel.png'];
    print ('-dpng','-r300',strrep(outputFileName,'\',''));
    if(saveFig)
      % savePlot.m not yet updated for subplots so save the file with
      % code in this script
      filenameFig = regexprep(strrep(outputFileName,'\',''),'.png','.fig');
      saveas(gcf,strrep(filenameFig,'\',''),'fig');
    end
       
    
    %% AINS Generate Attitude Errors Plot
    if(  (isnan(InputDataName{k}{20}) == 0) ...
      && (isnan(InputDataName{k}{21}) == 0) ...
      && (isnan(InputDataName{k}{22}) == 0) ...
      && (isnan(dRoll) == 0) ...
      && (ProcessAttitude == 1) )
    
      figure;
      subplot(3,1,1)
      plot(Merge(:,key.time), RollErr, 'k'); 
      hold all;
      plot(Merge(:,key.time), Merge(:,key.DUT_sigRoll),'r');
      plot(Merge(:,key.time),-Merge(:,key.DUT_sigRoll),'r');

      ylabel('Roll Error [Deg]');
      grid
      title(['Attitude Error - ',DataStr,' - ',Section]);
      Ax = axis;
      Ax(3) = -RPAxisMax;
      Ax(4) =  RPAxisMax;
      axis(Ax);

      subplot(3,1,2)
      plot(Merge(:,key.time), PitchErr, 'k'); 
      hold all;
      plot(Merge(:,key.time), Merge(:,key.DUT_sigPitch),'r');
      plot(Merge(:,key.time),-Merge(:,key.DUT_sigPitch),'r');

      ylabel('Pitch Error [Deg]');
      grid
      Ax = axis;
      Ax(3) = -RPAxisMax;
      Ax(4) =  RPAxisMax;
      axis(Ax);
    
      subplot(3,1,3)

      plot(Merge(:,key.time), HeadingErr, 'k'); 
      hold all;
      plot(Merge(:,key.time), Merge(:,key.DUT_sigHeading),'r');
      plot(Merge(:,key.time),-Merge(:,key.DUT_sigHeading),'r');

      ylabel('Heading Error [Deg]');
      grid
      Ax = axis;
      Ax(3) = -HeadAxisMax;
      Ax(4) =  HeadAxisMax;
      axis(Ax);
      xlabel('Time [GPS Seconds]');

      outputFileName = ['AINS_' Section '_' DataStr '_Att.png'];
      print ('-dpng','-r300',strrep(outputFileName,'\',''));
      if(saveFig)
        % savePlot.m not yet updated for subplots so save the file with
        % code in this script
        filenameFig = regexprep(strrep(outputFileName,'\',''),'.png','.fig');
        saveas(gcf,strrep(filenameFig,'\',''),'fig');
      end
    
    end
    
    %% Eric - Mean, std and rms
    % mean_std_rms = zeros(length(InputDataName),9,3);  
    
    % Pos - East, North and Up
    mean_std_rms(k,1,1) = mean(ENU(GoodTimeTag,1));
    mean_std_rms(k,2,1) = mean(ENU(GoodTimeTag,2));
    mean_std_rms(k,3,1) = mean(ENU(GoodTimeTag,3));
    mean_std_rms(k,1,2) = std(ENU(GoodTimeTag,1));
    mean_std_rms(k,2,2) = std(ENU(GoodTimeTag,2));
    mean_std_rms(k,3,2) = std(ENU(GoodTimeTag,3));
    mean_std_rms(k,1,3) = rms(ENU(GoodTimeTag,1));
    mean_std_rms(k,2,3) = rms(ENU(GoodTimeTag,2));
    mean_std_rms(k,3,3) = rms(ENU(GoodTimeTag,3));
    
    % Vel - East, North and Up
    eVelErr = Merge(GoodTimeTag,key.truth_eVel) - Merge(GoodTimeTag,key.DUT_eVel);
    nVelErr = Merge(GoodTimeTag,key.truth_nVel) - Merge(GoodTimeTag,key.DUT_nVel);
    uVelErr = Merge(GoodTimeTag,key.truth_uVel) - Merge(GoodTimeTag,key.DUT_uVel);
    mean_std_rms(k,4,1) = mean(eVelErr);
    mean_std_rms(k,5,1) = mean(nVelErr);
    mean_std_rms(k,6,1) = mean(uVelErr);
    mean_std_rms(k,4,2) = std(eVelErr);
    mean_std_rms(k,5,2) = std(nVelErr);
    mean_std_rms(k,6,2) = std(uVelErr);
    mean_std_rms(k,4,3) = rms(eVelErr);
    mean_std_rms(k,5,3) = rms(nVelErr);
    mean_std_rms(k,6,3) = rms(uVelErr);
    
    % Att - Roll, Pitch and Yaw
    if(  (isnan(InputDataName{k}{20}) == 0) ...
      && (isnan(InputDataName{k}{21}) == 0) ...
      && (isnan(InputDataName{k}{22}) == 0) ...
      && (isnan(dRoll) == 0) ...
      && (ProcessAttitude == 1) )
      mean_std_rms(k,7,1) = mean(RollErr);
      mean_std_rms(k,8,1) = mean(PitchErr);
      mean_std_rms(k,9,1) = mean(HeadingErr);
      mean_std_rms(k,7,2) = std(RollErr);
      mean_std_rms(k,8,2) = std(PitchErr);
      mean_std_rms(k,9,2) = std(HeadingErr);
      mean_std_rms(k,7,3) = rms(RollErr);
      mean_std_rms(k,8,3) = rms(PitchErr);
      mean_std_rms(k,9,3) = rms(HeadingErr);
    end
    %% Eric - Plots - Mean, std and rms - End

    if(plotMap)
      % Now plot a map with the data coloured based on the Height Track
      % Error
      figure
      current = colormap();
      colormap('jet');
      plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                     TimeLatLong{k}(1:DownSample:end,2), ...
                     ENU(1:DownSample:end,3), ...
                     [DataStr,' - ',Section,' Height Error [m]'], ...
                     VertAxisMax);

      outputFileName = [Section '_' DataStr '_HeightErrMap.png'];
      savePlot(strrep(outputFileName,'\',''),saveFig);
      colormap(current);
    end

    disp(['Mean Up Error: ',num2str(mean(ENU(GoodTimeTag,3)),'%.3f')]);

    if( (zoomLink ==1) && (length(InputDataName) <= 17) )
      % set up the call back function in MapZoomCallback.m, this will
      % get called whether the user releases the zoom on any of the
      % subplots of the E/N/U graph
      f = zoom(gcf);
      set(f,'ActionPostCallback',@(obj,event_obj)MapZoomCallback(obj,event_obj,MapFig,k ) );
    end

    figure
    plot(ENU(:,1),ENU(:,2))
    A = axis;
    A(1) = -ENMax;
    A(2) =  ENMax;
    A(3) = -ENMax;
    A(4) =  ENMax;
    axis(A)
    axis square;
    grid;
    xlabel('East Error [m]');
    ylabel('North Error [m]');
    title([DataStr,' - ',Section]);
    outputFileName = [Section '_' DataStr '_EN.png'];
    savePlot(strrep(outputFileName,'\',''),saveFig);

    if(ENMax < 1)
      Step = 0.025;
    else
      Step = 0.1;
    end

    if(exist('hist3') == 2)
      errEN(ENU(:,1),ENU(:,2),[Section '_' DataStr '_EN'], [DataStr,'-',Section] ,ENMax,Step,0);
      errIntegrityPlot( sqrt(ENU(:,1).^2 + ENU(:,2).^2), ...
                        EN_Sigma, ...
                        [Section '_' DataStr], ...
                        [DataStr,'-',Section], ...
                        ENMax * 4, ...
                        Step);
    end

    if(plotMap)
      % Now plot a map with the data coloured based on the Cross Track
      % Error
      figure
      current = colormap();
      colormap('jet');
      plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                     TimeLatLong{k}(1:DownSample:end,2), ...
                     sqrt(ENU(1:DownSample:end,1).^2 + ENU(1:DownSample:end,2).^2), ...
                     [DataStr,' - ',Section,' 2D Error [m]'], ...
                     HorzAxisMax);

      outputFileName = [Section '_' DataStr '_ENErrMap.png'];
      savePlot(strrep(outputFileName,'\',''),saveFig);
      colormap(current);
    end

    if(  (isnan(InputDataName{k}{14}) == 0) ...
      && (isnan(dEastVel) == 0) )
        
      VEast = Merge(:,key.truth_eVel) - Merge(:,key.DUT_eVel);
      VNorth = Merge(:,key.truth_nVel) - Merge(:,key.DUT_nVel);
      VUp = Merge(:,key.truth_uVel) - Merge(:,key.DUT_uVel);
      HzVel = sqrt(VEast.^2 + VNorth.^2)';
      UpV  = abs(VUp)';
        
      if(plotMap)
        % Now plot a map with the data coloured based on the Horizontal Velocity Error
        figure
        current = colormap();
        colormap('jet');
        plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                       TimeLatLong{k}(1:DownSample:end,2), ...
                       HzVel(1:DownSample:end), ...
                       [DataStr,' - ',Section,' Horizontal Velocity Error [m/s]'], ...
                       HorzAxisMax);
        outputFileName = [Section '_' DataStr '_HorzVelErrMap.png'];
        savePlot(strrep(outputFileName,'\',''),saveFig);
        colormap(current);

        % Now plot a map with the data coloured based on the Vertical Velocity Error
        figure
        current = colormap();
        colormap('jet');
        plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                       TimeLatLong{k}(1:DownSample:end,2), ...
                       UpV(1:DownSample:end), ...
                       [DataStr,' - ',Section,' Vertical Velocity Error [m/s]'], ...
                       HorzAxisMax);
        outputFileName = [Section '_' DataStr '_VertVelErrMap.png'];
        savePlot(strrep(outputFileName,'\',''),saveFig);
        colormap(current);
      end % if(plotMap)
    end

    % If heading is available...
    if(isnan(dHeading) == 0)
      
      % Compute the vehicle frame errors
      Xv = V2N11.*ENU(:,1) + V2N21.*ENU(:,2) + V2N31.*ENU(:,3);
      Yv = V2N12.*ENU(:,1) + V2N22.*ENU(:,2) + V2N32.*ENU(:,3);
      Zv = V2N13.*ENU(:,1) + V2N23.*ENU(:,2) + V2N33.*ENU(:,3);
      
      Xvrms = rms(Xv(GoodTimeTag));
      Yvrms = rms(Yv(GoodTimeTag));
      Zvrms = rms(Zv(GoodTimeTag));
      
      % Plot cross track error
      figure;
      subplot(3,1,1);
      plot(Merge(:,key.time),Yv); 
      ylabel('Cross track error [m]');
      grid;
      title([DataStr,' - ',Section]);
      A = axis;
      A(3) = -5*Yvrms;
      A(4) = -A(3);
      axis(A);
      disp(['Mean cross track error: ',num2str(mean(Yv(GoodTimeTag)),'%.3f')]);

      % Plot along track error
      subplot(3,1,2);
      plot(Merge(:,key.time),Xv); 
      ylabel('Along track error [m]');
      grid
      A = axis;
      A(3) = -5*Xvrms;
      A(4) = -A(3);
      axis(A);
      disp(['Mean along track error: ',num2str(mean(Xv(GoodTimeTag)),'%.3f')]);

      % Plot height error
      subplot(3,1,3);
      plot(Merge(:,key.time),Zv); 
      ylabel('Height error [m]');
      grid
      A = axis;
      A(3) = -5*Zvrms;
      A(4) = -A(3);
      axis(A);
      disp(['Mean height error: ',num2str(mean(Zv(GoodTimeTag)),'%.3f')]);
      xlabel('Time [GPS Seconds]');
      outputFileName = [Section '_' DataStr '_CAU.png'];
      % savePlot(strrep(outputFileName,'\',''),saveFig);

      print ('-dpng','-r300',strrep(outputFileName,'\',''));

      if(saveFig)
        % savePlot.m not yet updated for subplots so save the file with
        % code in this script
        filenameFig = regexprep(strrep(outputFileName,'\',''),'.png','.fig');
        saveas(gcf,strrep(filenameFig,'\',''),'fig');
      end

      if(plotMap)
        % Now plot a map with the data coloured based on the Cross Track
        % Error
        figure
        current = colormap();
        colormap('jet');
        plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                       TimeLatLong{k}(1:DownSample:end,2), ...
                       Yv(1:DownSample:end), ...
                       [DataStr,' - ',Section,' Cross Track [m]'], ...
                       HorzAxisMax);

        outputFileName = [Section '_' DataStr '_CrossTrackMap.png'];
        savePlot(strrep(outputFileName,'\',''),saveFig);
        colormap(current);

        % Now plot a map with the data coloured based on the Along Track
        % Error
        figure
        current = colormap();
        colormap('jet');
        plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                       TimeLatLong{k}(1:DownSample:end,2), ...
                       Xv(1:DownSample:end), ...
                       [DataStr,' - ',Section,' Along Track [m]'], ...
                       HorzAxisMax);

        outputFileName = [Section '_' DataStr '_AlongTrackMap.png'];
        savePlot(strrep(outputFileName,'\',''),saveFig);
        colormap(current);
      end

      % Plot cross track / along track scatter plot
      figure;
      plot(Yv,Xv)
      %CAMax = 5*max([Crms Arms]);
      CAMax = HorzAxisMax;
      A = axis;
      A(1) = -CAMax;
      A(2) =  CAMax;
      A(3) = -CAMax;
      A(4) =  CAMax;
      axis(A)
      axis square;
      grid;
      xlabel('Cross track error [m]');
      ylabel('Along track error [m]');
      title([DataStr,' - ',Section]);
      outputFileName = [Section '_' DataStr '_CA.png'];
      savePlot(strrep(outputFileName,'\',''),saveFig);

%
% Not sure why the vehicle frame errors are recomputed in this way. 
% Why not just use the absolute value of the previously computed
% errors?
%
% Original code:
%       CrossTrack = sqrt((SinHead(GoodTimeTag) .* ENU(GoodTimeTag,2)) .^ 2 ...
%                       + (CosHead(GoodTimeTag) .* ENU(GoodTimeTag,1)) .^2);
%       AlongTrack = sqrt((CosHead(GoodTimeTag) .* ENU(GoodTimeTag,2)) .^ 2 ...
%                       + (SinHead(GoodTimeTag) .* ENU(GoodTimeTag,1)) .^2);
%
% New code with full rotation matrix:
%       Xv = sqrt ( (V2N11(GoodTimeTag).*ENU(GoodTimeTag,1)) .^2 ...
%                 + (V2N21(GoodTimeTag).*ENU(GoodTimeTag,2)) .^2 ...
%                 + (V2N31(GoodTimeTag).*ENU(GoodTimeTag,3)) .^2 );
%       Yv = sqrt ( (V2N12(GoodTimeTag).*ENU(GoodTimeTag,1)) .^2 ...
%                 + (V2N22(GoodTimeTag).*ENU(GoodTimeTag,2)) .^2 ...
%                 + (V2N32(GoodTimeTag).*ENU(GoodTimeTag,3)) .^2 );
%       Zv = sqrt ( (V2N13(GoodTimeTag).*ENU(GoodTimeTag,1)) .^2 ...
%                 + (V2N23(GoodTimeTag).*ENU(GoodTimeTag,2)) .^2 ...
%                 + (V2N33(GoodTimeTag).*ENU(GoodTimeTag,3)) .^2 );
%
      % Compute the absolute value of the previously computed
      % vehicle frame errors
      Xv = abs(Xv);
      Yv = abs(Yv);
      Zv = abs(Zv);
      
      % Display the cross-track (Yv) percentiles
      CData = sort(Yv(GoodTimeTag));
      disp([DataStr, ' Epochs :',num2str(length(ENU(GoodTimeTag)))]);
      disp('Cross-Track Stats:');
      disp([DataStr, '( 68%):',num2str(CData( floor(0.68*length(CData))))]);
      disp([DataStr, '( 95%):',num2str(CData( floor(0.95*length(CData))))]);
      disp([DataStr, '(100%):',num2str(CData(end))]);
      disp(' ');
      if(length(InputDataName) <= 5)
        CrossTrackSummary{k} = [10 ' 68%:',num2str(CData( floor(0.68*length(CData))),2),'m' 10 ' 95%:',num2str(CData( floor(0.95*length(CData))),2),'m' 10 '100%:',num2str(CData(end),2),'m'];
      else
        CrossTrackSummary{k} = ' ';
      end

      % Display the along-track (Xv) percentiles
      AData = sort(Xv(GoodTimeTag));
      disp([DataStr, ' Epochs :',num2str(length(ENU(GoodTimeTag)))]);
      disp('Along-Track Stats:');
      disp([DataStr, '( 68%):',num2str(AData( floor(0.68*length(AData))))]);
      disp([DataStr, '( 95%):',num2str(AData( floor(0.95*length(AData))))]);
      disp([DataStr, '(100%):',num2str(AData(end))]);
      disp(' ');
      if(length(InputDataName) <= 5)
        AlongTrackSummary{k} = [10 ' 68%:',num2str(AData( floor(0.68*length(AData))),2),'m' 10 ' 95%:',num2str(AData( floor(0.95*length(AData))),2),'m' 10 '100%:',num2str(AData(end),2),'m'];
      else
        AlongTrackSummary{k} = ' ';
      end

      CrossStats(k,1) = CData(floor(0.50*length(CData)));
      CrossStats(k,2) = CData(floor(0.68*length(CData)));
      CrossStats(k,3) = CData(floor(0.95*length(CData)));
      CrossStats(k,4) = CData(floor(0.99*length(CData)));
      CrossStats(k,5) = CData(end);
      
      AlongStats(k,1) = AData(floor(0.50*length(AData)));
      AlongStats(k,2) = AData(floor(0.68*length(AData)));
      AlongStats(k,3) = AData(floor(0.95*length(AData)));
      AlongStats(k,4) = AData(floor(0.99*length(AData)));
      AlongStats(k,5) = AData(end);
      
      % Now plot a map with the data coloured based on the Roll Error
      if(max(isnan(RollStats(k,:)))==0 & plotMap & ProcessAttitude == 1)
          figure
          current = colormap();
          colormap('jet');
          plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                         TimeLatLong{k}(1:DownSample:end,2), ...
                         abs(RollErr(1:DownSample:end)), ...
                         [DataStr,' - ',Section,' Roll Error [deg]'], ...
                         HorzAxisMax);
          outputFileName = [Section '_' DataStr '_RollErrMap.png'];
          savePlot(strrep(outputFileName,'\',''),saveFig);
          colormap(current);
      end;
      
      % Now plot a map with the data coloured based on the Pitch Error
      if(max(isnan(PitchStats(k,:)))==0 & plotMap & ProcessAttitude == 1)
          figure
          current = colormap();
          colormap('jet');
          plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                         TimeLatLong{k}(1:DownSample:end,2), ...
                         abs(PitchErr(1:DownSample:end)), ...
                         [DataStr,' - ',Section,' Pitch Error [deg]'], ...
                         HorzAxisMax);
          outputFileName = [Section '_' DataStr '_PitchErrMap.png'];
          savePlot(strrep(outputFileName,'\',''),saveFig);
          colormap(current);
      end;
      
      % Now plot a map with the data coloured based on the Heading Error
      if(max(isnan(HeadingStats(k,:)))==0 & plotMap & ProcessAttitude == 1)
          figure
          current = colormap();
          colormap('jet');
          plotColourMap( TimeLatLong{k}(1:DownSample:end,3), ...
                         TimeLatLong{k}(1:DownSample:end,2), ...
                         abs(HeadingErr(1:DownSample:end)), ...
                         [DataStr,' - ',Section,' Heading Error [deg]'], ...
                         HorzAxisMax);
          outputFileName = [Section '_' DataStr '_HeadErrMap.png'];
          savePlot(strrep(outputFileName,'\',''),saveFig);
          colormap(current);
      end
    else
      CrossStats(k,1:5) = NaN;
      AlongStats(k,1:5) = NaN;
    end
    
    % Display the percentiles
    VData = sort(abs(ENU(GoodTimeTag,3)));
    disp([DataStr, ' Epochs :',num2str(length(ENU(GoodTimeTag)))]);
    disp('Vertical Stats:');
    disp([DataStr, '( 68%):',num2str(VData( floor(0.68*length(VData))))]);
    disp([DataStr, '( 95%):',num2str(VData( floor(0.95*length(VData))))]);
    disp([DataStr, '(100%):',num2str(VData(end))]);
    disp(' ');
    if(length(InputDataName) <= 5)
      VertSummary{k} = [10 ' 68%:',num2str(VData( floor(0.68*length(VData))),2),'m' 10 ' 95%:',num2str(VData( floor(0.95*length(VData))),2),'m' 10 '100%:',num2str(VData(end),2),'m'];
    else
      VertSummary{k} = ' ';
    end

    VertStats(k,1) = VData(floor(0.50*length(VData)));
    VertStats(k,2) = VData(floor(0.68*length(VData)));
    VertStats(k,3) = VData(floor(0.95*length(VData)));
    VertStats(k,4) = VData(floor(0.99*length(VData)));
    VertStats(k,5) = VData(end);
    
    Err2D = sqrt(ENU(GoodTimeTag,1).^2 + ENU(GoodTimeTag,2).^2);
    % Store 2D error distribution, pre sort to make the CDF plotting easier
    TwoD = sort(Err2D);
    
    TwoDStats(k,1) = TwoD(floor(0.50*length(TwoD)));
    TwoDStats(k,2) = TwoD(floor(0.68*length(TwoD)));
    TwoDStats(k,3) = TwoD(floor(0.95*length(TwoD)));
    TwoDStats(k,4) = TwoD(floor(0.99*length(TwoD)));
    TwoDStats(k,5) = TwoD(end);
    
    % Eric - Add Position Stats East and North
    EastP = sort(abs(ENU(GoodTimeTag,1)));    
    EastStats(k,1) = EastP(floor(0.50*length(EastP)));
    EastStats(k,2) = EastP(floor(0.68*length(EastP)));
    EastStats(k,3) = EastP(floor(0.95*length(EastP)));
    EastStats(k,4) = EastP(floor(0.99*length(EastP)));
    EastStats(k,5) = EastP(end);
    
    NorthP = sort(abs(ENU(GoodTimeTag,2)));    
    NorthStats(k,1) = NorthP(floor(0.50*length(NorthP)));
    NorthStats(k,2) = NorthP(floor(0.68*length(NorthP)));
    NorthStats(k,3) = NorthP(floor(0.95*length(NorthP)));
    NorthStats(k,4) = NorthP(floor(0.99*length(NorthP)));
    NorthStats(k,5) = NorthP(end);
    
    % Now convert the 2D outliers into a Google Earth file
    if(ProcessOutliers == 1)
      KMLOut = 1;
      % Error is more than the defined threshold or more than 99%
      if((d2DOutlierThres == -1) && (d2DOutlierPer < 100) && (d2DOutlierPer > 0) )
        per = TwoD(floor(length(TwoD)*d2DOutlierPer/100));
        i = find( Err2D > per);
        Label = [Section,'_',DataStr,' _2D_Outliers > ',num2str(d2DOutlierPer),'% error'];
      elseif( (d2DOutlierPer == -1) && (d2DOutlierThres > 0) )
        i = find(Err2D > d2DOutlierThres);
        Label = [Section,'_',DataStr,' _2D_Outliers > ', num2str(d2DOutlierThres,'%.0f'),'m error'];
      elseif((d2DOutlierThres > 0) && (d2DOutlierPer < 100) && (d2DOutlierPer > 0) )
        per = TwoD(floor(length(TwoD)*d2DOutlierPer/100));
        i = find( (Err2D > d2DOutlierThres) | (Err2D > per));
        Label = [Section,'_',DataStr,' _2D_Outliers > ', num2str(d2DOutlierThres,'%.0f'),'m | > ',num2str(d2DOutlierPer),'% error'];
      else
        disp('Bad 2D outlier threshold');
        KMLOut = 0;
      end

      if( (KMLOut == 1) && (~isempty(i)) )
        writeKML([Section,'_',DataStr, '_2D_Outliers'], ...
                  DataStr, ...
                  Label, ...
                  Err2D(i), ...
                  Merge(i,key.time), ...
                  Merge(i,key.truth_lat), ...
                  Merge(i,key.truth_lon), ...
                  Merge(i,key.truth_hgt), ...
                  Merge(i,key.DUT_lat), ...
                  Merge(i,key.DUT_lon), ...
                  Merge(i,key.DUT_hgt));
      end
    end


    % Format the stats for the legend
    TwoDSummary{k} = [10 ' 68%:',num2str(TwoDStats(k,2),2),'m' 10 ' 95%:',num2str(TwoDStats(k,3),2),'m' 10 '100%:',num2str(TwoDStats(k,5),2),'m'];

    Data2D{k} = TwoD';

    if(isnan(dHeading) == 0)
      % Store the Cross error distribution ...
      CrossTrackData{k} = CData';
      % Store the Along error distribution ...
      AlongTrackData{k} = AData';
    end

    DUTData{k} = [ Merge(:,key.DUT_lat),   ...
                   Merge(:,key.DUT_lon),   ...
                   Merge(:,key.DUT_hgt),   ...
                   Merge(:,key.truth_lat), ...
                   Merge(:,key.truth_lon), ...
                   Merge(:,key.truth_hgt), ...
                   ENU(:,1), ENU(:,2),ENU(:,3)];

    % Store the height error distribution;
    HgtMod = abs(ENU(GoodTimeTag,3));
    HgtSort = sort(HgtMod)';
    if(length(InputDataName) <= 5)
      HeightSummary{k} = [10 ' 68%:',num2str(HgtSort( floor(0.68*length(HgtSort))),2),'m' 10 ' 95%:',num2str(HgtSort( floor(0.95*length(HgtSort))),2),'m' 10 '100%:',num2str(HgtSort(end),2),'m'];
    else
      HeightSummary{k} = ' ';
    end
    HeightData{k} = HgtSort;

    % Now convert the Hgt outliers into a Google Earth file
    if(ProcessOutliers == 1)
      KMLOut = 1;
      % Error is more than the defined threshold or more than 99%
      if((dHgtOutlierThres == -1) && (dHgtOutlierPer < 100) && (dHgtOutlierPer > 0) )
        per = HgtSort(floor(length(HgtSort)*dHgtOutlierPer/100));
        i = find( HgtMod > per);
        Label = [Section,'_',DataStr,' _Hgt_Outliers > ',num2str(dHgtOutlierPer),'% error'];
      elseif( (dHgtOutlierPer == -1) && (dHgtOutlierThres > 0) )
        i = find(HgtMod > dHgtOutlierThres);
        Label = [Section,'_',DataStr,' _Hgt_Outliers > ', num2str(dHgtOutlierThres,'%.0f'),'m error'];
      elseif((dHgtOutlierThres > 0) && (dHgtOutlierPer < 100) && (dHgtOutlierPer > 0) )
        per = HgtSort(floor(length(HgtSort)*dHgtOutlierPer/100));
        i = find( (HgtMod > dHgtOutlierThres) | (HgtMod > per));
        Label = [Section,'_',DataStr,' _Hgt_Outliers > ', num2str(dHgtOutlierThres,'%.0f'),'m | > ',num2str(dHgtOutlierPer),'% error'];
      else
        disp('Bad Hgt outlier threshold');
        KMLOut = 0;
      end

      if( (KMLOut == 1) && (~isempty(i)) )
        writeKML([Section,'_',DataStr, '_Hgt_Outliers'], ...
                  DataStr, ...
                  Label, ...
                  HgtMod(i), ...
                  Merge(i,key.time), ...
                  Merge(i,key.truth_lat), ...
                  Merge(i,key.truth_lon), ...
                  Merge(i,key.truth_hgt), ...
                  Merge(i,key.DUT_lat), ...
                  Merge(i,key.DUT_lon), ...
                  Merge(i,key.DUT_hgt));
      end
    end

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

    % If requested close all the figures
    if(closeOnComplete == 1)
      close all;
    end
end %for k=1:length(InputDataName)

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

if(isnan(dHeading) == 0)
  % Cross-Track Cumulative Distribution 
  createCDF(CrossTrackData, ...
            NameArray, ...
            NaN, ... %CrossTrackSummary, ...
            colorArray, ...
            HorzAxisMax, ...
            0.0, ...
            'Cross Track Error [m]', ...
            ['Cumulative Distribution : Cross Track CDF (', Section,')'], ...
            [Section, '_CrossTrackCDF.'], ...
            saveFig, ...
            addPercentiles);

  createCDF(CrossTrackData, ...
            NameArray, ...
            NaN, ... %CrossTrackSummary, ...
            colorArray, ...
            HorzAxisMax*5, ...
            0.95, ...
            'Cross Track Error [m]', ...
            ['Cumulative Distribution : Cross Track CDF (', Section,')'], ...
            [Section, '_CrossTrackCDF_95pcnt.'], ...
            saveFig, ...
            addPercentiles);

  % Along-Track Cumulative Distribution 
  createCDF(AlongTrackData, ...
            NameArray, ...
            NaN, ... %AlongTrackSummary, ...
            colorArray, ...
            HorzAxisMax, ...
            0.0, ...
            'Along Track Error [m]', ...
            ['Cumulative Distribution : Along Track CDF (', Section,')'], ...
            [Section, '_AlongTrackCDF.'], ...
            saveFig, ...
            addPercentiles);

  createCDF(AlongTrackData, ...
            NameArray, ...
            NaN, ... %AlongTrackSummary, ...
            colorArray, ...
            HorzAxisMax*5, ...
            0.95, ...
            'Along Track Error [m]', ...
            ['Cumulative Distribution : Along Track CDF (', Section,')'], ...
            [Section, '_AlongTrackCDF_95pcnt.'], ...
            saveFig, ...
            addPercentiles);



  % Now generate the penalized CDF. For data sets with yield less than
  % 100% we pad the end of the array with NaNs so the CDF can never get
  % to 100%.

  for i=1:k
    if( (abs(Yield(i) - 100) * 1e6) > 1)
      % Need to pad as the yield is not 100%
      Pad = floor(length(CrossTrackData{i}) / (Yield(i)/100) - length(CrossTrackData{i}));
      if(Pad > 0)
        CrossTrackData{i} = [CrossTrackData{i} ones(1,Pad)*NaN];
        AlongTrackData{i} = [AlongTrackData{i} ones(1,Pad)*NaN];
      end
    end
  end

 
  createCDF(CrossTrackData, ...
            NameArray, ...
            NaN, ... %CrossTrackSummary, ...
            colorArray, ...
            HorzAxisMax, ...
            0.0, ...
            'Cross Track Error [m]', ...
            ['Penalized CDF : Cross Track (', Section,')'], ...
            [Section, '_PenalizedCrossTrackCDF.'], ...
            saveFig, ...
            addPercentiles);

  createCDF(AlongTrackData, ...
            NameArray, ...
            NaN, ... %AlongTrackSummary, ...
            colorArray, ...
            HorzAxisMax, ...
            0.0, ...
            'Along Track Error [m]', ...
            ['Penalized CDF : Along Track (', Section,')'], ...
            [Section, '_PenalizedAlongTrackCDF.'], ...
            saveFig, ...
            addPercentiles);


end

% 2D cumulative Distribution
createCDF(Data2D, ...
          NameArray, ...
          NaN, ... %TwoDSummary, ...
          colorArray, ...
          HorzAxisMax, ...
          0.0, ...
          '2D Error [m]', ...
          ['Cumulative Distribution : 2D CDF (', Section,')'], ...
          [Section, '_2DCDF.'], ...
          saveFig, ...
          addPercentiles);

% 2D cumulative Distribution - zoom in on 95%
createCDF(Data2D, ...
          NameArray, ...
          NaN, ... %TwoDSummary, ...
          colorArray, ...
          HorzAxisMax*5, ...
          0.95, ...
          '2D Error [m]', ...
          ['Cumulative Distribution : 2D CDF (', Section,')'], ...
          [Section, '_2DCDF_95pcnt.'], ...
          saveFig, ...
          addPercentiles);


% Height CDF
createCDF(HeightData, ...
          NameArray, ...
          NaN, ... %HeightSummary, ...
          colorArray, ...
          VertAxisMax, ...
          0.0, ...
          'abs(Height Error) [m]', ...
          ['Cumulative Distribution : abs(Height) CDF (', Section,')'], ...
          [Section, '_HgtCDF.'], ...
          saveFig, ...
          addPercentiles);


% Now generate the penalized CDF. For data sets with yield less than
% 100% we pad the end of the array with NaNs so the CDF can never get
% to 100%.

for i=1:k
  if( (abs(Yield(i) - 100) * 1e6) > 1)
    % Need to pad as the yield is not 100%
    Pad = floor(length(Data2D{i}) / (Yield(i)/100) - length(Data2D{i}));
    if(Pad > 0)
      Data2D{i}     = [Data2D{i} ones(1,Pad)*NaN];
      HeightData{i} = [HeightData{i} ones(1,Pad)*NaN];
    end
  end
end

% 2D cumulative Distribution
createCDF(Data2D, ...
          NameArray, ...
          NaN, ... %TwoDSummary, ...
          colorArray, ...
          HorzAxisMax, ...
          0.0, ...
          '2D Error [m]', ...
          ['Penalized CDF : 2D (', Section,')'], ...
          [Section, '_Penalized2DCDF.'], ...
          saveFig, ...
          addPercentiles);

% Height CDF
createCDF(HeightData, ...
          NameArray, ...
          NaN, ... %HeightSummary, ...
          colorArray, ...
          VertAxisMax, ...
          0.0, ...
          'abs(Height Error) [m]', ...
          ['Penalized CDF : abs(Height) CDF (', Section,')'], ...
          [Section, '_PenalizedHgtCDF.'], ...
          saveFig, ...
          addPercentiles);




if( (isnan(dEastVel) == 0) && exist('EastVel','var') )
  % Plot East/North/Up Velocity CDFs
  createCDF(EastVel, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            VelAxisMax, ...
            0.0, ...
            'abs(East Velocity Error) [m/s]', ...
            ['Cumulative Distribution : abs(East Vel) CDF (', Section,')'], ...
            [Section, '_EastVelCDF.'], ...
            saveFig, ...
            addPercentiles);
end            

if(   (isnan(dEastVel) == 0) ...
   && (isnan(dNorthVel) == 0) ...
   && exist('HorzVel','var') )
  % Plot Horizontal Velocity CDF
  createCDF(HorzVel, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            VelAxisMax, ...
            0.0, ...
            'Horizontal Velocity Error [m/s]', ...
            ['Cumulative Distribution : Horz Vel CDF (', Section,')'], ...
            [Section, '_HorzVelCDF.'], ...
            saveFig, ...
            addPercentiles);
end

if( (isnan(dNorthVel) == 0) && exist('NorthVel','var') )
  createCDF(NorthVel, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            VelAxisMax, ...
            0.0, ...
            'abs(North Velocity Error) [m/s]', ...
            ['Cumulative Distribution : abs(North Vel) CDF (', Section,')'], ...
            [Section, '_NorthVelCDF.'], ...
            saveFig, ...
            addPercentiles);
end


if( (isnan(dUpVel) == 0) && exist('UpVel','var') )
  createCDF(UpVel, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            VelAxisMax, ...
            0.0, ...
            'abs(Up Velocity Error) [m/s]', ...
            ['Cumulative Distribution : abs(Height Vel) CDF (', Section,')'], ...
            [Section, '_HgtVelCDF.'], ...
            saveFig, ...
            addPercentiles);
end


% Plot Roll/Pitch/Heading CDFs
if( (isnan(dRoll) == 0) && (exist('RollErrSort','var')) && (ProcessAttitude == 1) )
  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','var')) && (ProcessAttitude == 1) )
  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','var')) && (ProcessAttitude == 1) )
  createCDF(HeadingErrSort, ...
            NameArray, ...
            NaN, ...
            colorArray, ...
            HeadingCDFAxisMax, ...
            0.0, ...
            'abs(Heading Error) [deg]', ...
            ['Cumulative Distribution : abs(Heading) CDF (', Section,')'], ...
            [Section, '_HeadingCDF.'], ...
            saveFig, ...
            addPercentiles);
end

diary OFF;

%
% Collate the key stats as HTML (XXX_stats.html):-
%
HTML_name = [Section,'_stats.html'];
fid = fopen(HTML_name,'w');
fprintf(fid,'<!DOCTYPE html>\n');
fprintf(fid,'<html>\n');
fprintf(fid,'<head>\n');
fprintf(fid,'<title>%s - 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);

AttitudeAvailable = length(find(isnan(RollStats) == 0));

fprintf(fid,'<h3>Quick Links</h3>');
fprintf(fid,'<ul>\n');
fprintf(fid,'<li><a href="#POSPac">Reference POSPac</a></li>\n');
fprintf(fid,'<li><a href="#PvaResults">PVA Results</a></li>\n');
fprintf(fid,'<li><a href="#ResultPlots">Result Plot Table</a></li>\n');
fprintf(fid,'<li><a href="#Yield">Yield Information</a></li>\n');
fprintf(fid,'<li><a href="#CrossTrack">Cross Track Statistics</a></li>\n');
fprintf(fid,'<li><a href="#AlongTrack">Along Track Statistics</a></li>\n');
fprintf(fid,'<li><a href="#2D">2D Statistics</a></li>\n');
fprintf(fid,'<li><a href="#Hgt">Height Statistics</a></li>\n');
fprintf(fid,'<li><a href="#HorzVel">Horizontal Velocity Statistics</a></li>\n');
fprintf(fid,'<li><a href="#HgtVel">Vertical Velocity Statistics</a></li>\n');
if(AttitudeAvailable & (ProcessAttitude == 1) )
  fprintf(fid,'<li><a href="#Roll">Roll Statistics</a></li>\n');
  fprintf(fid,'<li><a href="#Pitch">Pitch Statistics</a></li>\n');
  fprintf(fid,'<li><a href="#Heading">Heading Statistics</a></li>\n');
end
fprintf(fid,'</ul>\n');
fprintf(fid,'<hr/>\n');

%% Eric - Add POSPac Head Title
fprintf(fid,'<a name="POSPac"/><h2>Reference - POSPac</h2>\n');

if(plotMap)
    File = [Section,'_POSPac - ',Section,'_Position_Map'];
    fprintf(fid,['<a href="',File,'.png"><img src="',File,'.png" height="400"/></a>']);
    if(saveFig)
      fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>&nbsp;']);
    end
    File = [Section,'_POSPac - ',Section,'_Position_satellite'];
    fprintf(fid,['<a href="',File,'.png"><img src="',File,'.png" height="400"/></a>']);
    if(saveFig)
      fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>&nbsp;']);
    end
    fprintf(fid,'<br/>');
end
File = [Section,'_POSPac_Error'];
fprintf(fid,['<a href="',File,'.png"><img src="',File,'.png" height="400"/></a>']);
if(saveFig)
  fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>&nbsp;']);
end
File = [Section,'_POSPac - ',Section,'_HorzVel'];
fprintf(fid,['<a href="',File,'.png"><img src="',File,'.png" height="400"/></a>']);
if(saveFig)
  fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>&nbsp;']);
end
fprintf(fid,'<br/>');

if(CreateKMZ == 1)
  fprintf(fid,['<a href="',Section,'-Truth.kmz">POSPac KMZ</a><br/>']);
end

File = [Section,'_POSPac - ',Section,'_Hgt'];
fprintf(fid,['<a href="',File,'.png">Height</a>']);
if(saveFig)
  fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
end
fprintf(fid,['<br/>']);

File = [Section,'_POSPac - ',Section,'_RollPitch'];
fprintf(fid,['<a href="',File,'.png">Roll-Pitch</a>']);
if(saveFig)
  fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
end
fprintf(fid,['<br/>']);

File = [Section,'_POSPac - ',Section,'_Heading'];
fprintf(fid,['<a href="',File,'.png">Heading</a>']);
if(saveFig)
  fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
end
fprintf(fid,['<br/>']);

File = [Section,'_POSPac - ',Section,'_Position'];
fprintf(fid,['<a href="',File,'.png">Position</a>']);
if(saveFig)
  fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
end
fprintf(fid,['<br/>']);

File = [Section,'_POSPac - ',Section,'_HorzVel_mph'];
fprintf(fid,['<a href="',File,'.png">HorzVel_mph</a>']);
if(saveFig)
  fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
end
fprintf(fid,['<br/>']);

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

fprintf(fid,'<a name="PvaResults"/><h2>PVA Results</h2>\n');

% Loop around for the 4 or 7 types of statistics we want to tabulate
if(AttitudeAvailable & (ProcessAttitude == 1) )
  MaxPlots = 9;
else
  MaxPlots = 6;
end

fprintf(fid,'<a name="ResultPlots"/><h3>Result Plots</h3>\n');
fprintf(fid,'<table>\n');
fprintf(fid,'<tr><th>Data Set</th>');
if(ProcessOutliers ~= 0)
  fprintf(fid,'<th>2D Outliers*</th><th>Hgt Outliers**</th>');
end

if(plotMap)
  fprintf(fid,'<th>Map</th>');
end
fprintf(fid,'<th>ENU</th><th>EN</th><th>CA</th><th>CAU</th>');
fprintf(fid,'<th>Vel</th>');
%fprintf(fid,'<th>East Vel</th><th>North Vel</th><th>Up Vel</th>');

for i=1:k
 if (~isnan(Yield(i)) && dataRate(i) > 0)
  printableName = strrep(InputDataName{i}{2},'\','');
  StatAxis{i} = printableName;

  fprintf(fid,'<tr>\n');

  fprintf(fid,['<td>',printableName,'</td>\n']);

  if(ProcessOutliers ~= 0)
    Str = [Section,'_',printableName,'_2D_Outliers.kmz'];
    fprintf(fid,['<td><a href="',Str,'">KMZ</a></td>\n']);
    Str = [Section,'_',printableName,'_Hgt_Outliers.kmz'];
    fprintf(fid,['<td><a href="',Str,'">KMZ</a></td>\n']);
  end

  if(plotMap)
    StrPNG = [Section,'_',printableName,'_EN_MapView.png'];
    fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>']);
    if(saveFig) 
      StrFIG = [Section,'_',printableName,'_EN_MapView.fig'];
      fprintf(fid,['&nbsp;<a href="',StrFIG,'"><font size="1"><sup>FIG</sup></font></a>']);
    end
    fprintf(fid,['</td>']);
  end

  StrPNG = [Section,'_',printableName,'_ENU.png'];
  fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>']);
  if(saveFig) 
    StrFIG = [Section,'_',printableName,'_ENU.fig'];
    fprintf(fid,['&nbsp;<a href="',StrFIG,'"><font size="1"><sup>FIG</sup></font></a>']);
  end
  fprintf(fid,['</td>']);

  StrPNG = [Section,'_',printableName,'_EN.png'];
  fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>']);
  if(saveFig) 
    StrFIG = [Section,'_',printableName,'_EN.fig'];
    fprintf(fid,['&nbsp;<a href="',StrFIG,'"><font size="1"><sup>FIG</sup></font></a>']);
  end
  fprintf(fid,['</td>']);

  if(~isnan(CrossStats(i,1)))
    StrPNG = [Section,'_',printableName,'_CA.png'];
    fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>']);
    if(saveFig) 
      StrFIG = [Section,'_',printableName,'_CA.fig'];
      fprintf(fid,['&nbsp;<a href="',StrFIG,'"><font size="1"><sup>FIG</sup></font></a>']);
    end
    fprintf(fid,['</td>']);
  else
    fprintf(fid,'<td>-</td>');
  end

  if(~isnan(CrossStats(i,1)))
    StrPNG = [Section,'_',printableName,'_CAU.png'];
    fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>']);
    if(saveFig) 
      StrFIG = [Section,'_',printableName,'_CAU.fig'];
      fprintf(fid,['&nbsp;<a href="',StrFIG,'"><font size="1"><sup>FIG</sup></font></a>']);
    end
    fprintf(fid,['</td>']);
  else
    fprintf(fid,'<td>-</td>');
  end
  
  StrPNG = [Section,'_',printableName,'_Vel.png'];
  fprintf(fid,['<td><a href="',StrPNG,'">PNG</a>']);
  if(saveFig) 
    StrFIG = [Section,'_',printableName,'_Vel.fig'];
    fprintf(fid,['&nbsp;<a href="',StrFIG,'"><font size="1"><sup>FIG</sup></font></a>']);
  end
  fprintf(fid,['</td>']);
  fprintf(fid,'</tr>\n');
 end
end

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

% add a link to the statistics output file
fprintf(fid,['<a href="',DiaryFileName,'">',DiaryFileName,'</a><br/>']);

if(ProcessOutliers ~= 0)
  fprintf(fid,'<h4>Outlier Definition:</h4>');
  if( (d2DOutlierThres > 0) && (d2DOutlierPer > 0) )
    fprintf(fid,'<h5>* 2D Outliers are points with the 2D position error greater than %.3f', d2DOutlierThres);
    fprintf(fid,'m or the worst %.3f',100-d2DOutlierPer);
    fprintf(fid,'%% versus truth</h5>\n');
  elseif( (d2DOutlierThres > 0) && (d2DOutlierPer < 0) )
    fprintf(fid,'<h5>* 2D Outliers are points with the 2D position error greater than %.3f', d2DOutlierThres);
    fprintf(fid,'m  versus truth</h5>\n');
  elseif( (d2DOutlierThres < 0) && (d2DOutlierPer > 0) )
    fprintf(fid,'<h5>* 2D Outliers are points with the 2D position in the worst %.3f',100-d2DOutlierPer);
    fprintf(fid,'%% versus truth</h5>\n');
  end
  
  if( (dHgtOutlierThres > 0) && (dHgtOutlierPer > 0) )
    fprintf(fid,'<h5>** Hgt Outliers are points with the Hgt position error greater than %.3f', dHgtOutlierThres);
    fprintf(fid,'m or the worst %.3f',100-dHgtOutlierPer);
    fprintf(fid,'%% versus truth</h5>\n');
  elseif( (dHgtOutlierThres > 0) && (dHgtOutlierPer < 0) )
    fprintf(fid,'<h5>** Hgt Outliers are points with the Hgt position error greater than %.3f', dHgtOutlierThres);
    fprintf(fid,'m  versus truth</h5>\n');
  elseif( (dHgtOutlierThres < 0) && (dHgtOutlierPer > 0) )
    fprintf(fid,'<h5>** Hgt Outliers are points with the Hgt position in the worst %.3f',100-dHgtOutlierPer);
    fprintf(fid,'%% versus truth</h5>\n');
  end
end
% fprintf(fid,'<hr/>\n');

%% Eric - Add AINS Section - Table and Plots
for i_table = 1:length(InputDataName)
    
  printableName = strrep(InputDataName{i_table}{2},'\','');
  fprintf(fid,['<a name="AINS"/><h3>',printableName,'</h3>\n']);

  % AINS Table
  % add a table of the results
  fprintf(fid,'<table>\n');
  fprintf(fid,'<tr><th>Errors</th>');
  fprintf(fid,'<th>Mean</th><th>STD</th><th>RMS</th>');
  fprintf(fid,'<th>50%%</th><th>68%%</th><th>95%%</th><th>99%%</th><th>100%%</th></tr>\n');

  clear StatAxis;

  if(ProcessAttitude == 1)
    pva_length = 3;
  else
    pva_length = 2; % No Attitude
  end
  
  % Pos, vel and att
  for i=1:pva_length*3
      
    fprintf(fid,'<tr>\n');  
      
    % Name
    if (i == 1)
      pva_name = 'Position - East [m]'; 
      pva_stats = EastStats;
    elseif (i == 2)
      pva_name = 'Position - North [m]'; 
      pva_stats = NorthStats;
    elseif (i == 3)
      pva_name = 'Position - Up [m]';  
      pva_stats = VertStats;
    elseif (i == 4)
      pva_name = 'Velocity - East [m/s]';
      pva_stats = EastVelStats;
    elseif (i == 5)
      pva_name = 'Velocity - North [m/s]';
      pva_stats = NorthVelStats;
    elseif (i == 6)
      pva_name = 'Velocity - Up  [m/s]';
      pva_stats = UpVelStats;
    elseif (i == 7)
      pva_name = 'Attitude - Roll [deg]';
      pva_stats = RollStats;
    elseif (i == 8)
      pva_name = 'Attitude - Pitch [deg]';
      pva_stats = PitchStats;
    elseif (i == 9)
      pva_name = 'Attitude - Heading [deg]';
      pva_stats = HeadingStats;
    end
          
    fprintf(fid,'<td>%s</td>\n',pva_name);
    
    % Mean, std and rms 
    if ((i <= 6) || (isnan(InputDataName{i_table}{20}) == 0) ...
     && (isnan(InputDataName{i_table}{21}) == 0) ...
     && (isnan(InputDataName{i_table}{22}) == 0) ...
     && (isnan(dRoll) == 0) )
      for j=1:3     
        % Mean, Sigma and RMS    
        fprintf(fid,'<td>%.3f</td>\n',mean_std_rms(i_table,i,j));
      end
      
      for j=1:5
        % 50%, 68%, 95%, 99%, 100%
        fprintf(fid,'<td>%.3f</td>\n',pva_stats(i_table,j));
      end
    else
      for j=1:3     
        % Mean, Sigma and RMS    
        fprintf(fid,'<td>%s</td>\n','-');
      end
      
      for j=1:5
        % 50%, 68%, 95%, 99%, 100%
        fprintf(fid,'<td>%s</td>\n','-');
      end
    end  
    
  
    fprintf(fid,'</tr>\n');
  end
  fprintf(fid,'</table>\n');
end

% AINS Plots
for i = 1:length(InputDataName)
  printableName = strrep(InputDataName{i}{2},'\','');

  fprintf(fid,['<a name="AINS"/><h3>',printableName,'</h3>\n']);

  if(plotMap)
    StrHTML = [Section,'_',printableName,'_EN_MapView.html'];
    StrPNG  = [Section,'_',printableName,'_EN_MapView.png'];
    fprintf(fid,['<a href="',StrHTML,'"><img src="',StrPNG,'" height="400"></a>']);
    fprintf(fid,'<figcaption>Click on the map for a dynamic Open Street View version</figcaption>');
  end

  % Position
  File = ['AINS_', Section,'_',printableName,'_ENU'];
  fprintf(fid,['<a href="',File,'.png"><img src="',File,'.png" height="400"/></a>']);
  if(saveFig)
    fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>&nbsp;']);
  end

  % Velocity
  File = ['AINS_', Section,'_',printableName,'_Vel'];
  fprintf(fid,['<a href="',File,'.png"><img src="',File,'.png" height="400"/></a>']);
  if(saveFig)
    fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>&nbsp;']);
  end

  % Attitude
  if(  (isnan(InputDataName{i}{20}) == 0) ...
    && (isnan(InputDataName{i}{21}) == 0) ...
    && (isnan(InputDataName{i}{22}) == 0) ...
    && (isnan(dRoll) == 0) ...
    && (ProcessAttitude == 1) )

    File = ['AINS_', Section,'_',printableName,'_Att'];
    fprintf(fid,['<a href="',File,'.png"><img src="',File,'.png" height="400"/></a>']);
    if(saveFig)
      fprintf(fid,['&nbsp;<a href="',File,'.fig"><font size="1"><sup>FIG</sup></font></a>&nbsp;']);
    end

  end
end

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

%% Eric - End of Add AINS Section - Table and Plots
first = 0;

for i=1:k
  if(isnan(InputDataName{i}{19}) == 0)
   if(first == 0) 
      fprintf(fid,'<a name="Yield"/><h2>Yield Info</h2>\n');
      first = 1;
    end
    fprintf(fid,'<table>\n');
    printableName = strrep(InputDataName{i}{2},'\','');
  
    fprintf(fid,'<tr><th colspan="2">%s</th></tr>',printableName);

    fprintf(fid,'<tr><th>Position Type</th><th>Yield [%%]</th></tr>');
    if( (abs(Yield(i) - 100) * 1e6) > 1)
      fprintf(fid,'<tr><td>No Position</td><td>%.3f</td>',100 - Yield(i));
    end

    [row col] = size(PosTypeArray{i});
    for type=1:row
      if(PosTypeArray{i}(type,2) > 0)
        fprintf(fid,'<tr><td>%s</td><td>%.3f</td>',PosStrings{PosTypeArray{i}(type,1) + 1}, 100* PosTypeArray{i}(type,2) * Yield(i)/100 );
      end
    end
    fprintf(fid,'</table>\n');
    fprintf(fid,'<br/>\n');
  end
end

if(first ~= 0) 
  fprintf(fid,'<hr/>\n');
end

fprintf(fid,'<a name="Stats"/><h2>Other Stats</h2>\n');

for s=1:MaxPlots
  % add header
  if(s == 1)
    fprintf(fid,'<a name="CrossTrack"/><h3>Cross-Track Stats:</h3>\n');
    fprintf(fid,['<a href="', Section, '_PenalizedCrossTrackCDF.png">Penalized CDF***</a>']);
    stats = CrossStats;
  elseif(s == 2)
    fprintf(fid,'<a name="AlongTrack"/><h3>Along-Track Stats:</h3>\n');
    fprintf(fid,['<a href="', Section, '_PenalizedAlongTrackCDF.png">Penalized CDF***</a>']);
    stats = AlongStats;
  elseif(s == 3)
    fprintf(fid,'<a name="2D"/><h3>2D Stats:</h3>\n');
    fprintf(fid,['<a href="', Section, '_Penalized2DCDF.png">Penalized CDF***</a>']);
    stats = TwoDStats;
  elseif(s == 4)
    fprintf(fid,'<a name="Hgt"/><h3>Vert Stats:</h3>\n');
    fprintf(fid,['<a href="', Section, '_PenalizedHgtCDF.png">Penalized CDF***</a>']);
    stats = VertStats;
  elseif(s == 5)
    fprintf(fid,'<a name="HorzVel"/><h3>Horizontal Velocity Stats:</h3>\n');
    stats = HorzVelStats;
  elseif(s == 6)
    fprintf(fid,'<a name="HgtVel"/><h3>Vertical Velocity Stats:</h3>\n');
    stats = UpVelStats;
  elseif(s == 7)
    fprintf(fid,'<a name="Roll"/><h3>Roll Stats:</h3>\n');
    stats = RollStats;
  elseif(s == 8)
    fprintf(fid,'<a name="Pitch"/><h3>Pitch Stats:</h3>\n');
    stats = PitchStats;
  elseif(s == 9)
    fprintf(fid,'<a name="Heading"/><h3>Heading Stats:</h3>\n');
    stats = HeadingStats;
  end
  
  % add a table of the results
  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><th>Primary [%%]**</th><th>Err Map</th>\n');

  if( (s==3) & (exist('hist3')==2))
    fprintf(fid,'<th>E/N Heat Map</th>\n');
    fprintf(fid,'<th>2D vs Est Err </th>\n');
  end
  fprintf(fid,'</tr>\n');

  RowCol = size(stats);
  clear StatAxis;
  for i=1:k
   if (~isnan(Yield(i)) && dataRate(i) > 0 && ~isnan(stats(i,1)))
    printableName = strrep(InputDataName{i}{2},'\','');
    StatAxis{i} = printableName;

    fprintf(fid,'<tr>\n');

    % Provide the name for each row of the table ... for some rows link
    % to a plot
    if( (s == 1) || (s==2) )
      Str = [Section,'_',printableName,'_CAU'];
      fprintf(fid,['<td><a href="',Str,'.png">',printableName,'</a>']);
      if(saveFig)
        fprintf(fid,['&nbsp;<a href="',Str,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
      end
      fprintf(fid,['</td>\n']);
    elseif(s == 3)
      Str = [Section,'_',printableName,'_ENErrVsEst'];
      fprintf(fid,['<td><a href="',Str,'.png">',printableName,'</a>']);
      if(saveFig)
        fprintf(fid,['&nbsp;<a href="',Str,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
      end
      fprintf(fid,['</td>\n']);
    elseif(s == 4)
      Str = [Section,'_',printableName,'_ENU'];
      fprintf(fid,['<td><a href="',Str,'.png">',printableName,'</a>']);
      if(saveFig)
        fprintf(fid,['&nbsp;<a href="',Str,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
      end
      fprintf(fid,['</td>\n']);
    %elseif(s == 5)
        % TODO: generate & link to a horizontal velocity time series plot
    %elseif(s == 6)
        % TODO: generate & link to a vertical velocity time series plot
    elseif(s == 7)
      Str = [Section,'_',printableName,'_RollErr'];
      fprintf(fid,['<td><a href="',Str,'.png">',printableName,'</a>']);
      if(saveFig)
        fprintf(fid,['&nbsp;<a href="',Str,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
      end
      fprintf(fid,['</td>\n']);
    elseif(s == 8)
      Str = [Section,'_',printableName,'_PitchErr'];
      fprintf(fid,['<td><a href="',Str,'.png">',printableName,'</a>']);
      if(saveFig)
        fprintf(fid,['&nbsp;<a href="',Str,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
      end
      fprintf(fid,['</td>\n']);
    elseif(s == 9)
      Str = [Section,'_',printableName,'_HeadingErr'];
      fprintf(fid,['<td><a href="',Str,'.png">',printableName,'</a>']);
      if(saveFig)
        fprintf(fid,['&nbsp;<a href="',Str,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
      end
      fprintf(fid,['</td>\n']);
    else
      fprintf(fid,'<td>%s</td>\n',printableName);
    end

    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

    if( isnan(PrimarySoln(i)) == 0)
      fprintf(fid,'<td>%s</td>\n',PosStrings{PrimarySoln(i) + 1});
      fprintf(fid,'<td>%3f</td>\n',PrimarySolnPercent(i));
    else
      fprintf(fid,'<td>-</td><td>-</td>\n');
    end
    
    if (s == 1)
        Str = [Section,'_',printableName,'_CrossTrackMap'];
    elseif s==2
        Str = [Section,'_',printableName,'_AlongTrackMap'];
    elseif s==3
        Str = [Section,'_',printableName,'_ENErrMap'];
    elseif s==4
        Str = [Section,'_',printableName,'_HeightErrMap'];
    elseif s==5
        Str = [Section,'_',printableName,'_HorzVelErrMap'];
    elseif s==6
        Str = [Section,'_',printableName,'_VertVelErrMap'];
    elseif s==7
        Str = [Section,'_',printableName,'_RollErrMap'];
    elseif s==8
        Str = [Section,'_',printableName,'_PitchErrMap'];
    elseif s==9
        Str = [Section,'_',printableName,'_HeadErrMap'];
    end;
    fprintf(fid,['<td><a href="',Str,'.png">PNG</a>']);
    if(saveFig)
      fprintf(fid,['&nbsp;<a href="',Str,'.fig"><font size="1"><sup>FIG</sup></font></a>']);
    end
    fprintf(fid,['</td>\n']);
    
    if( (s==3) & (exist('hist3')==2))
      Str = [Section,'_',printableName,'_EN-ENerrPlot.png'];
      fprintf(fid,['<td><a href="',Str,'">PNG</a></td>\n']);
      Str = [Section,'_',printableName,'-errPlot.png'];
      fprintf(fid,['<td><a href="',Str,'">PNG</a></td>\n']);
    end 

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

  % add CDF
  if(s == 1)
    File = [Section,'_CrossTrackCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
    File = [Section,'_CrossTrackCDF_95pcnt.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  elseif(s == 2)
    File = [Section,'_AlongTrackCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
    File = [Section,'_AlongTrackCDF_95pcnt.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  elseif(s == 3)
    File = [Section,'_2DCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
    File = [Section,'_2DCDF_95pcnt.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  elseif(s == 4)
    File = [Section,'_HgtCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  elseif(s == 5)
    File = [Section,'_HorzVelCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  elseif(s == 6)
    File = [Section,'_HgtVelCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  elseif(s == 7)
    File = [Section,'_RollCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  elseif(s == 8)
    File = [Section,'_PitchCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  elseif(s == 9)
    File = [Section,'_HeadingCDF.png'];
    fprintf(fid,['<a href="',File,'"><img src="',File,'" height="400"/></a>']);
  end
  
  % add bar graphs
  [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
    bar(stats(:,1:4));
    A = axis;
    if(s == 1)
      title([Section,' Cross-Track Stats']);
      name =  [Section,'_CrossTrack_summary.png'];
      A(4) = HorzAxisMax;
      ylabel('Error [m]');
    elseif(s == 2)
      title([Section,' Along-Track Stats']);
      name = [Section,'_AlongTrack_summary.png'];
      A(4) = HorzAxisMax;
      ylabel('Error [m]');
    elseif(s==3)
      title([Section,' 2D Stats']);
      name = [Section,'_2D_summary.png'];
      A(4) = HorzAxisMax;
      ylabel('Error [m]');
    elseif(s==4)
      title([Section,' Vert Stats']);
      name = [Section,'_Vertical_summary.png'];
      A(4) = VertAxisMax;
      ylabel('Error [m]');
    elseif(s==5)
      title([Section,' Horizontal Velocity Stats']);
      name = [Section,'_HorzVel_summary.png'];
      A(4) = VelAxisMax;
      ylabel('Error [m/s]');
    elseif(s==6)
      title([Section,' Vertical Velocity Stats']);
      name = [Section,'_VertVel_summary.png'];
      A(4) = VelAxisMax;
      ylabel('Error [m/s]');
    elseif(s==7)
      title([Section,' Roll Stats']);
      name = [Section,'_Roll_summary.png'];
      A(4) = RPAxisMax;
      ylabel('Error [deg]');
    elseif(s==8)
      title([Section,' Pitch Stats']);
      name = [Section,'_Pitch_summary.png'];
      A(4) = RPAxisMax;
      ylabel('Error [deg]');
    elseif(s==9)
      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,['&nbsp;<a href="',name,'"><img src="',name,'" height="400"/></a>']);
  end
  
  fprintf(fid,'<br/>');
  
  fprintf(fid,'<hr/>\n');
end

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. The primary percent is the percentage the primary ' ...
    'solution occurs 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 and the primary percentage would be 50%% (400 out of ' ...
    '800 logged).</h5>\n']);

fprintf(fid,['<h5>*** - The Penalized CDF takes into account the yield. ' ...
             'For example if a platform only has 80%% yield the CDF will ' ...
             'never go beyond 80%']);
             
fprintf(fid,'</body>\n');
fprintf(fid,'</html>\n');

fclose(fid);

% open the HTML file in a browser
if (exist(HTML_name,'file')==2)
    status=web(HTML_name,'-browser');
    if (status==1)
        fprintf('"%s" could not be opened!\n',HTML_name);
    end;
else
    fprintf('HTML file does not exist!  Skipping open.\n');
end;
