function pos(PREFIX,PosType,Compiled,HorzFilt,VertFilt,AxisType,FileType,EpochRate,LocalOffset);

%-----------------------------------------------------------------------
% $Id: pos.m,v 1.10 2002/05/22 17:13:01 stuart Exp $
%
% Takes a parsed output from viewdat or a parsed version of the RTK-PC
% output. The data must be the format:
%
% FileType = 0 (default)
% Time, dE, dN, dU, #SVs, PDOP, Pos-type
%
% FileType = 1
% Time, dE, dN, dU, #SVs, PDOP, Pos-type, RMS, sigE, sigN, sigU
%
% Produces a set of graphs which are saved in PNG format. Also creates a
% file stats.csv with the position domain statistics. This can be easily
% parsed to provide HTML etc.
%
% PREFIX   - String prefixed to the name of the output .png files
% PosType  - The data is filter to only analyze this data type.
% Compiled - (=0) Called by Matlab; (=1) Standalone application.
%            If the compiled flag is ommitted it will assume the
%            application is not compiled
% HorzFilt 
% VertFilt - Both in Metres, allow the user to produce plots of data
%            filtering out positions above the threshold. If both are
%            set to 1000 it is ignored (this allows us to use the
%            additional options without filtering the data).
% AxisType - (=0) time axis is in seconds, (=1) time axis is in HH:MM
% FileType - (=0) parsed .dat file, (=1) parsed rtk-pc file. The fields
%            are given above.
% EpochRate - the epoch rate in seconds, used to provide pen-up/pen-down
%             functionallity when there are gaps greated then 4 x
%             epoch rate.
% LocalOffset - Number of hours the time is offset from UTC. This is an
%               optional parameter and only used when AxisType=1.
% 
% For the compiled StandAlone application, these passed parameters are
% strings, and converted to numbers prior to first use.
% To 'simulate' running the compiled version from Matlab, the parameters
% must be passed as strings, for example:
%
%   pos('Prefix',PosType,'1' )  
%
% The capabilities have been expanded to better support fixed RTK data.
% If the position type is 76 (fixed RTK) the bin size for the histograms
% is reduced to 1mm and the x-axis of the graphs scaled to the range of
% the data +/-1mm. For non-fixed RTK data the bin step size is 5cm with
% the x-axis data range set to +/- 1m.
%
% Copyright Trimble Navigation 2000-2002
% Stuart Riley 
% 1st Oct 2000
% 
%-----------------------------------------------------------------------

%PREFIX = 'Test';
%PosType = 76;
%Compiled = 0;
%HorzFilt = 1000;
%VertFilt = 1000;
%AxisType = 1;
%FileType = 1;
%EpochRate = 1;
%LocalOffset = -7;
%nargin = 9;

FALSE = 0;
TRUE  = 1;
ASCII_ZERO = 48 ;
ASCII_ONE  = 49 ;

%
% First must check if using compiled StandAlone version. This passes
% parameters as strings (parameters from Matlab are numbers).
% Note: if the Compiled version uses no parameters, this is an error
% that cannot be detected (defaults to "not compiled").
%
if( nargin > 2 )
  if( Compiled == ASCII_ZERO )
    % Compiled application passing "not compiled". This is an error.
    fprintf( 'ERROR! Compiled Application must set first parameter to 1\n' ) ;
    return ;
  elseif( Compiled == ASCII_ONE )
    Compiled = TRUE ;
    % Compiled Version passes parameters as strings.
  else
    % Running Matlab and passing 0 or 1 (as a number).
    % To 'simulate' running the Compiled version from Matlab, the
    % parameter 'Compiled' must be passed as the string '1'. See header.
    Compiled = FALSE ;
  end
else
  % Default to not compiled if there are 2 or fewer input parameters
  Compiled = FALSE;
end


if(nargin==0)
  PREFIX='';
  PosType = 4;
  Filter= FALSE;
  AxisType = 0;
  FileType = 0;
  EpochRate = NaN;
  LocalOffset = NaN;
elseif(nargin==1)
  PosType = 4;
  Filter= FALSE;
  AxisType = 0;
  FileType = 0;
  EpochRate = NaN;
  LocalOffset = NaN;
elseif(nargin <= 4)
  Filter= FALSE;
  AxisType = 0;
  FileType = 0;
  EpochRate = NaN;
  LocalOffset = NaN;
  if(Compiled == TRUE)
    PosType = str2num( PosType ) ;
  end
elseif(nargin == 5)
  Filter= TRUE;
  AxisType = 0;
  FileType = 0;
  EpochRate = NaN;
  LocalOffset = NaN;
  if(Compiled == TRUE)
    PosType = str2num( PosType ) ;
  end
elseif(nargin == 6)
  Filter= TRUE;
  FileType = 0;
  EpochRate = NaN;
  LocalOffset = NaN;
  if(Compiled == TRUE)
    PosType  = str2num( PosType ) ;
    AxisType = str2num( AxisType ) ;
  end
elseif(nargin == 7)
  Filter= TRUE;
  EpochRate = NaN;
  LocalOffset = NaN;
  if(Compiled == TRUE)
    PosType = str2num( PosType ) ;
    AxisType = str2num( AxisType ) ;
    FileType = str2num( FileType ) ;
  end
elseif(nargin == 8)
  Filter= TRUE;
  LocalOffset = NaN;
  if(Compiled == TRUE)
    PosType = str2num( PosType ) ;
    AxisType = str2num( AxisType ) ;
    FileType = str2num( FileType ) ;
    EpochRate = str2num(EpochRate);
  end
else
  Filter= TRUE;
  if(Compiled == TRUE)
    PosType = str2num( PosType ) ;
    AxisType = str2num( AxisType ) ;
    FileType = str2num( FileType ) ;
    EpochRate = str2num(EpochRate);
    LocalOffset = str2num(LocalOffset);
  end
end  

if(Filter == TRUE)
  if(Compiled == TRUE)
    HorzFilt = str2num(HorzFilt);
    VertFilt = str2num(VertFilt);
  end

  % Turn off the filtering if Vert=Horz=1000
  if( (HorzFilt == 1000) & (VertFilt == 1000) )
    Filter = FALSE;
  end
end

% For test/development provide a way to stop Matlab from saving the
% plots to PNG files... this speeds up the script a little, especially
% on my prehistoric clock-work PC ;-)
savePlots = TRUE;
% savePlots = FALSE;

OutputGraphicsType = '-dpng';
FileExtension = 'png';

if( Compiled == FALSE )
  load test.enu
  fileError = FALSE;
  %if(length(test) == 0)
  %  fileError = TRUE;
  %else
  %  fileError = FALSE;
  %end
else
  %
  % The compiled version of Matlab does not support 
  % the "load" command. The code below solves this problem.
  %
  fid = fopen( 'test.enu', 'rt' ) ;
  if(fid == -1)
    fileError = TRUE;
  else
    fileError = FALSE;
    if(FileType  == 0)
      test = fscanf( fid, '%f', [7, inf] ) ;
    else
      test = fscanf( fid, '%f', [11, inf] ) ;
    end
    fclose( fid ) ;
    test = test' ;  % Transpose as fscanf() fills by column.
  end
end

disp('$Id: pos.m,v 1.10 2002/05/22 17:13:01 stuart Exp $');

if(fileError == TRUE)
  disp('Error opening data file test.enu');
else
  % If the data comes from a .dat file it will have week wrap-overs.
  % Detect for and correct for this.... assume that there is only one
  % jump in the data, ie could not span two weeks!
  DataDiff = diff(test(:,1));
  i = find(DataDiff < -3.5*86400); % all jumps backwards greater than 3.5 days
  if(length(i) > 0)
    test( (i(1)+1):length(test),1) = test( (i(1)+1):length(test),1) + 7*86400;
  end
  
  dataLength     = length(test);
  expectedLength = test(dataLength,1) - test(1,1);
  
  if(PosType == 76)
    %----------------------------------------------------------------
    %  RTK fixed requested - in addition we want to sort the data and
    %                        compute number of fixed etc positions.
    %----------------------------------------------------------------
   
    % 3D - Autonomous
    j = find(test(:,7) == 4);
    Type4 = length(j);
    
    % 3D - Differential
    j = find(test(:,7) == 20);
    Type20 = length(j);
    
    % 3D - Float
    j = find(test(:,7) == 68);
    Type68 = length(j);
  
    % This needs to be last as we use this info next
    j = find(test(:,7) == PosType);
    Type76 = length(j);
  
  else
    j = find(test(:,7) == PosType);
  end                 
  
  %
  % Remove any positions that are not the requested type... yes this could
  % probably be done more efficiently!
  %
  i = find(test(:,7) ~= PosType);
    
  if(Filter == TRUE)
    % Good inits (if this is RTK)
    k = find(      ( (test(j,2).^2+test(j,3).^2) < (HorzFilt*HorzFilt) ) ...
                 & (  test(j,4) < VertFilt) );
    NumGood = length(k);
      
    % Bad inits (if this is RTK)
    m = find(      ( (test(j,2).^2+test(j,3).^2) >= (HorzFilt*HorzFilt) ) ...
                 | (  test(j,4) >= VertFilt) );
  
    testStats = test(j(k),:);
    test(j(m),2) = NaN;
    test(j(m),3) = NaN;
    test(j(m),4) = NaN;
    test(j(m),5) = NaN;
    test(j(m),6) = NaN;
    test(j(m),7) = NaN;
  else
    NumGood = NaN;
    testStats = test(j,:);
  end
    
  len = length(testStats);
    
  test(i,2) = NaN;
  test(i,3) = NaN;
  test(i,4) = NaN;
  test(i,5) = NaN;
  test(i,6) = NaN;
  test(i,7) = NaN;
    
  clear i;
  
  
  % Create a CSV format file with the statistical results.
  ftxt = fopen('stats.csv','w');
  
  if( (dataLength > 0 ) & (length(j) == 0) )
    % No data of the requested type is available.
    
    fprintf(ftxt,' ,1-Sigma,RMS,Mean,50%%,68%%,95%%,Min,Max\n');
    fprintf(ftxt,'East,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN\n');
    fprintf(ftxt,'North,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN\n');
    fprintf(ftxt,'Height,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN\n');
    fprintf(ftxt,' ,RMS,CEP (50%%),95%%,Max\n');
    fprintf(ftxt,'2D,NaN,NaN,NaN,NaN\n');
    fprintf(ftxt,'3D,NaN,NaN,NaN,NaN\n');
  
    if( (PosType == 76) & (Filter == TRUE) )
     % Provide some info on the position types
     fprintf(ftxt,'Data Types, Num Expected, Num Epochs, Num 3D-Fixed (Good), Num 3D-Fixed (Bad), Num 3D-Float, Num 3D-DGPS, Num 3D-Autonomous\n');
     fprintf(ftxt,' ,%d,%d,%d,%d,%d,%d,%d\n', ...
                expectedLength, ...
                dataLength, ...
                NumGood, ...
                Type76 - NumGood, ...
                Type68, ...
                Type20, ...
                Type4);
    end            
  else
    
    
    %
    % Calculate the RMS
    % Setting the Truth to be the mean will give us the
    % standard deviation.
    % Truth_E= mean(test(:,2));
    % Truth_N= mean(test(:,3));
    % Truth_U= mean(test(:,4));
    %
    
    Truth_E= 0.0;
    Truth_N= 0.0;
    Truth_U= 0.0;
    
    SumAll_E = 0.0;
    SumAll_N = 0.0;
    SumAll_U = 0.0;
    
    for i=1:length(testStats)
      SumAll_E = SumAll_E + (testStats(i,2) - Truth_E)^2;
      SumAll_N = SumAll_N + (testStats(i,3) - Truth_N)^2;
      SumAll_U = SumAll_U + (testStats(i,4) - Truth_U)^2;
    end
    
    RMS_E = sqrt( SumAll_E / (length(testStats) - 1));
    RMS_N = sqrt( SumAll_N / (length(testStats) - 1));
    RMS_U = sqrt( SumAll_U / (length(testStats) - 1));
    
    TRUE      = 1;
    FALSE     = 0;
    
    figure
    subplot(2,1,1);
    if(PosType == 76) % fixed RTK
      range = -1:0.001:1;
    else  
      range = -100:0.05:100;
    end
    lastValue  = length(range);
    startValue = floor(length(range)/2) + 1;
    
    %
    % Process the East data
    %
    
    N=hist(testStats(:,2),range);
    bar(range,100*N./len,1);
    xlabel('East Position Error [m]');
    ylabel('Number of Epochs [%]');
    
    Sigma_E =  std(testStats(:,2));
    Mean_E  = mean(testStats(:,2));
    title([PREFIX,': East - Sigma = ',num2str(Sigma_E,2), ...
          'm; RMS = ',num2str(RMS_E,2), ...
          'm; Mean Error = ',num2str(Mean_E,2),...
          'm'], ...
          'Interpreter','none');
    A = axis;
    
    % Extract some percentiles...
    data = sort(abs(testStats(:,2)));
    len  = length(testStats(:,2));
    
    E_50 = data(floor(0.50 * len));
    E_68 = data(floor(0.68 * len));
    E_95 = data(floor(0.95 * len));
    
    
    fprintf(ftxt,' ,1-Sigma,RMS,Mean,50%%,68%%,95%%,Min,Max\n');
    fprintf(ftxt,'East,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f\n', ...
                   Sigma_E, ...
                   RMS_E, ...
                   Mean_E, ...
                   E_50,E_68,E_95, ...
                   min(testStats(:,2)), ...
                   max(testStats(:,2)));
    
    % 
    % Find the first and last zero so we can figure out the best axis.
    %
    S = length(N);
    START = TRUE;
    
    for i=1:S
      if(N(i) ~= 0)
        if(START == TRUE)
          FirstNonZero = i;
          START = FALSE;
        end
        LastZero = i;
      end    
    end
    
    if(PosType == 76)
      % Closest mm for fixed RTK
      MinX = floor(range(FirstNonZero)*1000)/1000;
      MaxX = ceil(range(LastZero)*1000)/1000;
    else
      % Closest m for everything else
      MinX = floor(range(FirstNonZero));
      MaxX = ceil(range(LastZero));
    end
    
    A(1) = MinX;
    A(2) = MaxX;
    axis(A);
    grid;
    
    
    subplot(2,1,2)
    N=hist(abs(testStats(:,2)),range(startValue:lastValue));
    X = cumsum(N);
    plot(range(startValue:lastValue),X.*100/length(testStats));
    i = find(X < length(testStats));
    if(length(i) < (lastValue - startValue) )
      fullScale = range(startValue+length(i)+1);
    else
      fullScale = range(lastValue);
    end
    grid
    axis([0 fullScale 0 100]);
    xlabel('East Position Error [m]');
    ylabel('Cumulative Proability of Absolute Error [%]');
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_E_Dist.',FileExtension]);
    end
    
    %
    % Process the North Data
    %
    
    
    figure
    subplot(2,1,1);
    N=hist(testStats(:,3),range);
    bar(range,100*N./len,1);
    xlabel('North Position Error [m]');
    ylabel('Number of Epochs [%]');
    
    Sigma_N =  std(testStats(:,3));
    Mean_N  = mean(testStats(:,3));
    title([PREFIX,': North - Sigma = ',num2str(Sigma_N,2), ...
          'm; RMS = ',num2str(RMS_N,2), ...
          'm; Mean Error = ',num2str(Mean_N,2),...
          'm'], ...
          'Interpreter','none');
    A = axis;
    
    % Extract some percentiles...
    data = sort(abs(testStats(:,3)));
    len  = length(testStats(:,3));
    
    N_50 = data(floor(0.50 * len));
    N_68 = data(floor(0.68 * len));
    N_95 = data(floor(0.95 * len));
    
    
    fprintf(ftxt,'North,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f\n', ...
                   Sigma_N, ...
                   RMS_N, ...
                   Mean_N, ...
                   N_50,N_68,N_95, ...
                   min(testStats(:,3)), ...
                   max(testStats(:,3)));
    
    % 
    % Find the first and last zero so we can figure out the best axis.
    %
    S = length(N);
    START = TRUE;
    
    for i=1:S
      if(N(i) ~= 0)
        if(START == TRUE)
          FirstNonZero = i;
          START = FALSE;
        end
        LastZero = i;
      end    
    end
    
    if(PosType == 76)
      % Closest mm for fixed RTK
      MinX = floor(range(FirstNonZero)*1000)/1000;
      MaxX = ceil(range(LastZero)*1000)/1000;
    else
      % Closest m for everything else
      MinX = floor(range(FirstNonZero));
      MaxX = ceil(range(LastZero));
    end
    A(1) = MinX;
    A(2) = MaxX;
    axis(A);
    grid
    
    
    subplot(2,1,2)
    N=hist(abs(testStats(:,3)),range(startValue:lastValue));
    X = cumsum(N);
    plot(range(startValue:lastValue),X.*100/length(testStats));
    i = find(X < length(testStats));
    if(length(i) < (lastValue - startValue) )
      fullScale = range(startValue+length(i)+1);
    else
      fullScale = range(lastValue);
    end
    grid
    axis([0 fullScale 0 100]);
    xlabel('North Position Error [m]');
    ylabel('Cumulative Proability of Absolute Error [%]');
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_N_Dist.',FileExtension]);
    end
    
    %
    % Process the Height Data
    %
    
    figure
    subplot(2,1,1);
    N=hist(testStats(:,4),range);
    bar(range,100*N./len,1);
    xlabel('Height Position Error [m]');
    ylabel('Number of Epochs [%]');
    
    Sigma_U =  std(testStats(:,4));
    Mean_U  = mean(testStats(:,4));
    title([PREFIX,': Height - Sigma = ',num2str(Sigma_U,2), ...
          'm; RMS = ',num2str(RMS_U,2), ...
          'm; Mean Error = ',num2str(Mean_U,2),...
          'm'], ...
          'Interpreter','none');
    A = axis;
    
    % Extract some percentiles...
    data = sort(abs(testStats(:,4)));
    len  = length(testStats(:,4));
    
    U_50 = data(floor(0.50 * len));
    U_68 = data(floor(0.68 * len));
    U_95 = data(floor(0.95 * len));
    
    fprintf(ftxt,'Height,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f\n', ...
                   Sigma_U, ...
                   RMS_U, ...
                   Mean_U, ...
                   U_50,U_68,U_95, ...
                   min(testStats(:,4)), ...
                   max(testStats(:,4)));
    
    % 
    % Find the first and last zero so we can figure out the best axis.
    %
    S = length(N);
    START = TRUE;
    
    for i=1:S
      if(N(i) ~= 0)
        if(START == TRUE)
          FirstNonZero = i;
          START = FALSE;
        end
        LastZero = i;
      end    
    end
    
    if(PosType == 76)
      % Closest mm for fixed RTK
      MinX = floor(range(FirstNonZero)*1000)/1000;
      MaxX = ceil(range(LastZero)*1000)/1000;
    else
      % Closest m for everything else
      MinX = floor(range(FirstNonZero));
      MaxX = ceil(range(LastZero));
    end
    A(1) = MinX;
    A(2) = MaxX;
    axis(A);
    grid
    
    
    subplot(2,1,2)
    N=hist(abs(testStats(:,2)),range(startValue:lastValue));
    X = cumsum(N);
    plot(range(startValue:lastValue),X.*100/length(testStats));
    i = find(X < length(testStats));
    if(length(i) < (lastValue - startValue) )
      fullScale = range(startValue+length(i)+1);
    else
      fullScale = range(lastValue);
    end
    grid
    axis([0 fullScale 0 100]);
    xlabel('Height Position Error [m]');
    ylabel('Cumulative Proability of Absolute Error [%]');
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_U_Dist.',FileExtension]);
    end
    
    %
    % Process the Horizontal Data
    %
    
    figure
    subplot(2,1,1);
    range = 0:0.05:100;
    if(PosType == 76) % fixed RTK
      range = 0:0.001:1;
    else  
      range = 0:0.05:100;
    end
    
    HorzErr = sqrt(testStats(:,2).^2 + testStats(:,3).^2); 
    N=hist(HorzErr,range);
    bar(range,100*N./len,1);
    xlabel('Horizontal Position Error [m]');
    ylabel('Number of Epochs [%]');
    A = axis;
    
    
    % 
    % Find the first and last zero so we can figure out the best axis.
    %
    S = length(N);
    START = TRUE;
    
    for i=1:S
      if(N(i) ~= 0)
        if(START == TRUE)
          FirstNonZero = i;
          START = FALSE;
        end
        LastZero = i;
      end    
    end
    
    if(PosType == 76)
      % Closest mm for fixed RTK
      MinX = floor(range(FirstNonZero)*1000)/1000;
      MaxX = ceil(range(LastZero)*1000)/1000;
    else
      % Closest m for everything else
      MinX = floor(range(FirstNonZero));
      MaxX = ceil(range(LastZero));
    end
    A(1) = MinX;
    A(2) = MaxX;
    axis(A);
    grid
    
  
    subplot(2,1,2)
    X = cumsum(N);
    plot(range,X.*100/length(testStats));
    grid
    axis([A(1:2) 0 100]);
    xlabel('Horizontal Position Error [m]');
    ylabel('Cumulative Proability [%]');
    
    %
    % Find the closest point to 50%
    % 
    i = find (X/length(testStats) < .5);
    i_max = find( i == max(i));
    
    %
    % Now interpolate to provide a better estimate of
    % the 95 percentile.
    % 
    
    if(i_max ~= size(range))
      % We are not on the last element therefore we can interpolate
      Y_horz(1) = X(i_max)/length(testStats);
      Y_horz(2) = X(i_max+1)/length(testStats);
      X_horz(1) = range(i_max);
      X_horz(2) = range(i_max+1);
      m         = (Y_horz(2) - Y_horz(1)) / (X_horz(2) - X_horz(1));
      c         = Y_horz(1) - m*X_horz(1);
      Horz_50   = (0.5 - c)/m;
    else
      Horz_50 = range(i_max);
    end
    
    
    %
    % Find the closest point to 95%
    % 
    i = find (X/length(testStats) < .95);
    i_max = find( i == max(i));
    
    %
    % Now interpolate to provide a better estimate of
    % the 95 percentile.
    % 
    
    if(i_max ~= size(range))
      % We are not on the last element therefore we can interpolate
      Y_horz(1) = X(i_max)/length(testStats);
      Y_horz(2) = X(i_max+1)/length(testStats);
      X_horz(1) = range(i_max);
      X_horz(2) = range(i_max+1);
      m         = (Y_horz(2) - Y_horz(1)) / (X_horz(2) - X_horz(1));
      c         = Y_horz(1) - m*X_horz(1);
      Horz_95   = (0.95 - c)/m;
    else
      Horz_95 = range(i_max);
    end
    
    Horz_RMS = sqrt(sum(HorzErr.^2) / length(HorzErr)); 
    
    subplot(2,1,1)
    title([PREFIX,': Horizontal Error: 95%=',num2str(Horz_95,'%.3f'),'m'],'Interpreter','none');
    fprintf(ftxt,' ,RMS,CEP (50%%),95%%,Max\n');
    fprintf(ftxt,'2D,%.3f,%.3f,%.3f,%.3f\n', ...
                   Horz_RMS,...
                   Horz_50,...
                   Horz_95,...
                   max(HorzErr));
    
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_HORZ_Dist.',FileExtension]);
    end
    
    %
    % Now plot in the time domain
    %
    
    % 
    % First find gaps in the data and set NaN's so we can get pen-up /
    % pen-down plotting around the gaps.
    %
    if(isnan(EpochRate))
      % don't have a defined epoch rate so we can't easily do anything...
      % yes we could try to estimate to EpochRate but there would always
      % be a corner case....
      %
      % Therefore do not apply any filtering, i.e. don't look for gaps.
      out = test;
    else
      ii = diff(test(:,1));
      % To prevent a gap of a couple of epochs (which happens a lot with
      % some TRS installs) triggering a pen-up / pen-down, the gap must be
      % greater than 4 epochs.
      i  = find(ii > EpochRate*4);
  
      if(length(i) > 0)
        % Complete the last segment so we get all data from the input
        i(length(i)+1) = length(test);
        out = test(1:i(1),:);
        outIndex = i(1);  % last element written into out
        for j = 2:length(i)
          out(outIndex + 1,1) = out(outIndex,1);
          if(FileType  == 0)
            out(outIndex + 1,2:7) = ones(1,6).*NaN;
          else
            out(outIndex + 1,2:11) = ones(1,10).*NaN;
          end
          outIndex = outIndex + 1; % last element written into out
    
          % length of the current contiguous segment
          segmentLength = i(j) - i(j-1);
    
          out( (outIndex+1) : (outIndex+segmentLength),:) = test(i(j-1)+1:i(j),:);
          outIndex = outIndex + segmentLength;
        end
      else
        out = test;
      end
    end
      
    figure;
    subplot(3,1,1);
    hold;
  
    if(FileType == 0)
      % Compute the min/max y values so we can optimize the axis
      % scaling.
      maxE = max(test(:,2));
      minE = min(test(:,2));
      
      X= [testStats(1,1), testStats(length(testStats),1) ];
      Y= [std(testStats(:,2)) + mean(testStats(:,2)), ...
         std(testStats(:,2)) + mean(testStats(:,2))];
         
      if(Y(1) > maxE)
        maxE = Y(1);
      end

      % Plot the red sigma line first so the legend command picks it up.
      % However, replot it later as the data may cover it!
      plot(X,Y,'--r');
      plot(out(:,1),out(:,2));
  
      % Number greater than 3sigma + mean
      i = find(test(:,2) > 3*std(testStats(:,2)) + mean(testStats(:,2)));
      num = length(i);
  
      plot(X,Y,'--r');
      Y= [-std(testStats(:,2)) + mean(testStats(:,2)), ...
        -std(testStats(:,2)) + mean(testStats(:,2))];

      if(Y(1) < minE)
        minE = Y(1);
      end
      
      plot(X,Y,'--r');
      
      % Number less than -3sigma + mean
      i = find(test(:,2) < -3*std(testStats(:,2)) + mean(testStats(:,2)));
      num = num + length(i);
      % As a percentage
      conf = 100*(length(test(:,1)) - num) / length(test(:,1));
      legend('+/- 1-Sigma Confidence',0);
      title([PREFIX,': East/North/Height Error [East: ',num2str(conf,'%.3f'),'% measurements within +/- 3 Sigma]'],'Interpreter','none');
      
    else
      % Compute the min/max y values so we can optimize the axis
      % scaling.
      maxE = max(  [max(test(:,2)) ...
                    max(test(:,9)+ mean(testStats(:,2)))] );
      
      minE = min(  [min(test(:,2)) ...
                    min(-test(:,9)+ mean(testStats(:,2)))] );
      
      plot(out(:,1), out(:,9) + mean(testStats(:,2)),'--r');
      plot(out(:,1), out(:,2));
      plot(out(:,1), out(:,9) + mean(testStats(:,2)),'--r');
      plot(out(:,1), -out(:,9) + mean(testStats(:,2)),'--r');
  
      % Greater than the +ve threshold
      i = test(:,2) > ( test(:,9) + mean(testStats(:,2)) );
      j = find(i==1);
      num = length(j);
  
      % Greater than the -ve threshold
      i = test(:,2) < ( -test(:,9) + mean(testStats(:,2)) );
      j = find(i==1);
      num = num + length(j);
  
      % As a percentage
      conf = 100*(length(test(:,1)) - num) / length(test(:,1));
  
      legend('Reported 99% Confidence',0);
      title([PREFIX,': East/North/Height Error [East ',num2str(conf,'%.3f'),'% measurements within 99% Confidence]'],'Interpreter','none');
    end
  
    grid;
    ylabel('East Error [m]');
    
    subplot(3,1,2);
    plot(out(:,1),out(:,3));
    hold;
    if(FileType == 0)
      % Compute the min/max y values so we can optimize the axis
      % scaling.
      maxN = max(test(:,3));
      minN = min(test(:,3));
      
      Y= [std(testStats(:,3)) + mean(testStats(:,3)), ...
          std(testStats(:,3)) + mean(testStats(:,3))];
      
      if(Y(1) > maxN)
        maxN = Y(1);
      end

      plot(X,Y,'--r');
      Y= [-std(testStats(:,3)) + mean(testStats(:,3)), ...
          -std(testStats(:,3)) + mean(testStats(:,3))];

      if(Y(1) < minN)
        minN = Y(1);
      end
      
      plot(X,Y,'--r');
      
      % Number greater than 3sigma + mean
      i = find(test(:,3) > 3*std(testStats(:,3)) + mean(testStats(:,3)));
      num = length(i);
      % Number less than -3sigma + mean
      i = find(test(:,3) < -3*std(testStats(:,3)) + mean(testStats(:,3)));
      num = num + length(i);
      % As a percentage
      conf = 100*(length(test(:,1)) - num) / length(test(:,1));
      title(['[North: ',num2str(conf,'%.3f'),'% measurements within +/- 3 Sigma']);
    else
      % Compute the min/max y values so we can optimize the axis
      % scaling.
      maxN = max(  [max(test(:,3)) ...
                    max(test(:,10)+ mean(testStats(:,3)))] );
      
      minN = min(  [min(test(:,3)) ...
                    min(-test(:,10)+ mean(testStats(:,3)))] );
                    
      plot(out(:,1), out(:,10) + mean(testStats(:,3)),'--r');
      plot(out(:,1), -out(:,10) + mean(testStats(:,3)),'--r');
      
      % Greater than the +ve threshold
      i = test(:,3) > ( test(:,10) + mean(testStats(:,3)) );
      j = find(i==1);
      num = length(j);
  
      % Greater than the -ve threshold
      i = test(:,3) < ( -test(:,10) + mean(testStats(:,3)) );
      j = find(i==1);
      num = num + length(j);
  
      % As a percentage
      conf = 100*(length(test(:,1)) - num) / length(test(:,1));
  
      title(['North ',num2str(conf,'%.3f'),'% measurements within 99% Confidence']);
    end
    grid;
    ylabel('North Error [m]');
    
    subplot(3,1,3);
    plot(out(:,1),out(:,4));
    hold;
    if(FileType == 0)
      % Compute the min/max y values so we can optimize the axis
      % scaling.
      maxU = max(test(:,4));
      minU = min(test(:,4));
      
      Y= [std(testStats(:,4)) + mean(testStats(:,4)), ...
          std(testStats(:,4)) + mean(testStats(:,4))];
      
      if(Y(1) > maxU)
        maxU = Y(1);
      end

      plot(X,Y,'--r');
      Y= [-std(testStats(:,4)) + mean(testStats(:,4)), ...
          -std(testStats(:,4)) + mean(testStats(:,4))];

      if(Y(1) < minU)
        minU = Y(1);
      end
      
      plot(X,Y,'--r');
      
      % Number greater than 3sigma + mean
      i = find(test(:,4) > 3*std(testStats(:,4)) + mean(testStats(:,4)));
      num = length(i);
      % Number less than -3sigma + mean
      i = find(test(:,4) < -3*std(testStats(:,4)) + mean(testStats(:,4)));
      num = num + length(i);
      % As a percentage
      conf = 100*(length(test(:,1)) - num) / length(test(:,1));
      title(['[Height: ',num2str(conf,'%.3f'),'% measurements within +/- 3 Sigma']);
    else
      % Compute the min/max y values so we can optimize the axis
      % scaling.
      maxU = max(  [max(test(:,4)) ...
                    max(test(:,11)+ mean(testStats(:,4)))] );
      
      minU = min(  [min(test(:,4)) ...
                    min(-test(:,11)+ mean(testStats(:,4)))] );
                    
      plot(out(:,1), out(:,11) + mean(testStats(:,4)),'--r');
      plot(out(:,1), -out(:,11) + mean(testStats(:,4)),'--r');
      
      % Greater than the +ve threshold
      i = test(:,4) > ( test(:,11) + mean(testStats(:,4)) );
      j = find(i==1);
      num = length(j);
  
      % Greater than the -ve threshold
      i = test(:,4) < ( -test(:,11) + mean(testStats(:,4)) );
      j = find(i==1);
      num = num + length(j);
  
      % As a percentage
      conf = 100*(length(test(:,1)) - num) / length(test(:,1));
  
      title(['Height ',num2str(conf,'%.3f'),'% measurements within 99% Confidence']);
    end
    grid;
    ylabel('Height Error [m]');
  
    datalength = test(end,1) - test(1,1);

    % Now go through each plot and set the y-axis scale to the optimal
    % setting. Only do this for fixed RTK.
    if(PosType == 76)
      % East Graph
      subplot(3,1,1)
      maxE =  ceil(maxE*40)/40;  % closest 2.5cm
      minE = floor(minE*40)/40;  % closest 2.5cm
      A = axis;
      axis([ A(1:2) minE maxE]);
      
      % North Graph
      subplot(3,1,2)
      maxN =  ceil(maxN*40)/40;  % closest 2.5cm
      minN = floor(minN*40)/40;  % closest 2.5cm
      A = axis;
      axis([ A(1:2) minN maxN]);

      % Height Graph
      subplot(3,1,3)
      maxU =  ceil(maxU*40)/40;  % closest 2.5cm
      minU = floor(minU*40)/40;  % closest 2.5cm
      A = axis;
      axis([ A(1:2) minU maxU]);
    end
  
    if(AxisType == 0)
      labelAxes([],[],[],[],1);
      xlabel(['Time [GPS Seconds] - data length = ',num2str(datalength,'%.0f'),' Seconds']);
    else
      % Want a HH:MM axis, first figure out how long the data is in hours
  
      hours = datalength / 3600;
  
      if(hours <= 1)
        % 10 minute ticks
        minTime = floor(test(1,1) / 600) * 600;
        maxTime = ceil(test(end,1) / 600) * 600;
        Ticks = (maxTime - minTime)/600;
      elseif(hours <= 2)
        % 15 minute ticks
        minTime = floor(test(1,1) / 900) * 900;
        maxTime = ceil(test(end,1) / 900) * 900;
        Ticks = (maxTime - minTime)/900;
      elseif(hours <= 4)
        % 30 minute ticks
        minTime = floor(test(1,1) / 1800) * 1800;
        maxTime = ceil(test(end,1) / 1800) * 1800;
        Ticks = (maxTime - minTime)/1800;
      elseif(hours <= 8)
        % 1 hour ticks
        minTime = floor(test(1,1) / 3600) * 3600;
        maxTime = ceil(test(end,1) / 3600) * 3600;
        Ticks = (maxTime - minTime)/3600;
      elseif(hours <= 16)
        % 2 hour ticks
        minTime = floor(test(1,1) / 7200) * 7200;
        maxTime = ceil(test(end,1) / 7200) * 7200;
        Ticks = (maxTime - minTime)/7200;
      else
        minTime = floor(test(1,1) / 3600) * 3600;
        maxTime = ceil(test(end,1) / 3600) * 3600;
        hours = (maxTime - minTime)/3600;
  
        if(hours == 17)
          maxTime = maxTime + 3600;
          Ticks = 6;  % every 3 hours
        elseif(hours == 18)
          Ticks = 6;  % every 3 hours
        elseif(hours == 19)
          maxTime = maxTime + 3600;
          Ticks = 5;  % every 4 hours
        elseif(hours == 20)
          Ticks = 5;  % every 4 hours
        elseif(hours == 21)
          Ticks = 7;  % every 3 hours
        elseif(hours == 22)
          minTime = minTime - 3600;
          maxTime = maxTime + 3600;
          Ticks = 6;  % every 4 hours
        elseif(hours == 23)
          maxTime = maxTime + 3600;
          Ticks = 6;  % every 4 hours
        elseif(hours == 24)
          Ticks = 6;  % every 4 hours
        elseif(hours == 25)
          Ticks = 5;  % every 5 hours
        elseif(hours == 26)
          Ticks = 8;  % every 3.25 hours
        elseif(hours == 27)
          maxTime = maxTime + 3600;
          Ticks = 7;  % every 4 hours
        elseif(hours == 28)
          Ticks = 7;  % every 4 hours
        elseif(hours == 29)
          maxTime = maxTime + 3600;
          Ticks = 6;  % every 5 hours
        elseif(hours == 30)
          Ticks = 6;  % every 4 hours
        elseif(hours == 31)
          maxTime = maxTime + 3600;
          Ticks = 8;  % every 4 hours
        elseif(hours == 32)
          Ticks = 8;  % every 4 hours
        elseif(hours == 33)
          maxTime = maxTime + 3600;
          minTime = minTime - 3600;
          Ticks = 7;  % every 5 hours
        elseif(hours == 34)
          maxTime = maxTime + 3600;
          Ticks = 7;  % every 5 hours
        elseif(hours == 35)
          Ticks = 7;  % every 5 hours
        else
          % 8 hour ticks
          minTime = floor(test(1,1) / 28800) * 28800;
          maxTime = ceil(test(end,1) / 28800) * 28800;
          Ticks = (maxTime - minTime)/28800;
        end
      end
  
      subplot(3,1,1);
      A = axis;
      axis([minTime maxTime A(3:4)]);
      subplot(3,1,2);
      A = axis;
      axis([minTime maxTime A(3:4)]);
      subplot(3,1,3);
      A = axis;
      axis([minTime maxTime A(3:4)]);
      if(isnan(LocalOffset))
        labelAxes(Ticks,[],[],[],0);
        xlabel(['Time UTC [HH:MM] - data length = ',num2str(datalength,'%.0f'),' Seconds']);
      else
        labelAxes(Ticks,[],[],[],0,LocalOffset);
        if(LocalOffset >= 0)
          xlabel(['Time Local (UTC + ',num2str(LocalOffset,'%.0f'),' Hours) [HH:MM] - data length = ',num2str(datalength,'%.0f'),' Seconds']);
        else
          xlabel(['Time Local (UTC - ',num2str(abs(LocalOffset),'%.0f'),' Hours) [HH:MM] - data length = ',num2str(datalength,'%.0f'),' Seconds']);
        end
      end
    end
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_ENU.',FileExtension]);
    end
    
    %
    % Plot the number of satellites and PDOP
    % 
    
    figure
    if(FileType == 0)
      subplot(2,1,1)
    else
      subplot(3,1,1)
    end
    plot(out(:,1),out(:,6));
    grid
    ylabel('PDOP');
   
    if(FileType == 0)
      subplot(2,1,2)
    else
      subplot(3,1,2)
    end
    plot(out(:,1),out(:,5));
    grid
    ylabel('# Satellites in Solution');
    
    if(FileType == 1)
      % Plot the RMS
      subplot(3,1,3);
      plot(out(:,1),out(:,8));
      grid
      ylabel('Reported RMS [m]');
      subplot(3,1,1);
    else
      subplot(2,1,1);
    end
      
    title([PREFIX,': PDOP/Number of satellites in solution'],'Interpreter','none');
    if(AxisType == 0)
      if(FileType == 0)
        subplot(2,1,2)
      else
        subplot(3,1,3)
      end
      labelAxes([],[],[],[],1);
      xlabel('Time [GPS Seconds]');
    else
      if(FileType == 0)
        subplot(2,1,1)
        A = axis;
        axis([minTime maxTime A(3:4)]);
        subplot(2,1,2)
        A = axis;
        axis([minTime maxTime A(3:4)]);
      else
        subplot(3,1,1)
        A = axis;
        axis([minTime maxTime A(3:4)]);
        subplot(3,1,2)
        A = axis;
        axis([minTime maxTime A(3:4)]);
        subplot(3,1,3)
        A = axis;
        axis([minTime maxTime A(3:4)]);
      end
  
      if(isnan(LocalOffset))
        labelAxes(Ticks,[],[],[],0);
        xlabel('Time UTC [HH:MM]');
      else
        labelAxes(Ticks,[],[],[],0,LocalOffset);
        if(LocalOffset >= 0)
          xlabel(['Time Local (UTC + ',num2str(LocalOffset,'%.0f'),' Hours) [HH:MM]']);
        else
          xlabel(['Time Local (UTC - ',num2str(abs(LocalOffset),'%.0f'),' Hours) [HH:MM]']);
        end
      end
    end
  
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_DOP.',FileExtension]);
    end
    
    %
    % Plot the east vs north error
    % 
  
    figure
    plot(testStats(:,2),testStats(:,3),'x');
    xlabel('East Error [m]');
    ylabel('North Error [m]');
    grid
    title([PREFIX,': Horizontal Error'],'Interpreter','none');
    
    %
    % Create a 1-sigma ellipse
    %
    rectangle('Position', [  mean(testStats(:,2)) - std(testStats(:,2)), ...
                             mean(testStats(:,3)) - std(testStats(:,3)), ...
                             2*std(testStats(:,2)), ...
                             2*std(testStats(:,3))], ...
              'Curvature',[1 1]);
    
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_HORZ.',FileExtension]);
    end
    
    %
    % Process the 3-D Data
    %
    
    figure
    subplot(2,1,1);
    ThreeDErr = sqrt(testStats(:,2).^2 + testStats(:,3).^2 + testStats(:,4).^2); 
    N=hist(ThreeDErr,range);
    bar(range,100*N./len,1);
    xlabel('3-D Position Error [m]');
    ylabel('Number of Epochs [%]');
    A = axis;
    
    % 
    % Find the first and last zero so we can figure out the best axis.
    %
    S = length(N);
    START = TRUE;
    
    for i=1:S
      if(N(i) ~= 0)
        if(START == TRUE)
          FirstNonZero = i;
          START = FALSE;
        end
        LastZero = i;
      end    
    end
    
    if(PosType == 76)
      % Closest mm for fixed RTK
      MinX = floor(range(FirstNonZero)*1000)/1000;
      MaxX = ceil(range(LastZero)*1000)/1000;
    else
      % Closest m for everything else
      MinX = floor(range(FirstNonZero));
      MaxX = ceil(range(LastZero));
    end
    A(1) = MinX;
    A(2) = MaxX;
    axis(A);
    grid
    
    
    subplot(2,1,2)
    X = cumsum(N);
    plot(range,X.*100/length(testStats));
    grid
    axis([A(1:2) 0 100]);
    xlabel('3D Position Error [m]');
    ylabel('Cumulative Proability [%]');
    
    %
    % Find the closest point to 50% == CEP
    % 
    i = find (X/length(testStats) < .5);
    i_max = find( i == max(i));
    
    %
    % Now interpolate to provide a better estimate of
    % the 95 percentile.
    % 
    
    if(i_max ~= size(range))
      % We are not on the last element therefore we can interpolate
      Y_3D(1) = X(i_max)/length(testStats);
      Y_3D(2) = X(i_max+1)/length(testStats);
      X_3D(1) = range(i_max);
      X_3D(2) = range(i_max+1);
      m         = (Y_3D(2) - Y_3D(1)) / (X_3D(2) - X_3D(1));
      c         = Y_3D(1) - m*X_3D(1);
      Spherical_50= (0.5 - c)/m;
    else
      Spherical_50 = range(i_max);
    end
    
    %
    % Find the closest point to 95%
    % 
    i = find (X/length(testStats) < .95);
    i_max = find( i == max(i));
    
    %
    % Now interpolate to provide a better estimate of
    % the 95 percentile.
    % 
    
    if(i_max ~= size(range))
      % We are not on the last element therefore we can interpolate
      Y_3D(1) = X(i_max)/length(testStats);
      Y_3D(2) = X(i_max+1)/length(testStats);
      X_3D(1) = range(i_max);
      X_3D(2) = range(i_max+1);
      m         = (Y_3D(2) - Y_3D(1)) / (X_3D(2) - X_3D(1));
      c         = Y_3D(1) - m*X_3D(1);
      Spherical_95= (0.95 - c)/m;
    else
      Spherical_95 = range(i_max);
    end
    
    Spherical_RMS = sqrt(sum(ThreeDErr.^2) / length(ThreeDErr)); 
    
    subplot(2,1,1);
    title([PREFIX,': 3-D Error: 95%=',num2str(Spherical_95,'%.3f'),'m'],'Interpreter','none');
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_3D_Dist.',FileExtension]);
    end
    
    fprintf(ftxt,'3D,%.3f,%.3f,%.3f,%.3f\n', ...
                   Spherical_RMS, ...
                   Spherical_50, ...
                   Spherical_95, ...
                   max(ThreeDErr));
                   
    if( (PosType == 76) & (Filter == TRUE) )
     % Provide some info on the position types
     fprintf(ftxt,'Data Types, Num Expected, Num Epochs, Num 3D-Fixed (Good), Num 3D-Fixed (Bad), Num 3D-Float, Num 3D-DGPS, Num 3D-Autonomous\n');
     fprintf(ftxt,' ,%d,%d,%d,%d,%d,%d,%d\n', ...
                expectedLength, ...
                dataLength, ...
                NumGood, ...
                Type76 - NumGood, ...
                Type68, ...
                Type20, ...
                Type4);
    end            
    
    %
    % Now produce a 3-D plot
    %
    figure
    plot3(testStats(:,2), testStats(:,3), testStats(:,4),'rx');
    xlabel('East Error [m]');
    ylabel('North Error [m]');
    zlabel('Height Error [m]');
    title([PREFIX,': 3D Error'],'Interpreter','none');
    grid
    if(savePlots == TRUE)
      print(OutputGraphicsType,[PREFIX,'_3D.',FileExtension]);
    end
    
    TESTCode = FALSE;
    if(TESTCode == TRUE) % Never Get Here
      %
      % Now plot in the time domain
      %
    
      StartTime = test(1,1);
      figure
      subplot(3,1,1)
      hold;
      X= (test(:,1)-StartTime)./3600;
      %
      % Plot the red sigma line first so the legend command picks it up.
      % However, replot it later as the data may cover it!
      %
      plot(X, sqrt(test(:,8)),'--r');
      plot(X,test(:,2));
      plot(X, sqrt(test(:,8)),'--r');
      plot(X,-sqrt(test(:,8)),'--r');
      legend('+/- 1-Sigma Error Estimates',0);
      title([PREFIX,': East/North/Height Error'],'Interpreter','none');
      grid
      ylabel('East Error [m]');
    
      subplot(3,1,2)
      plot(X,test(:,3));
      hold;
      plot(X, sqrt(test(:,9)),'--r');
      plot(X,-sqrt(test(:,9)),'--r');
      grid
      ylabel('North Error [m]');
    
      subplot(3,1,3)
      plot(X,test(:,4));
      hold;
      plot(X, sqrt(test(:,10)),'--r');
      plot(X,-sqrt(test(:,10)),'--r');
      grid
      ylabel('Height Error [m]');
      xlabel('Time [GPS Seconds]');
      if(savePlots == TRUE)
        print(OutputGraphicsType,[PREFIX,'_ENU_ERR.',FileExtension]);
      end
    end
    
    if( Compiled == TRUE)
      % Graphs displayed during processing, but must close now to enable the
      % compiled version to terminate back to DOS.
      % Note: Automated use from Matlab appends an 'exit' to the end of the
      % m-file ('exit' cannot be compiled, MEX files only).
      close all
    end
  end  
  fclose(ftxt);
end  
