function load_pi(Compiled,  ...
                 CodeLimit, ...
                 CarrLimit, ...
                 PlotPNGIn,   ...
                 doFFTPlotIn, ...
                 doACFPlotIn, ...
                 QuickModeIn, ...
                 CNRTestElMaskIn, ...
                 CNRTestSigmaMaskIn, ...
                 SkipCNRTestsIn, ...
                 dropLegendIn, ...
                 decimateIn, ...
                 startTimeIn, ...
                 stopTimeIn, ...
                 matFileConfigIn ...
                 ); % , ...
               % saveFileLabel);

               
%
% $Id: load_pi.m,v 1.64 2008/02/25 16:35:30 stuart Exp $
% 
% load_pi(Compiled,CodeLimit,CarrLimit,PlogPNG,doFFTPlot,doACFPlot);
%
% Where:
%   Compiled  - (=0) Called by Matlab; (=1) StandAlone application.
%   CodeLimit - Maximum Y-axis for pseudorange
%   CarrLimit - Maximum Y-axis for carrier phase
%   PlotPNGIn   - (=0) Don't Create .png files; (=1) Create .png files
%   doFFTPlotIn - (=0) No FFT Plots; (=1) Plots the FFT of the Double
%               difference residual, useful for pikcing out the 
%               dominant multipath frequencies.
%   doACFPlotIn - (=0) No ACF Plots;  (=1) Finds the longest contiguous
%               segment of DD data for each satellite and computes the
%               ACF, if and only if the data segment is longer than
%               1000 seconds. Note the ACF function is written to a text
%               file for future analysis.
%   QuickModeIn - Performs a sub-set of the processing.
%               0 = Full processing
%               1 = Sub-set of observable and C/No
%               2 = Full C/No processing only
%               3 = C/No Base-Rover processing only
%   CNRTestElMaskIn - used by the C/No data analysis code. Added to
%                   support Pogo testing. Data included in the
%                   Base-Rover versus time plots for L1 and L2, and the
%                   tabulated summary statistics must be above this
%                   elevation mask if set. If this parameter is missing
%                   all data present is included. Units are degrees.
%   CNRTestSigmaMaskIn - used by the C/No data analysis code. Added to 
%                      support Pogo testing. The standard deviation of 
%                      the Base-Rover difference must be lower than this
%                      sigma on a per satellite basis in order for the
%                      satellite data to be rolled up into the overall
%                      average displayed at the bottom of the C/No
%                      difference HTML table.
%   SkipCNRTestsIn - (=0) Process all tests, (=1) do not process any of
%                       the CNR tests. Useful for .dat file processing
%                       where the original data was RINEX with no C/No
%                       field.
%   dropLegendIn  -  (=0) Output legends as normal, (=1) do not plot the
%                       legends. In some of our testing when there is
%                       lots of data we have the problem that the legend
%                       over-writes some of the data. To prevent this
%                       turn off the legend.
%   decimateIn - (default = 0 - No decimation). Provides the interval 
%                in seconds that the input data should be decimated to
%                when it is read in. 
%   startTime  - (default = [] - No start time). GPS seconds of week.
%   stopTime   - (default = [] - No stop time).  GPS seconds of week.
% 
%   matFileConfigIn - Configure the use of MAT-files to capture output of 
%                   the pilineMex.dll for reuse. 
%        0 = don't use mat files
%        1 = save output of pilineMex.dll in piProData.mat
%        2 = use piProData.mat instead of running pilineMex.dll
%        
% 
% 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:
%   load_pi( '1' )  Sets other parameters to defaults.
% 
%   The Defaults are:
%
%       Compiled           = 0;
%       CodeLimit          = 4;     % [Metres]
%       CarrLimit          = 100;   % [Millicycles]
%       PlotPNGIn          = 1;
%       doFFTPlotIn        = 0;
%       doACFPlotIn        = 0;
%       QuickModeIn        = 0;
%       CNRTestElMaskIn    = 0;     % [Degrees]
%       CNRTestSigmaMaskIn = 10;    % [dB]
%       SkipCNRTestsIn     = 0;
%       dropLegendIn       = 0;
%       decimateIn         = 0;
%
% Description:
%  Calls pilineMex.dll and creates the double difference observables
%  which are plotted in a variety of ways. Using the command line
%  options the axis limits on the double difference graphs may be set
%  and whether a PNG file is saved may be controlled. However, unless
%  you are debugging this file note various figures are deleted as soon
%  as they are created so PNG files should always be saved. 
%
% Output:
%  In addition to the .PNG files statistics for each satellite are saved
%  to two files, a tracking plot is created and polar plots of the C/No
%  against elevation/azimuth.
%
%  stats.txt:
%    A space delimited file of the format
%    COL1  - PRN #
%    COL2  - Number L1 epochs
%    COL3  - Number L2 epochs
%    COL4  - Pseudorange L1 mean  [cm]
%    COL5  - Pseudorange L1 sigma [cm]
%    COL6  - Pseudorange L2 mean  [cm]
%    COL7  - Pseudorange L2 sigma [cm]
%    COL8  - Phase L1 mean  [mcycle]
%    COL9  - Phase L1 sigma [mcycle]
%    COL10 - Phase L2 mean  [mcycle]
%    COL11 - Phase L2 sigma [mcycle]
%
%  stats.htm:
%    HTML table containing the above data and HTML links to the .PNG
%    files
%
%    NOTE: Full support for .PNG is not available in Microsoft IE before
%    version 5.5. If you do not see the full graph and there are no
%    scroll bars you need to update your browser. The level of Netscape 
%    support is unkown.
%    
% Prerequisite:
%   There must be a valid piline.ini file (same format as the old
%   survey.dat file). If this file is edited it is possible to still use
%   the WIN32 console piline application. In which case Matlab would
%   expect the following files to be present in the working directory.
%       rg_dd_l1.out
%       rg_dd_l2.out
%       ph_dd_l1.out
%       ph_dd_l2.out
%
%   We must have used the -o64 Piline option; only need to know about
%   this if the command line piline is used.
%   This provides us with data as follows:
%   
%   ph_dd_l1 -> [DD EL AZ C/No]
%   rg_dd_l1 -> [DD EL C/No]
%   ph_dd_l2 -> [DD EL C/No]
%   rg_dd_l2 -> [DD EL C/No]
%   
% Other Matlab files needed:
%
%  do_xcorr.m
%  doLogo.m
%  excise_matrix.m
%  excise_vector.m
%  fftplot.m
%  labelaxes.m
%  loadPiHeaderHTM.m
%  PlotCNR.m
%  PlotCNRDiffElev.m
%  PlotCNRDiffTime.m
%  plotHist.m
%  plotIono.m
%  plotSVTrack.m
%  plotTrend.m
%  polar_expanded.m
%  runPiline.m
%  setColourCodes.m
%  suptitle.m
%  xcorr.m
%
%  pilineMEX.dll
% 
% Stuart Riley 6/8/2001
% Copyright Trimble Navigation Ltd. 2000 - 2001
%
% Updated to enable Matlab compiling as a StandAlone application
% Sept, 2000, Timo Allison. When used as StandAlone, the passed
% parameter 'Compiled' must be set to 1.


% Compiled  = 0;
% CodeLimit = 2;
% CarrLimit = 200;
% PlotPNGIn   = 1;
% nargin = 4;

% doFFTPlotIn = 0;
% doACFPlotIn = 0;
% QuickModeIn = 0;
% CNRTestElMaskIn    = 0;
% CNRTestSigmaMaskIn = 10;
% SkipSNRTestsIn = 0;
% nargin = 9;

%
% Print the version
% 
piProVer = pipro_version;

if strcmp(piProVer.buildType,'Test')
  fprintf('Running a %s build of PiPro derived from %s, which uses %s.\n', ...
    piProVer.buildType, piProVer.piproVersion, piProVer.pilineVersion);
else
  fprintf('Running a %s build of PiPro with CVS Tag %s, which uses %s.\n', ...
    piProVer.buildType, piProVer.piproVersion, piProVer.pilineVersion);
end

%
%********** LIST OF ALL DEFINES **********
%

L1 = 1;
L1_CA = 1; % GLONASS
L2 = 2;
L5 = 3;
L1_P = 3;  % GLONASS
LIGHT = 2.99792458e8;
L1_WAVELENGTH =  LIGHT / 1575.42e6;
L2_WAVELENGTH =  LIGHT / 1227.6e6;
L5_WAVELENGTH =  LIGHT / 1176.45e6;
colour = 0;
TRUE  = 1;
FALSE = 0;
ASCII_ZERO = 48 ;
ASCII_ONE  = 49 ;

% Enumeration for matfile use
MATFILE_NONE = 0;
MATFILE_SAVE = 1;
MATFILE_LOAD = 2; 

% 
% Load the RGB values we will use on the plots.
%
colour_LUT = setColourCodes;
disp('Defines set...');

% 
%********** END OF DEFINE LIST ***********
%

% Prevent Matlab from running out of resources.
% This fix was suggested by Will Lentz, Feb 6, 2001.
set(0,'defaultFigureToolbar','none') ;

% Set the defaults:
CodeMax = 4;   % Metres
CarrMax = 100; % Millicycles
PlotPNG = TRUE;
doFFTPlot = 0;
doACFPlot = 0;
QuickMode = 0;
CNRTestElMask    = 0;
CNRTestSigmaMask = 10;
SkipCNRTests = 1;
dropLegend = 0;
decimate = 0;
startTime = [] ;
stopTime = [] ;
matFileConfig = MATFILE_NONE;

% Threshold value to auto-scale the DD plots rather than
% applying the user defined values
CodeThres_AScale = CodeMax*2;
CarrThres_AScale = CarrMax*2;

% 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 > 0 )
  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 ;
  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
  Compiled = FALSE ;
end

if nargin >= 2
  if Compiled == TRUE
    arg = str2num( CodeLimit ) ;
  else
    arg = CodeLimit ;
  end
  if ~isempty( arg )
    CodeMax = arg ;
  end
end

if nargin >= 3
  if Compiled == TRUE
    arg = str2num( CarrLimit ) ;
  else
    arg = CarrLimit ;
  end
  if ~isempty( arg )
    CarrMax = arg ;
  end
end

if nargin >= 4
  if Compiled == TRUE
    arg = str2num( PlotPNGIn ) ;
  else
    arg = PlotPNGIn ;
  end
  if ~isempty( arg )
    PlotPNG = arg ;
  end
end

if nargin >= 5
  if Compiled == TRUE
    arg = str2num( doFFTPlotIn );
  else
    arg = doFFTPlotIn ;
  end
  if ~isempty( arg )
    doFFTPlot = arg ;
  end
end

if nargin >= 6
  if Compiled == TRUE
    arg = str2num( doACFPlotIn );
  else
    arg = doACFPlotIn ;
  end
  if ~isempty( arg )
    doACFPlot = arg ;
  end
end

if nargin >= 7
  if Compiled == TRUE
    arg = str2num( QuickModeIn );
  else
    arg = QuickModeIn ;
  end
  if ~isempty( arg )
    QuickMode = arg ;
  end
end

if nargin >= 8
  if Compiled == TRUE
    arg = str2num( CNRTestElMaskIn );
  else
    arg = CNRTestElMaskIn ;
  end
  if ~isempty( arg )
    CNRTestElMask = arg ;
  end
end

if nargin >= 9  
  if Compiled == TRUE
    arg = str2num( CNRTestSigmaMaskIn );
  else
    arg = CNRTestSigmaMaskIn ;
  end
  if ~isempty( arg )
    CNRTestSigmaMask = arg ;
  end
end

if nargin >= 10
  if Compiled == TRUE
    arg = str2num( SkipCNRTestsIn );
  else
    arg = SkipCNRTestsIn ;
  end
  if ~isempty( arg )
    SkipCNRTests = arg ;
  end
end

if nargin >= 11
  if Compiled == TRUE
    arg = str2num( dropLegendIn );
  else
    arg = dropLegendIn ;
  end
  if ~isempty( arg )
    dropLegend = arg ;
  end
end

if nargin >= 12
  if Compiled == TRUE
    arg = str2num( decimateIn );
  else
    arg = decimateIn ;
  end
  if ~isempty( arg )
    decimate = arg ;
  end
end

if nargin >= 13 
  if Compiled == TRUE 
    arg = str2num( startTimeIn ) ;
  else
    arg = startTimeIn;
  end
  if ~isempty( arg )
    startTime = arg;
  end
end

if nargin >= 14 
  if Compiled == TRUE 
    arg = str2num( stopTimeIn ) ;
  else
    arg = stopTimeIn;
  end
  if ~isempty( arg )
    stopTime = arg;
  end
end

if nargin >= 15
  if Compiled == TRUE
    arg = str2num( matFileConfigIn ) ;
  else
    arg = matFileConfigIn;
  end
  if ~isempty( arg )
    matFileConfig = arg ;
  end
end

disp('Command line arguments accepted...');
%*********************************************************************
%                  End of input argument Processing
%*********************************************************************

% Use pilineMEX.dll

%%%%%%%%%%%%%%%%%%%%%% Global data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global phaseDDiffL1  % Transferred from PiLine mex file.
global phaseDDiffL2  % Transferred from PiLine mex file.
global phaseDDiffL2C % Transferred from PiLine mex file.
global phaseDDiffL5  % Transferred from PiLine mex file.
global ionoDDiff     % Transferred from PiLine mex file.
global rangeDDiffL1  % Transferred from PiLine mex file.
global rangeDDiffL2  % Transferred from PiLine mex file.
global rangeDDiffL2C % Transferred from PiLine mex file.
global rangeDDiffL5  % Transferred from PiLine mex file.
global statRangeL1   % Transferred from PiLine mex file.
global statRangeL2   % Transferred from PiLine mex file.
global rangeDDiffLC  % Transferred from PiLine mex file.
global statRangeL5   % Transferred from PiLine mex file.
global rec17DataL1   % Transferred from PiLine mex file.
global rec17DataL2   % Transferred from PiLine mex file.
global rangeDDiffLC  % Transferred from PiLine mex file.
global rec17DataL5   % Transferred from PiLine mex file.

global phaseDDiffGL1_CA 
global phaseDDiffGL1_P
global phaseDDiffGL2
global rangeDDiffGL1_CA
global rangeDDiffGL1_P
global rangeDDiffGL2
  
% piLineParams structure fields:
%   refSV: 1x1 mxREAL array (C double)
%   stn1FileName: string (C string)
%   stn2FileName: string (C string)
global piLineParams  % Transferred from PiLine mex file.

if ( matFileConfig == MATFILE_LOAD )
  load piProData.mat;
else
  runPiLine(0,1,0,[],[],decimate,startTime,stopTime);
  if ( matFileConfig == MATFILE_SAVE )
    save piProData.mat;
  end
end

GLONASSInstalled = (length(rangeDDiffGL1_CA) > 0);
GL1_CA_present = GLONASSInstalled;
GL1_P_present = (length(rangeDDiffGL1_P) > 0);
GL2_P_present = (length(rangeDDiffGL2) > 0);
GL2_CA_present = 0; % Currently GL2_CA is not supported in PiPro (load_pi)

if ( GLONASSInstalled )

  % Align GLN L1 CA with GLN L1 P
  if ( GL1_P_present )
    fprintf('Realigning GLN L1-C/A and GLN L1-P arrays\n');
    [C,IA,IB] = intersect(phaseDDiffGL1_CA(:,1), phaseDDiffGL1_P(:,1));
    if ( ~isempty(C) )
      phaseDDiffGL1_CA = phaseDDiffGL1_CA(IA,:);
      phaseDDiffGL1_P  = phaseDDiffGL1_P(IB,:);
      rangeDDiffGL1_CA = rangeDDiffGL1_CA(IA,:);
      rangeDDiffGL1_P  = rangeDDiffGL1_P(IB,:);
    end
  end

  % Align GLN L1 CA with GLN L2(P)
  % Also, align GLN L1 P with the computed alignment. 
  if ( GL2_P_present )

      fprintf('Realigning GLN L1-C/A and GLN L2-P and GLN L1-P arrays\n');
      [C,IA,IB] = intersect(phaseDDiffGL1_CA(:,1), phaseDDiffGL2(:,1));
      if ( ~isempty(C) )
        
        rangeDDiffGL1_CA = rangeDDiffGL1_CA(IA,:);
        phaseDDiffGL1_CA = phaseDDiffGL1_CA(IA,:);
        if ( GL1_P_present )
          rangeDDiffGL1_P  = rangeDDiffGL1_P(IA,:);
          phaseDDiffGL1_P  = phaseDDiffGL1_P(IA,:);
        end

        rangeDDiffGL2    = rangeDDiffGL2(IB,:);
        phaseDDiffGL2    = phaseDDiffGL2(IB,:);
      end
  end
    
end

NumGLNFigs = GL1_CA_present + GL1_P_present + GL2_P_present;


%
% GLONASS Almanac to FDMA look-up-table (we'll eventually get this from
% runPiLine)
%
% The frequency channel mappings are available from:
% <http://www.glonass-ianc.rsa.ru/pls/htmldb/f?p=202:20:2674595394375172305::NO>
%
% The "slot" described on this web page is the same as what we call svId. 
% 
% The channel mappings below are written in the format:
%     FDMAChan(svId) = FrequencyChannel
% 
FDMAChan(1)  =    7;
FDMAChan(2)  =  NaN;
FDMAChan(3)  =  NaN;
FDMAChan(4)  =    6;
FDMAChan(5)  =  NaN;
FDMAChan(6)  =    1;
FDMAChan(7)  =    5;
FDMAChan(8)  =    6;
FDMAChan(9)  =   -2;
FDMAChan(10) =    4;
FDMAChan(11) =    0;
FDMAChan(12) =  NaN;
FDMAChan(13) =   -2;
FDMAChan(14) =    4;
FDMAChan(15) =    0;
FDMAChan(16) =  NaN;
FDMAChan(17) =   -1;
FDMAChan(18) =  NaN;
FDMAChan(19) =    3;
FDMAChan(20) =    2;
FDMAChan(21) =  NaN;
FDMAChan(22) =  NaN;
FDMAChan(23) =    3;
FDMAChan(24) =    2;


% ionoDDiff only needed for type 0 processing.
if( QuickMode > 0) 
  clear ionoDDiff;
end

Range_Sigma = ones(3,32).*NaN;
Range_Mean  = ones(3,32).*NaN;
Phase_Sigma = ones(3,32).*NaN;
Phase_Mean  = ones(3,32).*NaN;
NumEpochs   = ones(3,32).*NaN;

if(GLONASSInstalled)
  GLN_Range_Sigma = ones(3,24).*NaN;
  GLN_Range_Mean  = ones(3,24).*NaN;
  GLN_Phase_Sigma = ones(3,24).*NaN;
  GLN_Phase_Mean  = ones(3,24).*NaN;
  GLN_NumEpochs   = ones(3,24).*NaN;
end


NumFigs = 1; % Always have L1!

if(length(phaseDDiffL2) == 0)
  L2Present = FALSE;
else
  L2Present = TRUE;
  NumFigs = NumFigs + 1;
end

if(length(phaseDDiffL5) == 0)
  L5Present = FALSE;
else
  L5Present = TRUE;
  NumFigs = NumFigs + 1;
end

disp('All data read in and initialization complete...');

% %%%%% Save the Data and Break
% 
% save([saveFileLabel 'phaseDDiffL1.txt'],'-ASCII','-DOUBLE','phaseDDiffL1');
% save([saveFileLabel 'phaseDDiffL2.txt'],'-ASCII','-DOUBLE','phaseDDiffL2');
% save([saveFileLabel 'rangeDDiffL1.txt'],'-ASCII','-DOUBLE','rangeDDiffL1');
% save([saveFileLabel 'rangeDDiffL2.txt'],'-ASCII','-DOUBLE','rangeDDiffL2');
% 
% disp('Saved data. Exiting ...');
% return; 
% 
% %%%%% End of Saving Data


%
% ------------------------------------------------------------------
% End of initialization code. We now have all the data read into the
% Matlab environment and various variables initialized.
% ------------------------------------------------------------------
% 

% Get the start/stop times for axis scaling.
start = rangeDDiffL1(1,1);
stop  = rangeDDiffL1(size(rangeDDiffL1,1),1);

num = 2;
max_l1_phase = 0;
min_l1_phase = 0;
max_l2_phase = 0;
min_l2_phase = 0;
max_l5_phase = 0;
min_l5_phase = 0;

PHASE = 1;
RANGE = 2;

%
% Set up a string that we will use on the plots, provides the total
% amount time of the data set and the base/rover file names.
%

totalTime  = stop - start ;
totalHours = floor( totalTime / 3600 ) ;
totalMins  = round( mod( totalTime, 3600 ) / 60 ) ;
if(totalMins == 60)  % Can occur due to rounding.
  totalMins = 0 ;
  totalHours = totalHours + 1 ;
end

timeStr = [ num2str(totalHours,'%02d' ),'hrs:', ...
            num2str(totalMins,'%02d'),'mins - files Base: ', ...
            piLineParams.stn1FileName,' Rover: ', ...
            piLineParams.stn2FileName ] ;

% To make everything self contained the Trimble logo .png file is 
% created by a .m function.
doLogo;

if(QuickMode < 2)
  % Open the result HTML File
  fid  = fopen('stats.htm','w');
  % Start the HTML document
  fprintf(fid,'<HTML><HEAD><title>PiLine Analysis Results</title></head>\n');
  fprintf(fid,'<BODY background="../backgrnd.gif">');
  fprintf(fid,'<img src="logo.png">\n');
  fprintf(fid,'<font face="arial">\n');
end

if(QuickMode == 0)
  FigAll(PHASE) = figure('Visible','off');
  FigAll(RANGE) = figure('Visible','off');
  if (GLONASSInstalled)
    FigAllGln(PHASE) = figure('Visible','off');
    FigAllGln(RANGE) = figure('Visible','off');
  end

  AScale_GPSPhase_L1 = zeros(32,1);
  AScale_GPSPhase_L2 = zeros(32,1);
  AScale_GPSPhase_L5 = zeros(32,1);
  %
  % Loop around for each potential satellite (PRN #) and extract the L1 &
  % L2 double difference phase information. Plot this and compute the
  % statistics.
  %
  for sv=1:32
    l1_phase = find(isfinite(phaseDDiffL1(:,(sv+1))));
    if(L2Present == TRUE)
      l2_phase = find(isfinite(phaseDDiffL2(:,(sv+1))));
    else
      l2_phase = [];
    end
    
    if(L5Present == TRUE)
      l5_phase = find(isfinite(phaseDDiffL5(:,(sv+1))));
    else
      l5_phase = [];
    end
  
    if((size(l1_phase,1) > 0) | (size(l2_phase,1) > 0) | (size(l5_phase,1) > 0))
      num = num + 1;
      f(num) = figure('Visible','off');
      Sat(num) = sv;
      colour = colour + 1;
  
      % Make sure there is some L1 data to plot.
      if(size(l1_phase,1) > 0)
        l1_phase_mean  = mean(phaseDDiffL1(l1_phase,(sv+1)));
        l1_phase_sigma = std(phaseDDiffL1(l1_phase,(sv+1)));
        subplot(NumFigs,1,1);
        plot(phaseDDiffL1(:,1), phaseDDiffL1(:,(sv+1))./10);
        grid;
        title(['SV = ',num2str(sv),' L1 Carrier Phase DD \mu = ', ...
        num2str(l1_phase_mean / 10,'%.1f'),'mcycle; \sigma = ', ...
        num2str(l1_phase_sigma / 10,'%.1f'),'mcycle']);
          ylabel('L1 DD [mcycle]');
        if (max(abs(phaseDDiffL1(:,(sv+1))./10)) > CarrThres_AScale)
          AScale_GPSPhase_L1(sv) = 1;
        end
  
        % If there is no L2 data put the label under the L1 axis, otherwise
        % it gets placed at the bottom of the plot under the L2.
        if(NumFigs == 1)
          xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
        end
  
        % Contribute to the figure that plots all data superimposed
        % figure(FigAll(PHASE));
        set(0, 'CurrentFigure', FigAll(PHASE));

        subplot(NumFigs,1,1);

        plot(phaseDDiffL1(:,1), phaseDDiffL1(:,(sv+1))./10,'color', ...
           colour_LUT(colour,:));
        hold on;
        % figure(f(num));
        set(0, 'CurrentFigure', f(num));
      
        Phase_Sigma(L1,sv) = l1_phase_sigma/10;
        Phase_Mean (L1,sv) = l1_phase_mean/10;
      
        tmp = max(phaseDDiffL1(l1_phase,(sv+1)));
        if(tmp > max_l1_phase)
          max_l1_phase = tmp;
        end
        tmp = min(phaseDDiffL1(l1_phase,(sv+1)));
        if(tmp < min_l1_phase)
          min_l1_phase = tmp;
        end
      else
        l1_phase_mean  = NaN;
        l1_phase_sigma = NaN;
      end  % End of L1 processing loop
  
      % Make sure there is some L2 data to plot
      if(size(l2_phase,1) > 0)
        l2_phase_mean  = mean(phaseDDiffL2(l2_phase,(sv+1)));
        l2_phase_sigma = std(phaseDDiffL2(l2_phase,(sv+1)));

        subplot(NumFigs,1,2);
        plot(phaseDDiffL2(:,1), phaseDDiffL2(:,(sv+1))./10);
        grid;
        title(['SV = ',num2str(sv),' L2 Carrier Phase DD \mu = ', ...
        num2str(l2_phase_mean / 10,'%.1f'),'mcycle; \sigma = ', ...
        num2str(l2_phase_sigma / 10,'%.1f'),'mcycle']);
        ylabel('L2 DD [mcycle]');
        if (max(abs(phaseDDiffL2(:,(sv+1))./10)) > CarrThres_AScale)
          AScale_GPSPhase_L2(sv) = 1;
        end
        
        % If there is no L5 put the x-axis on the L2 graph
        if(NumFigs == 2)
          xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
        end
      
        % Contribute to the figure that plots all data superimposed
        % figure(FigAll(PHASE));
        set(0, 'CurrentFigure', FigAll(PHASE));
        subplot(NumFigs,1,2);
        plot(phaseDDiffL2(:,1), phaseDDiffL2(:,(sv+1))./10,'color', ...
            colour_LUT(colour,:));
        hold on;
        % figure(f(num));
        set(0, 'CurrentFigure', f(num));
      
        Phase_Sigma(L2,sv) = l2_phase_sigma/10;
        Phase_Mean (L2,sv) = l2_phase_mean/10;
  
        tmp = max(phaseDDiffL2(l2_phase,(sv+1)));
        if(tmp > max_l2_phase)
          max_l2_phase = tmp;
        end
        tmp = min(phaseDDiffL2(l2_phase,(sv+1)));
        if(tmp < min_l2_phase)
          min_l2_phase = tmp;
        end
        L2Available(num) = TRUE;
      else
        l2_phase_mean  = NaN;
        l2_phase_sigma = NaN;
        L2Available(num) = FALSE;
      end % End of L2 processing loop

      % Make sure there is some L5 data to plot
      if(size(l5_phase,1) > 0)
        l5_phase_mean  = mean(phaseDDiffL5(l5_phase,(sv+1)));
        l5_phase_sigma = std(phaseDDiffL5(l5_phase,(sv+1)));

        % L5 is the last data so even if there is no L2 this should
        % work!
        subplot(NumFigs,1,NumFigs);
        plot(phaseDDiffL5(:,1), phaseDDiffL5(:,(sv+1))./10);
        grid;
        title(['SV = ',num2str(sv),' L5 Carrier Phase DD \mu = ', ...
        num2str(l5_phase_mean / 10,'%.1f'),'mcycle; \sigma = ', ...
        num2str(l5_phase_sigma / 10,'%.1f'),'mcycle']);
        ylabel('L5 DD [mcycle]');
        
        xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');

        if (max(abs(phaseDDiffL5(:,(sv+1))./10)) > CarrThres_AScale)
          AScale_GPSPhase_L5(sv) = 1;
        end
      
        % Contribute to the figure that plots all data superimposed
        % figure(FigAll(PHASE));
        set(0, 'CurrentFigure', FigAll(PHASE));
        subplot(NumFigs,1,NumFigs);
        plot(phaseDDiffL5(:,1), phaseDDiffL5(:,(sv+1))./10,'color', ...
            colour_LUT(colour,:));
        hold on;
        % figure(f(num));
        set(0, 'CurrentFigure', f(num));
      
        Phase_Sigma(L5,sv) = l5_phase_sigma/10;
        Phase_Mean (L5,sv) = l5_phase_mean/10;
  
        tmp = max(phaseDDiffL5(l5_phase,(sv+1)));
        if(tmp > max_l5_phase)
          max_l5_phase = tmp;
        end
        tmp = min(phaseDDiffL5(l5_phase,(sv+1)));
        if(tmp < min_l5_phase)
          min_l5_phase = tmp;
        end
        L5Available(num) = TRUE;
      else
        l5_phase_mean  = NaN;
        l5_phase_sigma = NaN;
        L5Available(num) = FALSE;
      end % End of L5 processing loop

    else
      % No L1 or L2 data available
      l1_phase_mean  = NaN;
      l1_phase_sigma = NaN;
      l2_phase_mean  = NaN;
      l2_phase_sigma = NaN;
      l5_phase_mean  = NaN;
      l5_phase_sigma = NaN;
    end
  end  
  
  % Now go through the phase plots and set the axis up
  % so they are all the same.
  
  max_l1_phase =  ceil(max_l1_phase/10000) * 1000;
  min_l1_phase = floor(min_l1_phase/10000) * 1000;
  max_l2_phase =  ceil(max_l2_phase/10000) * 1000;
  min_l2_phase = floor(min_l2_phase/10000) * 1000;
  max_l5_phase =  ceil(max_l5_phase/10000) * 1000;
  min_l5_phase = floor(min_l5_phase/10000) * 1000;
  
  for i=3:num
    % figure(f(i));
    set(0, 'CurrentFigure', f(i));
    subplot(NumFigs,1,1);
    if (AScale_GPSPhase_L1(Sat(i)) == 0)
      axis([start stop -CarrMax CarrMax]);
    else
      % The max DD has significantly exceeded the axis limits
      % autoscale the Y-axis
      axisinfo = axis;
      axis([start stop axisinfo(3) axisinfo(4)]);
    end

    % L2 present is TRUE if there is a least one common L2 observable in
    % the base and rover files
    if((L2Present == TRUE) & (L2Available(i) == TRUE))
      % L2Available is TRUE if there is common data for the current plot.
      subplot(NumFigs,1,2);
      if (AScale_GPSPhase_L2(Sat(i)) == 0)
        axis([start stop -CarrMax CarrMax]);
      else
        % The max DD has significantly exceeded the axis limits
        % autoscale the Y-axis
        axisinfo = axis;
        axis([start stop axisinfo(3) axisinfo(4)]);
      end
    end
    
    % Now check for L5 data.
    if((L5Present == TRUE) & (L5Available(i) == TRUE))
      subplot(NumFigs,1,NumFigs);
      if (AScale_GPSPhase_L5(Sat(i)) == 0)
        axis([start stop -CarrMax CarrMax]);
      else
        % The max DD has significantly exceeded the axis limits
        % autoscale the Y-axis
        axisinfo = axis;
        axis([start stop axisinfo(3) axisinfo(4)]);
      end
    end

    labelAxes([],[],[],[],1);
    if(PlotPNG == TRUE)
      print('-dpng',['PRN',num2str(Sat(i)),'_PHASE.png']);
    end
  end
  
  % Close all the phase plots except the one where all satellites are
  % plotted on the same graph.
  for loop = 3:num
    close(f(loop));
  end
  
  num = 2;
  
  disp('Completed GPS Phase Double Differences...');
  
  max_l1_range = 0;
  min_l1_range = 0;
  max_l2_range = 0;
  min_l2_range = 0;
  max_l5_range = 0;
  min_l5_range = 0;
  num_range= num + 1;
  colour = 0;
  
  AScale_GPSRange_L1 = zeros(32,1);
  AScale_GPSRange_L2 = zeros(32,1);
  AScale_GPSRange_L5 = zeros(32,1);
  
  %
  % Loop around for each potential satellite (PRN #) and extract the L1 &
  % L2 double difference range information. Plot this and compute the
  % statistics.
  %
  for sv=1:32
    l1_range = find(isfinite(rangeDDiffL1(:,(sv+1))));
  
    if(L2Present == TRUE)
      l2_range = find(isfinite(rangeDDiffL2(:,(sv+1))));
    else
      l2_range = [];
    end
    
    if(L5Present == TRUE)
      l5_range = find(isfinite(rangeDDiffL5(:,(sv+1))));
    else
      l5_range = [];
    end
  
    if((size(l1_range,1) > 0) | (size(l2_range,1) > 0) | (size(l5_range,1) > 0))
      num = num + 1;
      f(num) = figure('Visible','off');
      Sat(num) = sv;
      colour = colour + 1;
  
      % Make sure there is some L1 data
      if(size(l1_range,1) > 0)
        l1_range_mean  = mean(rangeDDiffL1(l1_range,(sv+1)));
        l1_range_epochs = size(rangeDDiffL1(l1_range,(sv+1)),1);
        l1_range_sigma = std(rangeDDiffL1(l1_range,(sv+1)));
        subplot(NumFigs,1,1);
        plot(rangeDDiffL1(:,1), rangeDDiffL1(:,(sv+1))./10000);
        grid;
        title(['SV = ',num2str(sv),' L1 Pseudorange DD \mu = ', ...
        num2str(l1_range_mean / 100,'%.2f'),'cm; \sigma = ', ...
        num2str(l1_range_sigma / 100,'%.2f'),'cm']);
        ylabel('L1 DD [M]');
        if (max(abs(rangeDDiffL1(:,(sv+1))./10000)) > CodeThres_AScale)
          AScale_GPSRange_L1(sv) = 1;
        end
  
        % If there is no L2 data put the label under the L1 axis, otherwise
        % it gets placed at the bottom of the plot under the L2.
        if(NumFigs == 1)
          xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
        end
  
        % Contribute to the figure that plots all data superimposed
        % figure(FigAll(RANGE));
        set(0, 'CurrentFigure', FigAll(RANGE));        
        subplot(NumFigs,1,1);
        plot(rangeDDiffL1(:,1), rangeDDiffL1(:,(sv+1))./10000,'color', ...
           colour_LUT(colour,:));
        hold on;
        % figure(f(num));
        set(0, 'CurrentFigure', f(num));

        Range_Sigma(L1,sv) = l1_range_sigma/100;
        Range_Mean (L1,sv) = l1_range_mean/100;
        NumEpochs(L1,sv)   = l1_range_epochs;    
  
        tmp = max(rangeDDiffL1(l1_range,(sv+1)));
        if(tmp > max_l1_range)
          max_l1_range = tmp;
        end
        tmp = min(rangeDDiffL1(l1_range,(sv+1)));
        if(tmp < min_l1_range)
          min_l1_range = tmp;
        end
      else
        l1_range_mean  = NaN;
        l1_range_sigma = NaN;
      end
  
      % Make sure there is some L2 data
      if(size(l2_range,1) > 0)
        l2_range_mean   = mean(rangeDDiffL2(l2_range,(sv+1)));
        l2_range_epochs = size(rangeDDiffL2(l2_range,(sv+1)),1);
        l2_range_sigma  = std(rangeDDiffL2(l2_range,(sv+1)));
        subplot(NumFigs,1,2);
        plot(rangeDDiffL2(:,1), rangeDDiffL2(:,(sv+1))./10000);
        grid;
        title(['SV = ',num2str(sv),' L2 Pseudorange DD \mu = ', ...
        num2str(l2_range_mean / 100,'%.2f'),'cm; \sigma = ', ...
        num2str(l2_range_sigma / 100,'%.2f'),'cm']);
        ylabel('L2 DD [M]');
        if (max(abs(rangeDDiffL2(:,(sv+1))./10000)) > CodeThres_AScale)
          AScale_GPSRange_L2(sv) = 1;
        end

        % If there is no L5 put the x-axis on the L2 graph
        if(NumFigs == 2)
          xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
        end

        % Contribute to the figure that plots all data superimposed
        % figure(FigAll(RANGE));
        set(0, 'CurrentFigure', FigAll(RANGE));
        subplot(NumFigs,1,2);
        plot(rangeDDiffL2(:,1), rangeDDiffL2(:,(sv+1))./10000,'color', ...
           colour_LUT(colour,:));
        hold on;
        % figure(f(num));
        set(0, 'CurrentFigure', f(num));
  
        Range_Sigma(L2,sv) = l2_range_sigma/100;
        Range_Mean (L2,sv) = l2_range_mean/100;
        NumEpochs(L2,sv)   = l2_range_epochs;    
  
        tmp = max(rangeDDiffL2(l2_range,(sv+1)));
        if(tmp > max_l2_range)
          max_l2_range = tmp;
        end
        tmp = min(rangeDDiffL2(l2_range,(sv+1)));
        if(tmp < min_l2_range)
          min_l2_range = tmp;
        end
        L2Available(num) = TRUE;
      else
        l2_range_mean  = NaN;
        l2_range_sigma = NaN;
        L2Available(num) = FALSE;
      end % End of L2 range processing
      
      % Make sure there is some L5 data
      if(size(l5_range,1) > 0)
        l5_range_mean   = mean(rangeDDiffL5(l5_range,(sv+1)));
        l5_range_epochs = size(rangeDDiffL5(l5_range,(sv+1)),1);
        l5_range_sigma  = std(rangeDDiffL5(l5_range,(sv+1)));
        subplot(NumFigs,1,NumFigs);
        plot(rangeDDiffL5(:,1), rangeDDiffL5(:,(sv+1))./10000);
        grid;
        title(['SV = ',num2str(sv),' L5 Pseudorange DD \mu = ', ...
        num2str(l5_range_mean / 100,'%.2f'),'cm; \sigma = ', ...
        num2str(l5_range_sigma / 100,'%.2f'),'cm']);
        ylabel('L5 DD [M]');
        xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
        if (max(abs(rangeDDiffL5(:,(sv+1))./10000)) > CodeThres_AScale)
          AScale_GPSRange_L5(sv) = 1;
        end
  
        % Contribute to the figure that plots all data superimposed
        % figure(FigAll(RANGE));
        set(0, 'CurrentFigure', FigAll(RANGE));
        subplot(NumFigs,1,NumFigs);
        plot(rangeDDiffL5(:,1), rangeDDiffL5(:,(sv+1))./10000,'color', ...
           colour_LUT(colour,:));
        hold on;
        % figure(f(num));
        set(0, 'CurrentFigure', f(num));
  
        Range_Sigma(L5,sv) = l5_range_sigma/100;
        Range_Mean (L5,sv) = l5_range_mean/100;
        NumEpochs(L5,sv)   = l5_range_epochs;    
  
        tmp = max(rangeDDiffL5(l5_range,(sv+1)));
        if(tmp > max_l5_range)
          max_l5_range = tmp;
        end
        tmp = min(rangeDDiffL5(l5_range,(sv+1)));
        if(tmp < min_l5_range)
          min_l5_range = tmp;
        end
        L5Available(num) = TRUE;
      else
        l5_range_mean  = NaN;
        l5_range_sigma = NaN;
        L5Available(num) = FALSE;
      end % End of L5 range processing
    else
      % No L1 or L2 data available
      l1_range_mean  = NaN;
      l1_range_sigma = NaN;
      l2_range_mean  = NaN;
      l2_range_sigma = NaN;
      l5_range_mean  = NaN;
      l5_range_sigma = NaN;
    end
  end   
  
  % Now go through the phase plots and set the axis up
  % so they are all the same.
  
  max_l1_range =  ceil(max_l1_range/10000);
  min_l1_range = floor(min_l1_range/10000);
  if(L2Present == TRUE)
    max_l2_range =  ceil(max_l2_range/10000);
    min_l2_range = floor(min_l2_range/10000);
  end
  if(L5Present == TRUE)
    max_l5_range =  ceil(max_l5_range/10000);
    min_l5_range = floor(min_l5_range/10000);
  end
  
  for i=num_range:num
    % figure(f(i))
    set(0, 'CurrentFigure', f(i));
  
    subplot(NumFigs,1,1)
    if (AScale_GPSRange_L1(Sat(i)) == 0)
      axis([start stop -CodeMax CodeMax]);
    else
      % The max DD has significantly exceeded the axis limits
      % autoscale the Y-axis
      axisinfo = axis;
      axis([start stop axisinfo(3) axisinfo(4)]);
    end
    % L2 present is TRUE if there is a least one common L2 observable in
    % the base and rover files
    if((L2Present == TRUE) & (L2Available(i) == TRUE))
      % L2Available is TRUE if there is common data for the current plot.
      subplot(NumFigs,1,2)
      if (AScale_GPSRange_L2(Sat(i)) == 0)
        axis([start stop -CodeMax CodeMax]);
      else
        % The max DD has significantly exceeded the axis limits
        % autoscale the Y-axis
        axisinfo = axis;
        axis([start stop axisinfo(3) axisinfo(4)]);
      end
    end
    
    % Now test for L5
    if((L5Present == TRUE) & (L5Available(i) == TRUE))
      subplot(NumFigs,1,NumFigs)
      if (AScale_GPSRange_L5(Sat(i)) == 0)
        axis([start stop -CodeMax CodeMax]);
      else
        % The max DD has significantly exceeded the axis limits
        % autoscale the Y-axis
        axisinfo = axis;
        axis([start stop axisinfo(3) axisinfo(4)]);
      end
    end
    labelAxes([],[],[],[],1);
    if(PlotPNG == TRUE)
      print('-dpng',['PRN',num2str(Sat(i)),'_RANGE.png']);
    end
  end
  disp('Completed GPS Range Double Differences...');
  
  % Close all the range plots except the one where all satellites are
  % plotted on the same graph.
  for loop = num_range:num
    close(f(loop));
  end
  
  num = 2;


  num_range= num + 1;
  colour = 0;
  
  % DEBUG
  if(GLONASSInstalled)    
   
    %
    % GLONASS Phase processing
    %
    
    max_l1_phase = 0;
    min_l1_phase = 0;
    max_l2_phase = 0;
    min_l2_phase = 0;
    max_l1_p_phase = 0;
    min_l1_p_phase = 0;

    % GPS is plotted in millicycles, GLONASS in mm. 1mm ~= 5 millicycles.
    % To give similar scales divide the requested maximum by 5.
  
    CarrMaxGLN = CarrMax / 5;
    CarrThresGLN_AScale = CarrThres_AScale / 5;
    
    AScale_GLNPhase_L1CA = zeros(24,1);
    AScale_GLNPhase_L2   = zeros(24,1);
    AScale_GLNPhase_L1P  = zeros(24,1);
    
    for sv=1:24
      l1_phase = find(isfinite(phaseDDiffGL1_CA(:,(sv+1))));

      if(GL2_P_present == TRUE)
        l2_phase = find(isfinite(phaseDDiffGL2(:,(sv+1))));
      else
        l2_phase = [];
      end
    
      if(GL1_P_present == TRUE)
        l1_p_phase = find(isfinite(phaseDDiffGL1_P(:,(sv+1))));
      else
        l1_p_phase = [];
      end
  
      if((size(l1_phase,1) > 0) | (size(l2_phase,1) > 0) | (size(l1_p_phase,1) > 0))
        num = num + 1;
        f(num) = figure('Visible','off');
        set(0, 'CurrentFigure', f(num));
        Sat(num) = sv;
        colour = colour + 1;
  
        % Make sure there is some L1 data to plot.
        if(size(l1_phase,1) > 0)
          l1_phase_mean  = mean(phaseDDiffGL1_CA(l1_phase,(sv+1)));
          l1_phase_sigma = std(phaseDDiffGL1_CA(l1_phase,(sv+1)));
          subplot(NumGLNFigs,1,1);
          plot(phaseDDiffGL1_CA(:,1), phaseDDiffGL1_CA(:,(sv+1))./10);
          grid;
          title(['GLN SV = ',num2str(sv),' L1-C/A Carrier Phase DD \mu = ', ...
          num2str(l1_phase_mean / 10,'%.1f'),'mm; \sigma = ', ...
          num2str(l1_phase_sigma / 10,'%.1f'),'mm']);
          ylabel('L1-C/A DD [mm]');
          if (max(abs(phaseDDiffGL1_CA(:,(sv+1))./10)) > CarrThresGLN_AScale)
            AScale_GLNPhase_L1CA(sv) = 1;
          end

          % If there is no L2 data put the label under the L1 axis, otherwise
          % it gets placed at the bottom of the plot under the L2.
          if(NumGLNFigs == 1)
            xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
          end
          
          % Contribute to the figure that plots all data superimposed
          % figure(FigAll(PHASE));
          set(0, 'CurrentFigure', FigAllGln(PHASE));
          subplot(NumGLNFigs,1,1);
          plot(phaseDDiffGL1_CA(:,1), phaseDDiffGL1_CA(:,(sv+1))./10,'color', ...
            colour_LUT(colour,:));
          hold on;
          % figure(f(num));
          set(0, 'CurrentFigure', f(num));
          
          GLN_Phase_Sigma(L1_CA,sv) = l1_phase_sigma/10;
          GLN_Phase_Mean (L1_CA,sv) = l1_phase_mean/10;
      
          tmp = max(phaseDDiffGL1_CA(l1_phase,(sv+1)));
          if(tmp > max_l1_phase)
            max_l1_phase = tmp;
          end
          tmp = min(phaseDDiffGL1_CA(l1_phase,(sv+1)));
          if(tmp < min_l1_phase)
            min_l1_phase = tmp;
          end
        else
          l1_phase_mean  = NaN;
          l1_phase_sigma = NaN;
        end  % End of L1 processing loop
  
        % Make sure there is some L2 data to plot
        if(size(l2_phase,1) > 0)
          l2_phase_mean  = mean(phaseDDiffGL2(l2_phase,(sv+1)));
          l2_phase_sigma = std(phaseDDiffGL2(l2_phase,(sv+1)));

          subplot(NumGLNFigs,1,2);
          plot(phaseDDiffGL2(:,1), phaseDDiffGL2(:,(sv+1))./10);
          grid;
          title(['GLN SV = ',num2str(sv),' L2 Carrier Phase DD \mu = ', ...
          num2str(l2_phase_mean / 10,'%.1f'),'mm; \sigma = ', ...
          num2str(l2_phase_sigma / 10,'%.1f'),'mm']);
          ylabel('L2 DD [mm]');
          if (max(abs(phaseDDiffGL2(:,(sv+1))./10)) > CarrThresGLN_AScale)
            AScale_GLNPhase_L2(sv) = 1;
          end
        
          % If there is no L1-_ put the x-axis on the L2 graph
          if(NumGLNFigs == 2)
            xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
          end

          % Contribute to the figure that plots all data superimposed
          % figure(FigAll(PHASE));
          set(0, 'CurrentFigure', FigAllGln(PHASE));
          subplot(NumGLNFigs,1,2);
          plot(phaseDDiffGL2(:,1), phaseDDiffGL2(:,(sv+1))./10,'color', ...
            colour_LUT(colour,:));
          hold on;
          % figure(f(num));
          set(0, 'CurrentFigure', f(num));
                    
          GLN_Phase_Sigma(L2,sv) = l2_phase_sigma/10;
          GLN_Phase_Mean (L2,sv) = l2_phase_mean/10;
  
          tmp = max(phaseDDiffGL2(l2_phase,(sv+1)));
          if(tmp > max_l2_phase)
            max_l2_phase = tmp;
          end
          tmp = min(phaseDDiffGL2(l2_phase,(sv+1)));
          if(tmp < min_l2_phase)
            min_l2_phase = tmp;
          end
          L2Available(num) = TRUE;
        else
          l2_phase_mean  = NaN;
          l2_phase_sigma = NaN;
          L2Available(num) = FALSE;
        end % End of L2 processing loop

        % Make sure there is some L1-P data to plot
        if(size(l1_p_phase,1) > 0)
          l1_p_phase_mean  = mean(phaseDDiffGL1_P(l1_p_phase,(sv+1)));
          l1_p_phase_sigma = std(phaseDDiffGL1_P(l1_p_phase,(sv+1)));

          subplot(NumGLNFigs,1,NumGLNFigs);

          plot(phaseDDiffGL1_P(:,1), phaseDDiffGL1_P(:,(sv+1))./10);
          grid;
          title(['GLN SV = ',num2str(sv),' L1-P Carrier Phase DD \mu = ', ...
          num2str(l1_p_phase_mean / 10,'%.1f'),'mm; \sigma = ', ...
          num2str(l1_p_phase_sigma / 10,'%.1f'),'mm']);
          ylabel('L1-P DD [mm]');
          if (max(abs(phaseDDiffGL1_P(:,(sv+1))./10)) > CarrThresGLN_AScale)
            AScale_GLNPhase_L1P(sv) = 1;
          end
          fprintf('%d %f %f\n', sv, max(abs(phaseDDiffGL1_P(:,(sv+1))./10)), CarrThresGLN_AScale);
        
          xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
          
          % Contribute to the figure that plots all data superimposed
          % figure(FigAll(PHASE));
          set(0, 'CurrentFigure', FigAllGln(PHASE));
          subplot(NumGLNFigs,1,NumGLNFigs);
          plot(phaseDDiffGL1_P(:,1), phaseDDiffGL1_P(:,(sv+1))./10,'color', ...
            colour_LUT(colour,:));
          hold on;
          % figure(f(num));
          set(0, 'CurrentFigure', f(num));
      
          GLN_Phase_Sigma(L1_P,sv) = l1_p_phase_sigma/10;
          GLN_Phase_Mean (L1_P,sv) = l1_p_phase_mean/10;
  
          tmp = max(phaseDDiffGL1_P(l1_p_phase,(sv+1)));
          if(tmp > max_l1_p_phase)
            max_l1_p_phase = tmp;
          end
          tmp = min(phaseDDiffGL1_P(l1_p_phase,(sv+1)));
          if(tmp < min_l1_p_phase)
            min_l1_p_phase = tmp;
          end
          L1PAvailable(num) = TRUE;
        else
          l1_p_phase_mean  = NaN;
          l1_p_phase_sigma = NaN;
          L1PAvailable(num) = FALSE;
        end % End of L1P processing loop
      else
        % No L1 or L2 data available
        l1_phase_mean  = NaN;
        l1_phase_sigma = NaN;
        l2_phase_mean  = NaN;
        l2_phase_sigma = NaN;
        l1_p_phase_mean  = NaN;
        l1_p_phase_sigma = NaN;
      end
    end  
    
    % Now go through the phase plots and set the axis up
    % so they are all the same.
  
    max_l1_phase =  ceil(max_l1_phase/10000) * 1000;
    min_l1_phase = floor(min_l1_phase/10000) * 1000;
    if(GL2_P_present == TRUE)
      max_l2_phase =  ceil(max_l2_phase/10000) * 1000;
      min_l2_phase = floor(min_l2_phase/10000) * 1000;
    end
    if(GL1_P_present == TRUE)
      max_l1_p_phase =  ceil(max_l1_p_phase/10000) * 1000;
      min_l1_p_phase = floor(min_l1_p_phase/10000) * 1000;
    end
  
    for i=num_range:num
     % figure(f(i))
     set(0, 'CurrentFigure', f(i));
  
      subplot(NumGLNFigs,1,1)
      if (AScale_GLNPhase_L1CA(Sat(i)) == 0)
        axis([start stop -CarrMaxGLN CarrMaxGLN]);
      else
        % The max DD has significantly exceeded the axis limits
        % autoscale the Y-axis
        axisinfo = axis;
        axis([start stop axisinfo(3) axisinfo(4)]);
      end
      % L2 present is TRUE if there is a least one common L2 observable in
      % the base and rover files
      if((GL2_P_present == TRUE) & (L2Available(i) == TRUE))
        % L2Available is TRUE if there is common data for the current plot.
        subplot(NumGLNFigs,1,2)
        if (AScale_GLNPhase_L2(Sat(i)) == 0)
          axis([start stop -CarrMaxGLN CarrMaxGLN]);
        else
          % The max DD has significantly exceeded the axis limits
          % autoscale the Y-axis
          axisinfo = axis;
          axis([start stop axisinfo(3) axisinfo(4)]);
        end
      end
    
      % Now test for L1-P
      if((GL1_P_present == TRUE) & (L1PAvailable(i) == TRUE))
        subplot(NumGLNFigs,1,NumGLNFigs)
        if (AScale_GLNPhase_L1P(Sat(i)) == 0)
          axis([start stop -CarrMaxGLN CarrMaxGLN]);
        else
          % The max DD has significantly exceeded the axis limits
          % autoscale the Y-axis
          axisinfo = axis;
          axis([start stop axisinfo(3) axisinfo(4)]);
        end
      end
      labelAxes([],[],[],[],1);
      if(PlotPNG == TRUE)
        print('-dpng',['GLN',num2str(Sat(i)),'_PHASE.png']);
      end
    end
    disp('Completed GLONASS Phase Double Differences...');

    % Close all the phase plots except the one where all satellites are
    % plotted on the same graph.
    for loop = num_range:num
      close(f(loop));
    end


    %
    % GLONASS Range processing
    % 
    
    max_l1_range = 0;
    min_l1_range = 0;
    max_l2_range = 0;
    min_l2_range = 0;
    max_l1_p_range = 0;
    min_l1_p_range = 0;
    num = 2;
    num_range = num + 1;
    colour = 0;
    
    AScale_GLNRange_L1CA = zeros(24,1);
    AScale_GLNRange_L2   = zeros(24,1);
    AScale_GLNRange_L1P  = zeros(24,1);
    %
    % Loop around for each potential satellite (GLN #) and extract the L1 &
    % L2 double difference range information. Plot this and compute the
    % statistics.
    %
    for sv=1:24
      l1_range = find(isfinite(rangeDDiffGL1_CA(:,(sv+1))));
  
      if(GL2_P_present == TRUE)
        l2_range = find(isfinite(rangeDDiffGL2(:,(sv+1))));
      else
        l2_range = [];
      end
    
      if(GL1_P_present == TRUE)
        l1_p_range = find(isfinite(rangeDDiffGL1_P(:,(sv+1))));
      else
        l1_p_range = [];
      end
  
      if((size(l1_range,1) > 0) | (size(l2_range,1) > 0) | (size(l1_p_range,1) > 0))
        num = num + 1;
        f(num) = figure('Visible','off');
        Sat(num) = sv;
        colour = colour + 1;
  
        % Make sure there is some L1 data
        if(size(l1_range,1) > 0)
          l1_range_mean  = mean(rangeDDiffGL1_CA(l1_range,(sv+1)));
          l1_range_epochs = size(rangeDDiffGL1_CA(l1_range,(sv+1)),1);
          l1_range_sigma = std(rangeDDiffGL1_CA(l1_range,(sv+1)));
          subplot(NumGLNFigs,1,1);
          plot(rangeDDiffGL1_CA(:,1), rangeDDiffGL1_CA(:,(sv+1))./10000);
          grid;
          title(['GLN SV = ',num2str(sv),' L1-C/A Pseudorange DD \mu = ', ...
          num2str(l1_range_mean / 100,'%.2f'),'cm; \sigma = ', ...
          num2str(l1_range_sigma / 100,'%.2f'),'cm']);
          ylabel('L1-C/A DD [M]');
          if (max(abs(rangeDDiffGL1_CA(:,(sv+1))./10000)) > CodeThres_AScale)
            AScale_GLNRange_L1CA(sv) = 1;
          end
        
          % If there is no L2 data put the label under the L1 axis, otherwise
          % it gets placed at the bottom of the plot under the L2.
          if(NumGLNFigs == 1)
            xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
          end
  
          % Contribute to the figure that plots all data superimposed
          % figure(FigAll(RANGE));
          set(0, 'CurrentFigure', FigAllGln(RANGE));
          subplot(NumGLNFigs,1,1);
          plot(rangeDDiffGL1_CA(:,1), rangeDDiffGL1_CA(:,(sv+1))./10000,'color', ...
            colour_LUT(colour,:));
          hold on;
          % figure(f(num));
          set(0, 'CurrentFigure', f(num));
  
          GLN_Range_Sigma(L1_CA,sv) = l1_range_sigma/100;
          GLN_Range_Mean (L1_CA,sv) = l1_range_mean/100;
          GLN_NumEpochs(L1_CA,sv)   = l1_range_epochs;    
  
          tmp = max(rangeDDiffGL1_CA(l1_range,(sv+1)));
          if(tmp > max_l1_range)
            max_l1_range = tmp;
          end
          tmp = min(rangeDDiffGL1_CA(l1_range,(sv+1)));
          if(tmp < min_l1_range)
            min_l1_range = tmp;
          end
        else
          l1_range_mean  = NaN;
          l1_range_sigma = NaN;
        end
  
        % Make sure there is some L2 data
        if(size(l2_range,1) > 0)
          l2_range_mean   = mean(rangeDDiffGL2(l2_range,(sv+1)));
          l2_range_epochs = size(rangeDDiffGL2(l2_range,(sv+1)),1);
          l2_range_sigma  = std(rangeDDiffGL2(l2_range,(sv+1)));
          subplot(NumGLNFigs,1,2);
          plot(rangeDDiffGL2(:,1), rangeDDiffGL2(:,(sv+1))./10000);
          grid;
          title(['GLN SV = ',num2str(sv),' L2 Pseudorange DD \mu = ', ...
          num2str(l2_range_mean / 100,'%.2f'),'cm; \sigma = ', ...
          num2str(l2_range_sigma / 100,'%.2f'),'cm']);
          ylabel('L2 DD [M]');
 
          % If there is no L1-P put the x-axis on the L2 graph
          if(NumGLNFigs == 2)
            xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
          end

          % Contribute to the figure that plots all data superimposed
          % figure(FigAll(RANGE));
          set(0, 'CurrentFigure', FigAllGln(RANGE));
          subplot(NumGLNFigs,1,2);
          plot(rangeDDiffGL2(:,1), rangeDDiffGL2(:,(sv+1))./10000,'color', ...
            colour_LUT(colour,:));
          hold on;
          % figure(f(num));
          set(0, 'CurrentFigure', f(num));
          
          if (max(abs(rangeDDiffGL2(:,(sv+1))./10000)) > CodeThres_AScale)
            AScale_GLNRange_L2(sv) = 1;
          end
  
          GLN_Range_Sigma(L2,sv) = l2_range_sigma/100;
          GLN_Range_Mean (L2,sv) = l2_range_mean/100;
          GLN_NumEpochs(L2,sv)   = l2_range_epochs;    
  
          tmp = max(rangeDDiffGL2(l2_range,(sv+1)));
          if(tmp > max_l2_range)
            max_l2_range = tmp;
          end
          tmp = min(rangeDDiffGL2(l2_range,(sv+1)));
          if(tmp < min_l2_range)
            min_l2_range = tmp;
          end
          L2Available(num) = TRUE;
        else
          l2_range_mean  = NaN;
          l2_range_sigma = NaN;
          L2Available(num) = FALSE;
        end % End of L2 range processing
      
        % Make sure there is some L1-P data
        if(size(l1_p_range,1) > 0)
          l1_p_range_mean   = mean(rangeDDiffGL1_P(l1_p_range,(sv+1)));
          l1_p_range_epochs = size(rangeDDiffGL1_P(l1_p_range,(sv+1)),1);
          l1_p_range_sigma  = std(rangeDDiffGL1_P(l1_p_range,(sv+1)));
          subplot(NumGLNFigs,1,NumGLNFigs);
          plot(rangeDDiffGL1_P(:,1), rangeDDiffGL1_P(:,(sv+1))./10000);
          grid;
          title(['GLN SV = ',num2str(sv),' L1-P Pseudorange DD \mu = ', ...
          num2str(l1_p_range_mean / 100,'%.2f'),'cm; \sigma = ', ...
          num2str(l1_p_range_sigma / 100,'%.2f'),'cm']);
          ylabel('L1-P DD [M]');
          xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
          
          if (max(abs(rangeDDiffGL1_P(:,(sv+1))./10000)) > CodeThres_AScale)
            AScale_GLNRange_L1P(sv) = 1;
          end
  
          % Contribute to the figure that plots all data superimposed
          % figure(FigAll(RANGE));
          set(0, 'CurrentFigure', FigAllGln(RANGE));
          subplot(NumGLNFigs,1,NumGLNFigs);
          plot(rangeDDiffGL1_P(:,1), rangeDDiffGL1_P(:,(sv+1))./10000,'color', ...
            colour_LUT(colour,:));
          hold on;
          % figure(f(num));
          set(0, 'CurrentFigure', f(num));
  
          GLN_Range_Sigma(L1_P,sv) = l1_p_range_sigma/100;
          GLN_Range_Mean (L1_P,sv) = l1_p_range_mean/100;
          GLN_NumEpochs(L1_P,sv)   = l1_p_range_epochs;    
  
          tmp = max(rangeDDiffGL1_P(l1_p_range,(sv+1)));
          if(tmp > max_l1_p_range)
            max_l1_p_range = tmp;
          end
          tmp = min(rangeDDiffGL1_P(l1_p_range,(sv+1)));
          if(tmp < min_l1_p_range)
            min_l1_p_range = tmp;
          end
          L1PAvailable(num) = TRUE;
        else
          l1_p_range_mean  = NaN;
          l1_p_range_sigma = NaN;
          L1PAvailable(num) = FALSE;
        end % End of L5 range processing
      else
        % No L1 or L2 data available
        l1_range_mean  = NaN;
        l1_range_sigma = NaN;
        l2_range_mean  = NaN;
        l2_range_sigma = NaN;
        l1_p_range_mean  = NaN;
        l1_p_range_sigma = NaN;
      end
    end   
  
    % Now go through the phase plots and set the axis up
    % so they are all the same.
  
    max_l1_range =  ceil(max_l1_range/10000);
    min_l1_range = floor(min_l1_range/10000);
    if(GL2_P_present == TRUE)
      max_l2_range =  ceil(max_l2_range/10000);
      min_l2_range = floor(min_l2_range/10000);
    end
    if(GL1_P_present == TRUE)
      max_l1_p_range =  ceil(max_l1_p_range/10000);
      min_l1_p_range = floor(min_l1_p_range/10000);
    end
  
    for i=num_range:num
       % figure(f(i))
       set(0, 'CurrentFigure', f(i));
  
      subplot(NumGLNFigs,1,1)
      if (AScale_GLNRange_L1CA(Sat(i)) == 0)
        axis([start stop -CodeMax CodeMax]);
      else
        % The max DD has significantly exceeded the axis limits
        % autoscale the Y-axis
        axisinfo = axis;
        axis([start stop axisinfo(3) axisinfo(4)]);
      end
      % L2 present is TRUE if there is a least one common L2 observable in
      % the base and rover files
      if((GL2_P_present == TRUE) & (L2Available(i) == TRUE))
        % L2Available is TRUE if there is common data for the current plot.
        subplot(NumGLNFigs,1,2)
        if (AScale_GLNRange_L2(Sat(i)) == 0)
          axis([start stop -CodeMax CodeMax]);
        else
          % The max DD has significantly exceeded the axis limits
          % autoscale the Y-axis
          axisinfo = axis;
          axis([start stop axisinfo(3) axisinfo(4)]);
        end
      end
    
      % Now test for L1-P
      if((GL1_P_present == TRUE) & (L1PAvailable(i) == TRUE))
        subplot(NumGLNFigs,1,NumGLNFigs)
        if (AScale_GLNRange_L1P(Sat(i)) == 0)
          axis([start stop -CodeMax CodeMax]);
        else
          % The max DD has significantly exceeded the axis limits
          % autoscale the Y-axis
          axisinfo = axis;
          axis([start stop axisinfo(3) axisinfo(4)]);
        end
      end
      labelAxes([],[],[],[],1);
      if(PlotPNG == TRUE)
        print('-dpng',['GLN',num2str(Sat(i)),'_RANGE.png']);
      end
    end
    disp('Completed GLONASS Range Double Differences...');

    % Close all the range plots except the one where all satellites are
    % plotted on the same graph.
    for loop = num_range:num
      close(f(loop));
    end

  end % if(GLONASSInstalled)    
  
  disp('Processing Plots with all range/phase data included...');
  
  %
  % Now put some labels on the figures which contain data for
  % all the satellites
  % 
  % figure(FigAll(PHASE));
  set(0, 'CurrentFigure', FigAll(PHASE));
  subplot(NumFigs,1,NumFigs);
  axisinfo = axis;
  axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
  grid;
  xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
  if(L5Present == TRUE)
    % L5 is present
    ylabel('L5 DD [mcycle]');
    % Is L2 also present?
    if(L2Present == TRUE)
      subplot(NumFigs,1,2);
      ylabel('L2 DD [mcycle]');
      axisinfo = axis;
      axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
      grid;
    end
  elseif(L2Present == TRUE)
    % No L5 only L2
    ylabel('L2 DD [mcycle]');
  end
  subplot(NumFigs,1,1);
  axisinfo = axis;
  axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
  grid
  ylabel('L1 DD [mcycle]');
  title('DD Phase - Data from all satellites');
  labelAxes([],[],[],[],1);
  if(PlotPNG == TRUE)
    print('-dpng',['ALLSVS_PHASE_ASCALE.png']);
  end  
  % Now sort out the Y-axis limits to those the user specified
  subplot(NumFigs,1,NumFigs);
  axis([start stop -CarrMax CarrMax]);
  if(L5Present == TRUE)
    if(L2Present == TRUE)
      subplot(NumFigs,1,2);
      axis([start stop -CarrMax CarrMax]);
    end
  end
  subplot(NumFigs,1,1); 
  axis([start stop -CarrMax CarrMax]);
  if(PlotPNG == TRUE)
    print('-dpng',['ALLSVS_PHASE.png']);
  end  
  
  % figure(FigAll(RANGE));
  set(0, 'CurrentFigure', FigAll(RANGE));
  subplot(NumFigs,1,NumFigs);
  axisinfo = axis;
  axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
  grid;
  xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
  if(L5Present == TRUE)
    % L5 is present
    ylabel('L5 DD [M]');
    % Is L2 also present?
    if(L2Present == TRUE)
      subplot(NumFigs,1,2);
      ylabel('L2 DD [M]');
      axisinfo = axis;
      axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
      grid;
    end
  elseif(L2Present == TRUE)
    % No L5 only L2
    ylabel('L2 DD [M]');
  end
  subplot(NumFigs,1,1);
  axisinfo = axis;
  axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
  grid;
  ylabel('L1 DD [M]');
  title('DD Range - Data from all satellites');
  labelAxes([],[],[],[],1);
  if(PlotPNG == TRUE)
    print('-dpng',['ALLSVS_RANGE_ASCALE.png']);
  end  
  % Now sort out the Y-axis limits to those the user specified
  subplot(NumFigs,1,NumFigs);
  axis([start stop -CodeMax CodeMax]);
  if(L5Present == TRUE)
    if(L2Present == TRUE)
      subplot(NumFigs,1,2);
      axis([start stop -CodeMax CodeMax]);
    end
  end
  subplot(NumFigs,1,1); 
  axis([start stop -CodeMax CodeMax]);
  if(PlotPNG == TRUE)
    print('-dpng',['ALLSVS_RANGE.png']);
  end  
  
  if (GLONASSInstalled)
    % figure(FigAll(PHASE));
    set(0, 'CurrentFigure', FigAllGln(PHASE));
    subplot(NumGLNFigs,1,NumGLNFigs);
    axisinfo = axis;
    axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
    grid;
    xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
    if(GL1_P_present == TRUE)
      % L1P is present
      ylabel('L1-P DD [mm]');
      % Is L2 also present?
      if(GL2_P_present == TRUE)
        subplot(NumGLNFigs,1,2);
        ylabel('L2-P DD [mm]');
        axisinfo = axis;
        axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
        grid;
      end
    elseif(GL2_P_present == TRUE)
      % No L5 only L2
      ylabel('L2-P DD [mm]');
    end
    subplot(NumGLNFigs,1,1);
    axisinfo = axis;
    axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
    grid
    ylabel('L1-C/A DD [mm]');
    title('GLONASS DD Phase - Data from all satellites');
    labelAxes([],[],[],[],1);
    if(PlotPNG == TRUE)
      print('-dpng',['ALLSVS_PHASE_GLN_ASCALE.png']);
    end  
    % Now sort out the Y-axis limits to those the user specified
    subplot(NumGLNFigs,1,NumGLNFigs);
    axis([start stop -CarrMaxGLN CarrMaxGLN]);
    if(GL1_P_present == TRUE)
      if(GL2_P_present == TRUE)
        subplot(NumGLNFigs,1,2);
        axis([start stop -CarrMaxGLN CarrMaxGLN]);
      end
    end
    subplot(NumGLNFigs,1,1);
    axis([start stop -CarrMaxGLN CarrMaxGLN]);
    if(PlotPNG == TRUE)
      print('-dpng',['ALLSVS_PHASE_GLN.png']);
    end  
  
    % figure(FigAll(RANGE));
    set(0, 'CurrentFigure', FigAllGln(RANGE));
    subplot(NumGLNFigs,1,NumGLNFigs);
    axisinfo = axis;
    axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
    grid;
    xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
    if(GL1_P_present == TRUE)
      % L1P is present
      ylabel('L1-P DD [m]');
      % Is L2 also present?
      if(GL2_P_present == TRUE)
        subplot(NumGLNFigs,1,2);
        ylabel('L2-P DD [m]');
        axisinfo = axis;
        axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
        grid;
      end
    elseif(GL2_P_present == TRUE)
      % No L5 only L2
      ylabel('L2-P DD [m]');
    end
    subplot(NumGLNFigs,1,1);
    axisinfo = axis;
    axis([start stop axisinfo(3) axisinfo(4)]); % allow autoscale of the Y-axis
    grid;
    ylabel('L1 DD [m]');
    title('GLONASS DD Range - Data from all satellites');
    labelAxes([],[],[],[],1);
    if(PlotPNG == TRUE)
      print('-dpng',['ALLSVS_RANGE_GLN_ASCALE.png']);
    end  
    % Now sort out the Y-axis limits
    subplot(NumGLNFigs,1,NumGLNFigs);
    axis([start stop -CodeMax CodeMax]);
    if(GL1_P_present == TRUE)
      if(GL2_P_present == TRUE)
        subplot(NumGLNFigs,1,2);
        axis([start stop -CodeMax CodeMax]);
      end
    end
    subplot(NumGLNFigs,1,1);
    axis([start stop -CodeMax CodeMax]);
    if(PlotPNG == TRUE)
      print('-dpng',['ALLSVS_RANGE_GLN.png']);
    end  
  end

  %
  % Completed the basic double difference analysis. Now output the results
  % into a space delimited text file and an HTML page.
  % 
  loadPiHeaderHTM(fid, ...
                  L2Present, ...
                  L5Present, ...
                  PlotPNG, ...
                  doACFPlot, ...
                  doFFTPlot, ...
                  SkipCNRTests, ...
                  GLONASSInstalled, ...
                  GL1_P_present, ...
                  GL2_P_present); % Write the HTM header
  
  %
  % Create figures for the elevation/azimuth plots.
  % To preserve some memory first remove all the figures.
  %
  close all;

  % For now only do L1/L2 ACF analysis, 
  % TODO Add L5
  if(doACFPlot == TRUE)
    % Create the ACF graphs if requested.
    Processed = do_Xcorr(phaseDDiffL1, phaseDDiffL2, ...
                         rangeDDiffL1, rangeDDiffL2, ...
                         PlotPNG, 1000);
    close all;
  end
  
  f(1) = figure('Visible','off');
  f(2) = figure('Visible','off');
  f(3) = figure('Visible','off');
  colour = 0;
  
  %
  % For each satellite compute the correlation coefficients of the double
  % differences and plot the FFT if requested. Then build up a plot of
  % elevation vs azimuth.
  %
  for i=1:32
    if( isfinite(Range_Mean(L1,i)) | ... 
        isfinite(Range_Mean(L2,i)) | ...
        isfinite(Range_Mean(L5,i)) )
  
      % Output to an HTML file
      fprintf(fid,'<TR><td align="center">%d</td>\n',i);
      fprintf(fid,'<td align="center">%d</td>\n', NumEpochs(L1,i));
  
      if(L2Present == TRUE)
        fprintf(fid,'<td align="center">%d</td>\n', NumEpochs(L2,i));
      end
      
      if(L5Present == TRUE)
        fprintf(fid,'<td align="center">%d</td>\n', NumEpochs(L5,i));
      end
      
      if(PlotPNG == TRUE)
        fprintf(fid,'<td align="center"><a href="PRN%d_RANGE.png">%.1f</a></td>\n',i,Range_Mean(L1,i));
      else
        fprintf(fid,'<td align="center">%.1f</a></td>\n',Range_Mean(L1,i));
      end
      
      fprintf(fid,'<td align="center">%.1f</td>\n', Range_Sigma(L1,i));

      if(L2Present == TRUE)
        fprintf(fid,'<td align="center">%.1f</td>\n', Range_Mean(L2,i));
        fprintf(fid,'<td align="center">%.1f</td>\n', Range_Sigma(L2,i));
      end
      
      if(L5Present == TRUE)
        fprintf(fid,'<td align="center">%.1f</td>\n', Range_Mean(L5,i));
        fprintf(fid,'<td align="center">%.1f</td>\n', Range_Sigma(L5,i));
      end
      
      if(PlotPNG == TRUE)
        fprintf(fid,'<td align="center"><a href="PRN%d_PHASE.png">%.1f | %.3f</a></td>\n', ...
                                    i, ...
                                    Phase_Mean(L1,i), ...
                                    Phase_Mean(L1,i)*L1_WAVELENGTH);
      else
        fprintf(fid,'<td align="center">%.1f | %.3f</a></td>\n', ...
                                      Phase_Mean(L1,i), ...
                                      Phase_Mean(L1,i)*L1_WAVELENGTH);
      end
  
  
      fprintf(fid,'<td align="center">%.1f | %.3f</td>\n',  ...
                                      Phase_Sigma(L1,i), ...
                                      Phase_Sigma(L1,i)*L1_WAVELENGTH);
  
      if(L2Present == TRUE)
        fprintf(fid,'<td align="center">%.1f | %.3f</td>\n', ...
                                      Phase_Mean(L2,i), ...
                                      Phase_Mean(L2,i)*L2_WAVELENGTH);
  
        fprintf(fid,'<td align="center">%.1f | %.3f</td>\n', ...
                                      Phase_Sigma(L2,i), ...
                                      Phase_Sigma(L2,i)*L2_WAVELENGTH);
      end
      
      if(L5Present == TRUE)
        fprintf(fid,'<td align="center">%.1f | %.3f</td>\n', ...
                                      Phase_Mean(L5,i), ...
                                      Phase_Mean(L5,i)*L5_WAVELENGTH);
  
        fprintf(fid,'<td align="center">%.1f | %.3f</td>\n', ...
                                      Phase_Sigma(L5,i), ...
                                      Phase_Sigma(L5,i)*L5_WAVELENGTH);
      end
      

      if( (L2Present == TRUE) & (L5Present == TRUE) )
        if(  (size(phaseDDiffL1,1) ~= size(phaseDDiffL2,1)) & (size(phaseDDiffL1,1) ~= size(phaseDDiffL5,1)) )
          % The L1, L2 & L5 vectors are different lengths. This
          % occasionally happens when multiple .dat files are
          % concatenated. Just kludge this by adding some null
          % data to the end of the shorter array.
  
          len_l1 = size(phaseDDiffL1,1);
          len_l2 = size(phaseDDiffL2,1);
          len_l5 = size(phaseDDiffL5,1);
% TODO Handle L5    
          if(len_l1 > len_l2)
            for a=(len_l2+1):len_l1
              phaseDDiffL2(a,1:129) = ones(1,129).*NaN;
              rangeDDiffL2(a,1:97)  = ones(1,97).*NaN;
            end
          else
            for a=(len_l1+1):len_l2
              phaseDDiffL1(a,1:161) = ones(1,161).*NaN;
              rangeDDiffL1(a,1:97)  = ones(1,97).*NaN;
            end
          end
        end

        % Form an array of DD code/carrier observables
        DD = [phaseDDiffL1(:,i+1), ...
              phaseDDiffL2(:,i+1), ...
              phaseDDiffL5(:,i+1), ...
              rangeDDiffL1(:,i+1), ...
              rangeDDiffL2(:,i+1), ...
              rangeDDiffL5(:,i+1)];

        ddL1 = 1;
        ddL2 = 2;
        ddL5 = 3; 
        
      elseif(L2Present == TRUE)
        if( size(phaseDDiffL1,1) ~= size(phaseDDiffL2,1) )
          % The L1 and L2 vectors are different lengths. This
          % occasionally happens when multiple .dat files are
          % concatenated. Just kludge this by adding some null
          % data to the end of the shorter array.
  
          len_l1 = size(phaseDDiffL1,1);
          len_l2 = size(phaseDDiffL2,1);
    
          if(len_l1 > len_l2)
            for a=(len_l2+1):len_l1
              phaseDDiffL2(a,1:129) = ones(1,129).*NaN;
              rangeDDiffL2(a,1:97)  = ones(1,97).*NaN;
            end
          else
            for a=(len_l1+1):len_l2
              phaseDDiffL1(a,1:161) = ones(1,161).*NaN;
              rangeDDiffL1(a,1:97)  = ones(1,97).*NaN;
            end
          end
        end
        % Form an array of DD code/carrier observables
        DD = [phaseDDiffL1(:,i+1), ...
              phaseDDiffL2(:,i+1), ...
              rangeDDiffL1(:,i+1), ...
              rangeDDiffL2(:,i+1)];

        ddL1 = 1;
        ddL2 = 2;
        ddL5 = NaN; 
            
      elseif(L5Present == TRUE)
        
        if( size(phaseDDiffL1,1) ~= size(phaseDDiffL5,1) )
          % The L1 and L5 vectors are different lengths. This
          % occasionally happens when multiple .dat files are
          % concatenated. Just kludge this by adding some null
          % data to the end of the shorter array.
  
          len_l1 = size(phaseDDiffL1,1);
          len_l5 = size(phaseDDiffL5,1);
    
          if(len_l1 > len_l5)
            for a=(len_l5+1):len_l1
              phaseDDiffL5(a,1:129) = ones(1,129).*NaN;
              rangeDDiffL5(a,1:97)  = ones(1,97).*NaN;
            end
          else
            for a=(len_l1+1):len_l5
              phaseDDiffL1(a,1:161) = ones(1,161).*NaN;
              rangeDDiffL1(a,1:97)  = ones(1,97).*NaN;
            end
          end
        end
        % Form an array of DD code/carrier observables
        DD = [phaseDDiffL1(:,i+1), ...
              phaseDDiffL5(:,i+1), ...
              rangeDDiffL1(:,i+1), ...
              rangeDDiffL5(:,i+1)];

        ddL1 = 1;
        ddL2 = NaN;
        ddL5 = 2; 

      else
        % Form an array of DD code/carrier observables
        DD = [phaseDDiffL1(:,i+1), ...
              rangeDDiffL1(:,i+1)];

        ddL1 = 1;
        ddL2 = NaN;
        ddL5 = NaN; 

      end

      DD = excise_matrix(DD);

      % Make sure there is a reasonable amount of data
      % before computing the correlation coefficients.
      if(size(DD,1) > 100)
        
        if(L2Present == TRUE)
          CC = corrcoef(DD(:,ddL1),DD(:,ddL2)); 
          fprintf(fid,'<td align="center">%.3f</td>\n', ...
                        CC(1,2));
                                      
          CC = corrcoef(DD(:,NumFigs+ddL1),DD(:,NumFigs+ddL2)); 
          fprintf(fid,'<td align="center">%.3f</td>\n', ...
                        CC(1,2));
        end            
                      
        CC = corrcoef(DD(:,ddL1),DD(:,NumFigs+ddL1));
        fprintf(fid,'<td align="center">%.3f</td>\n', ...
                        CC(1,2));

        if (L2Present == TRUE)
          CC = corrcoef(DD(:,ddL2),DD(:,NumFigs+ddL2)); 
          fprintf(fid,'<td align="center">%.3f</td>\n', ...
                        CC(1,2));
        end

        if(L5Present == TRUE)

          CC = corrcoef(DD(:,ddL1),DD(:,ddL5)); 
          fprintf(fid,'<td align="center">%.3f</td>\n', ...
                        CC(1,2));
                                      
          CC = corrcoef(DD(:,NumFigs+ddL1),DD(:,NumFigs+ddL5)); 
          fprintf(fid,'<td align="center">%.3f</td>\n', ...
                        CC(1,2)); 
                                      
          CC = corrcoef(DD(:,ddL5),DD(:,NumFigs+ddL5)); 
          fprintf(fid,'<td align="center">%.3f</td>\n', ...
                        CC(1,2));
        end
          
        if(doFFTPlot == 1)
          % Analyze the frequency spectra for multipath
          fftplot(i+1,phaseDDiffL1,phaseDDiffL2,rangeDDiffL1,rangeDDiffL2);
        end
  
        if( (doFFTPlot == 1) & (PlotPNG == TRUE) )
          fprintf(fid,'<td align="center"><a href="PRN%d_FREQ.png">graph</a></td>\n',i);
          print('-dpng',['PRN',num2str(i),'_FREQ.png']);
        end
        
        if( (doACFPlot == 1) & (PlotPNG == TRUE) & (Processed(i) == 1))
          fprintf(fid,'<td align="center"><a href="ACF_Pseudorange_%d.png">graph</a></td>\n',i);
          fprintf(fid,'<td align="center"><a href="ACF_Phase_%d.png">graph</a></td>\n',i);
        end

        fprintf(fid,'</TR>\n'); % Terminate this row of the table
  
      else
        fprintf(fid,'<td></td></TR>\n');
      end
  
      fprintf(fid,'\n');
  
      ftxt = fopen('stats.txt','w');
      % Output the data to a space delimited text file
      fprintf(ftxt,'%d %d %d %.1f %.1f %.1f %.1f %.1f %.1f %.1f %.1f\n',...
                 i, ...
                 NumEpochs(L1,i), ...
                 NumEpochs(L2,i), ...
                 Range_Mean(L1,i), ...
                 Range_Sigma(L1,i), ...
                 Range_Mean(L2,i), ...
                 Range_Sigma(L2,i), ...
                 Phase_Mean(L1,i), ...
                 Phase_Sigma(L1,i), ...
                 Phase_Mean(L2,i), ...
                 Phase_Sigma(L2,i));
       fclose(ftxt);
    end  
  
    %
    % The previous if statement checked for valid double differences. We
    % now check for valid elevation angles. If we do not do this, and
    % the reference satellite has been constant throughout the data set
    % we will not plot the elevation and azimuth data for the
    % reference, as there is no DD for the reference.
    %
    if(sum(isfinite(phaseDDiffL1(:,i+1+32))) > 0)
      %
      % Now contribute data to the elevation/azimuth plot
      %
      % Format: [Time Elevation Azimuth]
      % 
      el_az = [phaseDDiffL1(:,1) ...
               phaseDDiffL1(:,i+1+32) ...
               phaseDDiffL1(:,i+1+64)];
  
      colour = colour + 1;
      
      % Elevation/Azimuth polar plot
      % figure(f(1));
      set(0, 'CurrentFigure', f(1));      
      polar_expanded(el_az(:,3).*pi/180, 90 - el_az(:,2),'color', ...
                     colour_LUT(colour,:),'linewidth',2);
      hold on;
      
      % Elevation versus time
      % figure(f(2));
      set(0, 'CurrentFigure', f(2));      
      plot(el_az(:,1), el_az(:,2),'color',colour_LUT(colour,:), ...
           'linewidth',2);
      hold on;
     
      % Azimuth versus time
      % figure(f(3));
      set(0, 'CurrentFigure', f(3));            
      plot(el_az(:,1), el_az(:,3),'color',colour_LUT(colour,:), ...
           'linewidth',2);
      clear el_az;
      hold on;
  
      % Array of PRN numbers used later to create the legend.
      SvPRN(colour) = i;
    end
    outBuffer = sprintf('Completed correlation coefficient el/az data processing PRN %2d...',i);
    disp(outBuffer);
  end
  
  % End the table
  fprintf(fid,'</table></center></div><HR>\n');
 
  if(GLONASSInstalled)
    analyzeGLONASS(fid, ...
                   GL2_P_present, ...
                   GL1_P_present, ...
                   PlotPNG, ...
                   L1_CA, ...
                   L2, ...
                   L1_P, ...
                   GLN_Range_Mean, ...
                   GLN_Range_Sigma, ...
                   GLN_Phase_Mean, ...
                   GLN_Phase_Sigma, ...
                   GLN_NumEpochs, ...
                   FDMAChan);
  end                 
                        
  % figure(f(1));
  set(0, 'CurrentFigure', f(1));

  % Form the legend text
  for i=1:colour
    StringMatrix(i) = {['PRN ',int2str(SvPRN(i))]};
  end
  legend(StringMatrix,0);
  
  title('Azimuth Versus Elevation, Origin = Zenith')
  if(PlotPNG == TRUE)
    print('-dpng','ElevAzi.png');
  end  
  
  disp('Processing Elevation Data');
  % figure(f(2));
  set(0, 'CurrentFigure', f(2));
  legend(StringMatrix,0);
  title('Elevation versus Time');
  xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
  ylabel('Elevation [Degrees]');
  % Set the axis
  A    = axis;
  A(3) = 0;
  A(4) = 90;
  axis(A);
  grid
  labelAxes([],[],[],[],1);
  if(PlotPNG == TRUE)
    print('-dpng','Elev.png');
    disp('Completed Printing');
  end  
  
  
  disp('Processing azimuth Data');
  % figure(f(3));
  set(0, 'CurrentFigure', f(3));
  legend(StringMatrix,0);
  title('Azimuth versus Time');
  xlabel(['Time [GPS Sec]    ',timeStr],'Interpreter','none');
  ylabel('Azimuth [Degrees]');
  % Set the axis
  A    = axis;
  A(3) = 0;
  A(4) = 360;
  axis(A);
  grid
  labelAxes([],[],[],[],1);
  if(PlotPNG == TRUE)
    print('-dpng','Azi.png');
    disp('Completed Printing');
  end

  % Free up some memory by closing all open figures.
  close all
else
  if( (PlotPNG == TRUE) & (QuickMode < 2) )
    fprintf(fid,'<h1>QuickMode - PiLine Analysis Results</h1>\n');
    fprintf(fid,'<p><ul>\n');

    fprintf(fid,'  <li><a href="ElevCNRL1.png">GPS: Base-Station L1 C/No versus Elevation</a></li>\n');
    fprintf(fid,'  <li><a href="RoverElevCNRL1.png">GPS: Rover-Station L1 C/No versus Elevation</a></li>\n');
    fprintf(fid,'  <li><a href="ElevCNRL1Diff.png">GPS: Base minus Rover L1 C/No versus Elevation</a></li>\n');
    fprintf(fid,'  <li><a href="CNRL1Diff.png">GPS: Base minus Rover L1 C/No versus time</a></li>\n');
    
    if(L2Present == TRUE)
      fprintf(fid,'  <li><a href="ElevCNRL2.png">GPS: Base-Station L2 C/No versus Elevation</a></li>\n');
      fprintf(fid,'  <li><a href="RoverElevCNRL2.png">GPS: Rover-Station L2 C/No versus Elevation</a></li>\n');
      fprintf(fid,'  <li><a href="ElevCNRL2Diff.png">GPS: Base minus Rover L2 C/No versus Elevation</a></li>\n');
      fprintf(fid,'  <li><a href="CNRL2Diff.png">GPS: Base minus Rover L2 C/No versus time</a></li>\n');
    end
  
    if(L2Present == TRUE)
      fprintf(fid,'  <li><a href="ElevCNRL5.png">GPS: Base-Station L5 C/No versus Elevation</a></li>\n');
      fprintf(fid,'  <li><a href="RoverElevCNRL5.png">GPS: Rover-Station L5 C/No versus Elevation</a></li>\n');
      fprintf(fid,'  <li><a href="ElevCNRL5Diff.png">GPS: Base minus Rover L5 C/No versus Elevation</a></li>\n');
      fprintf(fid,'  <li><a href="CNRL5Diff.png">GPS: Base minus Rover L5 C/No versus time</a></li>\n');
    end

    fprintf(fid,'  <li><a href="CNRDiff.htm">GPS/GLN Tabulated Base Minus Rover C/No</a></li>\n');

    fprintf(fid,'  <li><a href="tracking.png">GPS: Satellite Tracking plot</a></li>\n');
    if(GLONASSInstalled)
      fprintf(fid,'  <li><a href="GLNtracking.png">GLN: Satellite Tracking plot</a></li>\n');
    end
    fprintf(fid,'  <li><a href="DDVersusElev.png">GPS: Double Difference Residuals versus elevation angle</a></li>\n');
    fprintf(fid,'  <li><a href="DDVersusCNo.png">GPS: Double Difference Residuals versus C/No</a></li>\n');
    fprintf(fid,'  <li><a href="L1PhaseDist.png">GPS: Double Difference L1 Phase distribution</a></li>\n');
    
    if(L2Present == TRUE)
      fprintf(fid,'  <li><a href="L2PhaseDist.png">GPS: Double Difference L2 Phase distribution</a></li>\n');
    end
    
    if(L5Present == TRUE)
      fprintf(fid,'  <li><a href="L5PhaseDist.png">GPS: Double Difference L5 Phase distribution</a></li>\n');
    end

    if(GLONASSInstalled)
      fprintf(fid,'  <li><a href="GLNL1CAPhaseDist.png">GLN: Double Difference L1-C/A Phase distribution</a></li>\n');
     
      if(GL2_P_present == TRUE)
        fprintf(fid,'  <li><a href="GLNL2PhaseDist.png">GLN: Double Difference L2 Phase distribution</a></li>\n');
      end
     
      if(GL1_P_present == TRUE)
        fprintf(fid,'  <li><a href="GLNL1PPhaseDist.png">GLN: Double Difference L1-P Phase distribution</a></li>\n');
      end
    end

    
    fprintf(fid,'  <li><a href="L1RangeDist.png">GPS: Double Difference L1 Range distribution</a></li>\n');
    
    if(L2Present == TRUE)
      fprintf(fid,'  <li><a href="L2RangeDist.png">GPS: Double Difference L2 Range distribution</a></li>\n');
    end
    
    if(L5Present == TRUE)
      fprintf(fid,'  <li><a href="L5RangeDist.png">GPS: Double Difference L5 Range distribution</a></li>\n');
    end
    
    if(GLONASSInstalled)
      fprintf(fid,'  <li><a href="GLNL1CARangeDist.png">GLN: Double Difference L1-C/A Range distribution</a></li>\n');
     
      if(GL2_P_present == TRUE)
        fprintf(fid,'  <li><a href="GLNL2RangeDist.png">GLN: Double Difference L2 Range distribution</a></li>\n');
      end
     
      if(GL1_P_present == TRUE)
        fprintf(fid,'  <li><a href="GLNL1PRangeDist.png">GLN: Double Difference L1-P Range distribution</a></li>\n');
      end
    end

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

disp('Making Temporary arrays');
% For ease seperate the data out, each coloumn of the new matrix
% is a single PRN (a possible 32 in total). A NaN means there is no data
% available at that epoch for the PRN number.
Elev  = phaseDDiffL1(:,34:65);
Azi   = phaseDDiffL1(:,66:97);
if(GLONASSInstalled)
  ElevL1CA = phaseDDiffGL1_CA(:,26:49);
  ElevL2G  = phaseDDiffGL2(:,26:49);
end

if(SkipCNRTests == 0)
  CNRL1 = phaseDDiffL1(:,98:129);
else  
  CNRL1 = [];
end

if(L2Present == TRUE)
  CNRL2 = phaseDDiffL2(:,66:97);
else
  CNRL2 = [];
end

if(L5Present == TRUE)
  CNRL5 = phaseDDiffL5(:,66:97);
else
  CNRL5 = [];
end

disp('Made arrays');

if((QuickMode == 0) & (SkipCNRTests == 0))

  binSize = 0.5 ; % 0.5dBHz steps

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Loop around for L1 and L2 and create a polar plot of the C/No against
  % elevation/azimuth. The colour map represents the C/No, the code to
  % manipulate this and replace the elevation labels (which get hidden
  % when there is a large amount of data) was provided by Timo Allison.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for figNumber=1:NumFigs
    fig        = figure('Visible','off');
    if(figNumber == L1)
      % 25-60dBHz (L1)
      % binMaxEdge = 60 ;
      % binMinEdge = 25 ;
      binMaxEdge = max(1 + round(max(CNRL1)));
      binMinEdge = min(round(min(CNRL1)));
      
      % Clip the min/max CNR
      % i = find(CNRL1 < binMinEdge);
      % if(length(i) > 0)
      %   CNRL1(i) = ones(length(i),1).*binMinEdge; 
      % end
      % i = find(CNRL1 > binMaxEdge);
      % if(length(i) > 0)
      %   CNRL1(i) = ones(length(i),1).*binMaxEdge; 
      % end
    elseif( (figNumber == L2) & (L2Present == TRUE) )
      % 0-55dBHz (L2)
      % binMaxEdge = 55 ;
      % binMinEdge = 0 ;
      binMaxEdge = max(1 + round(max(CNRL2)));
      binMinEdge = min(round(min(CNRL2)));
      
      % Clip the min/max CNR
      % i = find(CNRL2 < binMinEdge);
      % if(length(i) > 0)
      %   CNRL2(i) = ones(length(i),1).*binMinEdge; 
      % end
      % i = find(CNRL2 > binMaxEdge);
      % if(length(i) > 0)
      %   CNRL2(i) = ones(length(i),1).*binMaxEdge; 
      % end
    else
      % 0-55dBHz (L5)
      % binMaxEdge = 55 ;
      % binMinEdge = 0 ;
      binMaxEdge = max(1 + round(max(CNRL5)));
      binMinEdge = min(round(min(CNRL5)));
      
      % Clip the min/max CNR
      % i = find(CNRL5 < binMinEdge);
      % if(length(i) > 0)
      %   CNRL5(i) = ones(length(i),1).*binMinEdge; 
      % end
      % i = find(CNRL5 > binMaxEdge);
      % if(length(i) > 0)
      %   CNRL5(i) = ones(length(i),1).*binMaxEdge; 
      % end
    end
    % Compute the number of bins and set the colour map appropriately.
    binNumber  = 1 + ( binMaxEdge - binMinEdge ) / binSize ;  
    Map        = jet(binNumber);
    colormap( Map ) ;  % Set custom map.
    % Set a dummy point in black ('k') and hold to plot so
    % we can loop around and plot each CNR colour seperately.
    polar_expanded(0, 0, 'k');
    hold on;
  
    if(figNumber == L5)
      disp('Processing L5 C/No Versus Elev/Azimuth - Please Wait...');
    else
      outBuffer = sprintf('Processing L%d C/No Versus Elev/Azimuth - Please Wait...',figNumber);
      disp(outBuffer);
    end

    for i=binMinEdge:binSize:binMaxEdge  % i == CNR
    
      % Find data for the current CNR
      if(figNumber == L1)
        j = find( round(CNRL1./binSize)*binSize == i);
      elseif( (figNumber == L2) & (L2Present == TRUE) )
        j = find( round(CNRL2./binSize)*binSize == i);
      else
        j = find( round(CNRL5./binSize)*binSize == i);
      end
  
      if(length(j) > 0)
        % Add the data to the polar plot 
        polar_expanded(Azi(j).*pi/180, 90 - Elev(j),'x', 'color', ...
           Map( ( (i-binMinEdge(1)) / binSize) + 1,:),'markersize',4);
      end

      if(round(i) == i)
        if(figNumber == L5)
          outBuffer = sprintf('Processed L5 %2d dBHz ....',i);
        else
          outBuffer = sprintf('Processed L%d %2d dBHz ....',figNumber,i);
        end
        disp(outBuffer); % Avoid the "ans =" string by using disp.
      end
    end
  
    if(figNumber == 1)
      title('Rover L1 C/No as a function of Elevation/Azimuth');
    elseif( (figNumber == L2) & (L2Present == TRUE) )
      title('Rover L2 C/No as a function of Elevation/Azimuth');
    else
      title('Rover L5 C/No as a function of Elevation/Azimuth');
    end
    
    colorBarHandle = colorbar;
 
    % Get handle for the color bar axes and polar axes.
    childrenHandles = get( gcf, 'children' ) ;
    % Matlab R2006a includes uimenu handles along with axes menus. We only
    % care about axes handles. So filter handles of type axes. Filter out
    % uimenu handles.
    childrenHandles = ...
      childrenHandles( strcmp( get(childrenHandles, 'type'), 'axes' ) ) ;

    % Get polar axes handle. The order of the axes handles is not
    % guaranteed, and has changed between MatLab R12 and R13, so must
    % compare with the colorbar handle.
    if colorBarHandle == childrenHandles( 1 )
      polarAxesHandle = childrenHandles( 2 ) ;
    else
      polarAxesHandle = childrenHandles( 1 ) ;
    end
  
    % Fully customize the colour bar 5 dB steps
    Min = round(binMinEdge/5)*5;
    Max = round(binMaxEdge/5)*5;
    colorTickLabelNum = Min:5:Max;
    % binMinEdge is element 1 of the colormap, binMaxEdge is the
    % last element of the map. The ytick is the labeling we want
    % to use and is computed according to the offset in the map.
    yTick = ( colorTickLabelNum - binSize/2 - binMinEdge) / binSize + 1 ;
  
    index = 0 ;
    for labelNum = colorTickLabelNum
      % Build label cell array. In this case, the cell array could be set
      % as in the following line. The loop shows the generic method.
      % colorTickLabel { '50', '100', .... } ;
      index = index + 1 ;
      colorTickLabel( index ) = { num2str( labelNum ) } ;
    end
  
    set( colorBarHandle, 'ytick', yTick ) ;  % Label every second bin.
    set( colorBarHandle, 'yticklabel', colorTickLabel ) ;

    % axes( colorBarHandle ) ;
    set(gcf, 'CurrentAxes', colorBarHandle);
    title( 'dBHz' ) ;

    % The polar plot "tracks" have now obliterated the underlying
    % elevation labels (on the radius).
    % These need to be redraw - not a simple process. The handles of the
    % labels (text objects) (which have been hidden) must first be found.
    % Then the text positions and strings. The original text objects are
    % deleted, and new labels written with bold font.
  
    % Update root property to unhide polar plot text objects.
    set( 0, 'showhiddenhandles', 'on')

    % Get all children handles of the current axes object.
    h = get( polarAxesHandle, 'children') ;
  
    % Filter out  text objects.
    h = h( strcmp( get(h, 'type'), 'text' ) ) ;
    elevTextHandles = h( end-8: end ) ;  % Assumes 9 labels '0'...'80'.
  
    % Store original values for rewriting.
    stringCellArray   = get( elevTextHandles, 'string' ) ;
    positionCellArray = get( elevTextHandles, 'position' ) ;
  
    % axes( polarAxesHandle ) ;
    set(gcf, 'CurrentAxes', polarAxesHandle);
    for i=1:length( elevTextHandles )
      % Useful if decide to update font size.
      delete( elevTextHandles( i ) ) ; 
      % Each cell of positionCellArray = [x y 0]. (z set to 0).
      text( positionCellArray{ i }(1), positionCellArray{ i }(2), ...
            stringCellArray( i ), 'verticalalignment', 'bottom', ...
            'FontWeight', 'bold' ) ;
    end
  
    % Set root property back to default.
    set( 0, 'showhiddenhandles', 'off')
    % End of updates to elevation labels.

    if(PlotPNG == TRUE)
      if(figNumber == 1)
        print('-dpng',['ElevAziCNRL1.png']);
        disp('L1 C/No Plot created ...');
      elseif( (figNumber == L2) & (L2Present == TRUE) )
        print('-dpng',['ElevAziCNRL2.png']);
        disp('L2 C/No Plot created ...');
      else
        print('-dpng',['ElevAziCNRL5.png']);
        disp('L5 C/No Plot created ...');
      end  
    end
  end

  close all;
end


ElevL2 = phaseDDiffL2(:,34:65);
ElevL5 = phaseDDiffL5(:,34:65);
if(SkipCNRTests == 0)
  %********************************************************************
  % Now plot C/No versus elevation angle for the base and 
  % rover-stations L1 & L2 data.
  %********************************************************************
  CNRDiff = ones(4,32).*NaN; 

  for freq=1:NumFigs
    if(freq == 1)
      CNRL1Base = phaseDDiffL1(:,130:161);
    elseif( (freq == L2) & (L2Present == TRUE) )
      CNRL2Base = phaseDDiffL2(:,98:129);
    else
      CNRL5Base = phaseDDiffL5(:,98:129);
    end

    if(QuickMode < 3)
      if(freq == 1)
        %**************************************************************
        %             Plot Base L1 versus Elevation
        %**************************************************************
        plotCNR(CNRL1Base,Elev,dropLegend,32);
        if(PlotPNG == TRUE)
          ylabel('L1 C/No [dBHz]');
          title('GPS: Elevation Versus Base L1 C/No');
          print('-dpng',['ElevCNRL1.png']);
          disp('GPS: L1 Base C/No versus Elev Plot created ...');
          close;
        end
        %**************************************************************
        %             Plot Rover L1 versus Elevation
        %**************************************************************
        plotCNR(CNRL1,Elev,dropLegend,32);
        if(PlotPNG == TRUE)
          ylabel('L1 C/No [dBHz]');
          title('GPS: Elevation Versus Rover L1 C/No');
          print('-dpng',['RoverElevCNRL1.png']);
          disp('GPS: L1 Rover C/No versus Elev Plot created ...');
          close;
        end
      elseif( (freq == L2) & (L2Present == TRUE) )
        %**************************************************************
        %             Plot Base L2 versus Elevation
        %**************************************************************
        plotCNR(CNRL2Base,ElevL2,dropLegend,32);
        if(PlotPNG == TRUE)
          ylabel('L2 C/No [dBHz]');
          title('GPS: Elevation Versus Base L2 C/No');
          print('-dpng',['ElevCNRL2.png']);
          disp('GPS: L2 Base C/No versus Elev Plot created ...');
          close;
        end
        %**************************************************************
        %             Plot Rover L2 versus Elevation
        %**************************************************************
        plotCNR(CNRL2,ElevL2,dropLegend,32);
        if(PlotPNG == TRUE)
          ylabel('L2 C/No [dBHz]');
          title('GPS: Elevation Versus Rover L2 C/No');
          print('-dpng',['RoverElevCNRL2.png']);
          disp('GPS: L2 Rover C/No versus Elev Plot created ...');
          close;
        end
      else  
        %**************************************************************
        %             Plot Base L5 versus Elevation
        %**************************************************************
        plotCNR(CNRL5Base,ElevL5,dropLegend,32);
        if(PlotPNG == TRUE)
          ylabel('L5 C/No [dBHz]');
          title('GPS: Elevation Versus Base L5 C/No');
          print('-dpng',['ElevCNRL5.png']);
          disp('GPS: L5 Base C/No versus Elev Plot created ...');
          close;
        end
        %**************************************************************
        %             Plot Rover L5 versus Elevation
        %**************************************************************
        plotCNR(CNRL5,ElevL5,dropLegend,32);
        if(PlotPNG == TRUE)
          ylabel('L5 C/No [dBHz]');
          title('GPS: Elevation Versus Rover L5 C/No');
          print('-dpng',['RoverElevCNRL5.png']);
          disp('GPS: L5 Rover C/No versus Elev Plot created ...');
          close;
        end
      end
    end  

    %**************************************************************
    % Finally plot the difference between base and rover 
    %**************************************************************
    if(freq == 1)
      %************************************************************
      %            Plot the L1 difference versus Elevation
      %************************************************************
      PlotCNRDiffElev(CNRL1Base,CNRL1,Elev,dropLegend,32);
      ylabel('GPS: Base - Rover L1 C/No [dB]');
      title('GPS: Elevation Versus Base-Rover L1 C/No');
      if(PlotPNG == TRUE)
        print('-dpng',['ElevCNRL1Diff.png']);
        disp('GPS: L1 Base - Rover C/No versus Elev Plot created ...');
      end
      close
    elseif( (freq == L2) & (L2Present == TRUE) )
      %************************************************************
      %            Plot the L2 difference versus Elevation
      %************************************************************
      PlotCNRDiffElev(CNRL2Base,CNRL2,ElevL2,dropLegend,32);
      ylabel('GPS: Base - Rover L2 C/No [dB]');
      title('GPS: Elevation Versus Base-Rover L2 C/No');
      if(PlotPNG == TRUE)
        print('-dpng',['ElevCNRL2Diff.png']);
        disp('GPS: L2 Base - Rover C/No versus Elev Plot created ...');
      end
      close
    else
      %************************************************************
      %            Plot the L5 difference versus Elevation
      %************************************************************
      PlotCNRDiffElev(CNRL5Base,CNRL5,ElevL5,dropLegend,32);
      ylabel('GPS: Base - Rover L5 C/No [dB]');
      title('GPS: Elevation Versus Base-Rover L5 C/No');
      if(PlotPNG == TRUE)
        print('-dpng',['ElevCNRL5Diff.png']);
        disp('GPS: L5 Base - Rover C/No versus Elev Plot created ...');
      end
      close
    end

    if(freq == 1)
      %************************************************************
      %            Plot the L1 difference versus Time
      %************************************************************
      CNRDiffL1 = PlotCNRDiffTime(CNRL1Base,CNRL1,phaseDDiffL1,timeStr,CNRTestElMask,dropLegend,32);
      ylabel('GPS: Base - Rover L1 C/No [dB]');
      if(PlotPNG == TRUE)
        print('-dpng',['CNRL1Diff.png']);
        close
      end
      disp('GPS: L1 Base - Rover C/No Plot created ...');
      % Now free some memory
      clear CNRL1Base;
    elseif( (freq == L2) & (L2Present == TRUE) )
      %************************************************************
      %            Plot the L2 difference versus Time
      %************************************************************
      CNRDiffL2 = PlotCNRDiffTime(CNRL2Base,CNRL2,phaseDDiffL2,timeStr,CNRTestElMask,dropLegend,32);
      ylabel('GPS: Base - Rover L2 C/No [dB]');
      if(PlotPNG == TRUE)
        print('-dpng',['CNRL2Diff.png']);
        close
      end
      disp('GPS: L2 Base - Rover C/No Plot created ...');
      % Now free some memory
      clear CNRL2Base;
    else
      %************************************************************
      %            Plot the L5 difference versus Time
      %************************************************************
      CNRDiffL5 = PlotCNRDiffTime(CNRL5Base,CNRL5,phaseDDiffL5,timeStr,CNRTestElMask,dropLegend,32);
      ylabel('GPS: Base - Rover L5 C/No [dB]');
      if(PlotPNG == TRUE)
        print('-dpng',['CNRL5Diff.png']);
        close
      end
      disp('GPS: L5 Base - Rover C/No Plot created ...');
      % Now free some memory
      clear CNRL5Base;
    end  
  end

  %
  % GLONASS
  %
  if(GLONASSInstalled)

    %********************************************************************
    % Now plot C/No versus elevation angle for the base and 
    % rover-stations L1 & L2 data.
    %********************************************************************
    CNRDiff = ones(4,24).*NaN; 

    ElevL1CA = phaseDDiffGL1_CA(:,26:49);
    ElevL2G   = phaseDDiffGL2(:,26:49);
    ElevL1P  = phaseDDiffGL1_P(:,26:49);

    CNRL1CA  = phaseDDiffGL1_CA(:,74:97);
    CNRL2    = phaseDDiffGL2(:,50:73);
    CNRL1P   = phaseDDiffGL1_P(:,74:97);

    for freq=1:NumGLNFigs
      if(freq == 1)
        CNRL1CABase = phaseDDiffGL1_CA(:,98:121);
      elseif( (freq == L2) & (GL2_P_present == TRUE) )
        CNRL2Base = phaseDDiffGL2(:,74:97);
      else
        CNRL1PBase = phaseDDiffGL1_P(:,98:121);
      end
  
      if(QuickMode < 3)
        if(freq == 1)
          %**************************************************************
          %             Plot Base L1-C/A versus Elevation
          %**************************************************************
          plotCNR(CNRL1CABase,ElevL1CA,dropLegend,24);
          if(PlotPNG == TRUE)
            ylabel('GLN: L1-C/A C/No [dBHz]');
            title('GLN: Elevation Versus Base L1-C/A C/No');
            print('-dpng',['ElevGLNCNRL1_CA.png']);
            disp('GLN: L1-C/A Base C/No versus Elev Plot created ...');
            close;
          end
          %**************************************************************
          %             Plot Rover L1-C/A versus Elevation
          %**************************************************************
          plotCNR(CNRL1CA,ElevL1CA,dropLegend,24);
          if(PlotPNG == TRUE)
            ylabel('GLN: L1-C/A C/No [dBHz]');
            title('GLN: Elevation Versus Rover L1-C/A C/No');
            print('-dpng',['RoverElevGLNCNRL1_CA.png']);
            disp('GLN: L1-C/A Rover C/No versus Elev Plot created ...');
            close;
          end
        elseif( (freq == L2) & (GL2_P_present == TRUE) )
          %**************************************************************
          %             Plot Base L2 versus Elevation
          %**************************************************************
          plotCNR(CNRL2Base,ElevL2G,dropLegend,24);
          if(PlotPNG == TRUE)
            ylabel('GLN: L2 C/No [dBHz]');
            title('GLN: Elevation Versus Base L2 C/No');
            print('-dpng',['ElevGLNCNRL2.png']);
            disp('GLN: L2 Base C/No versus Elev Plot created ...');
            close;
          end
          %**************************************************************
          %             Plot Rover L2 versus Elevation
          %**************************************************************
          plotCNR(CNRL2,ElevL2G,dropLegend,24);
          if(PlotPNG == TRUE)
            ylabel('GLN: L2 C/No [dBHz]');
            title('GLN: Elevation Versus Rover L2 C/No');
            print('-dpng',['RoverElevGLNCNRL2.png']);
            disp('GLN: L2 Rover C/No versus Elev Plot created ...');
            close;
          end
        else  
          %**************************************************************
          %             Plot Base L1-P versus Elevation
          %**************************************************************
          plotCNR(CNRL1PBase,ElevL1P,dropLegend,24);
          if(PlotPNG == TRUE)
            ylabel('GLN: L1-P C/No [dBHz]');
            title('GLN: Elevation Versus Base L1-P C/No');
            print('-dpng',['ElevGLNCNRL1P.png']);
            disp('GLN: L1-P Base C/No versus Elev Plot created ...');
            close;
          end
          %**************************************************************
          %             Plot Rover L1-P versus Elevation
          %**************************************************************
          plotCNR(CNRL1P,ElevL1P,dropLegend,24);
          if(PlotPNG == TRUE)
            ylabel('GLN: L1-P C/No [dBHz]');
            title('GLN: Elevation Versus Rover L1-P C/No');
            print('-dpng',['RoverElevGLNCNRL1P.png']);
            disp('GLN: L1-P Rover C/No versus Elev Plot created ...');
            close;
          end
        end
      end  

      %**************************************************************
      % Finally plot the difference between base and rover 
      %**************************************************************
      if(freq == 1)
        %************************************************************
        %            Plot the L1-C/A difference versus Elevation
        %************************************************************
        PlotCNRDiffElev(CNRL1CABase,CNRL1CA,ElevL1CA,dropLegend,24);
        ylabel('GLN: Base - Rover L1-C/A C/No [dB]');
        title('GLN: Elevation Versus Base-Rover L1-C/A C/No');
        if(PlotPNG == TRUE)
          print('-dpng',['ElevGLNCNRL1CADiff.png']);
          disp('GLN: L1-C/A Base - Rover C/No versus Elev Plot created ...');
        end
        close
      elseif( (freq == L2) & (GL2_P_present == TRUE) )
        %************************************************************
        %            Plot the L2 difference versus Elevation
        %************************************************************
        PlotCNRDiffElev(CNRL2Base,CNRL2,ElevL2G,dropLegend,24);
        ylabel('GLN: Base - Rover L2 C/No [dB]');
        title('GLN: Elevation Versus Base-Rover L2 C/No');
        if(PlotPNG == TRUE)
          print('-dpng',['ElevGLNCNRL2Diff.png']);
          disp('GLN: L2 Base - Rover C/No versus Elev Plot created ...');
        end
        close
      else
        %************************************************************
        %            Plot the L1-P difference versus Elevation
        %************************************************************
        PlotCNRDiffElev(CNRL1PBase,CNRL1P,ElevL1P,dropLegend,24);
        ylabel('GLN: Base - Rover L1-P C/No [dB]');
        title('GLN: Elevation Versus Base-Rover L1-P C/No');
        if(PlotPNG == TRUE)
          print('-dpng',['ElevGLNCNRL1PDiff.png']);
          disp('GLN: L1-P Base - Rover C/No versus Elev Plot created ...');
        end
        close
      end

      if(freq == 1)
        %************************************************************
        %            Plot the L1-C/A difference versus Time
        %************************************************************
        GLNCNRDiffL1CA = PlotCNRDiffTime(CNRL1CABase,CNRL1CA,phaseDDiffGL1_CA,timeStr,CNRTestElMask,dropLegend,24);
        ylabel('GLN: Base - Rover L1-C/A C/No [dB]');
        if(PlotPNG == TRUE)
          print('-dpng',['GLNCNRL1CADiff.png']);
          close
        end
        disp('GLN: L1-C/A Base - Rover C/No Plot created ...');
        % Now free some memory
        clear CNRL1CABase;
      elseif( (freq == L2) & (GL2_P_present == TRUE) )
        %************************************************************
        %            Plot the L2 difference versus Time
        %************************************************************
        GLNCNRDiffL2 = PlotCNRDiffTime(CNRL2Base,CNRL2,phaseDDiffGL2,timeStr,CNRTestElMask,dropLegend,24);
        ylabel('GLN: Base - Rover L2 C/No [dB]');
        if(PlotPNG == TRUE)
          print('-dpng',['GLNCNRL2Diff.png']);
          close
        end
        disp('GLN: L2 Base - Rover C/No Plot created ...');
        % Now free some memory
        clear CNRL2Base;
      else
        %************************************************************
        %            Plot the L1-P difference versus Time
        %************************************************************
        GLNCNRDiffL1P = PlotCNRDiffTime(CNRL1PBase,CNRL1P,phaseDDiffGL1_P,timeStr,CNRTestElMask,dropLegend,24);
        ylabel('GLN: Base - Rover L1-P C/No [dB]');
        if(PlotPNG == TRUE)
          print('-dpng',['GLNCNRL1PDiff.png']);
          close
        end
        disp('GLN: L1-P Base - Rover C/No Plot created ...');
        % Now free some memory
        clear CNRL1PBase;
      end  
    end
  end % if(GLONASSInstalled)

  fDiff     = fopen('CNRDiff.txt','w');
  if( (QuickMode == 2) | (QuickMode == 3) )
    fDiffHTML = fopen('stats.htm','w');
  else
    fDiffHTML = fopen('CNRDiff.htm','w');
  end

  % Define the table
  fprintf(fDiffHTML,'<HTML><HEAD><title>Base - Rover C/No difference Results</title></head>\n');
  fprintf(fDiffHTML,'<BODY background="../backgrnd.gif">');
  fprintf(fDiffHTML,'<img src="logo.png">\n');
  fprintf(fDiffHTML,'<font face="arial">\n');
  fprintf(fDiffHTML,'<hr>');
  fprintf(fDiffHTML,'<h1>Base Minus Rover C/No</h1>');
  fprintf(fDiffHTML,'<ul>\n');

  if(QuickMode < 3)
    fprintf(fDiffHTML,'  <li><a href="ElevCNRL1.png">GPS: Base-Station L1 C/No versus Elevation</a></li>\n');
    fprintf(fDiffHTML,'  <li><a href="RoverElevCNRL1.png">GPS: Rover-Station L1 C/No versus Elevation</a></li>\n');
    
    if(GLONASSInstalled)
      fprintf(fDiffHTML,'  <li><a href="ElevGLNCNRL1_CA.png">GLN: Base-Station L1-C/A C/No versus Elevation</a></li>\n');
      fprintf(fDiffHTML,'  <li><a href="RoverElevGLNCNRL1_CA.png">GLN: Rover-Station L1-C/A C/No versus Elevation</a></li>\n');
    end
  end  
  fprintf(fDiffHTML,'  <li><a href="ElevCNRL1Diff.png">GPS: Base minus Rover L1 C/No versus Elevation</a></li>\n');
  fprintf(fDiffHTML,'  <li><a href="CNRL1Diff.png">GPS: Base minus Rover L1 C/No versus time</a></li>\n');
  
  if(GLONASSInstalled)
    fprintf(fDiffHTML,'  <li><a href="ElevGLNCNRL1CADiff.png">GLN: Base minus Rover L1-C/A C/No versus Elevation</a></li>\n');
    fprintf(fDiffHTML,'  <li><a href="GLNCNRL1CADiff.png">GLN: Base minus Rover L1-C/A C/No versus time</a></li>\n');
  end


  if( (GL1_P_present == TRUE)  & (GLONASSInstalled) )
    if(QuickMode < 3)
      fprintf(fDiffHTML,'  <li><a href="ElevGLNCNRL1P.png">GLN: Base-Station L1-P C/No versus Elevation</a></li>\n');
      fprintf(fDiffHTML,'  <li><a href="RoverElevGLNCNRL1P.png">GLN: Rover-Station L1-P C/No versus Elevation</a></li>\n');
    end
    fprintf(fDiffHTML,'  <li><a href="ElevGLNCNRL1PDiff.png">GLN: Base minus Rover L1-P C/No versus Elevation</a></li>\n');
    fprintf(fDiffHTML,'  <li><a href="GLNCNRL1PDiff.png">GLN: Base minus Rover L1-P C/No versus time</a></li>\n');
  end

  if ( (QuickMode < 3) && (L2Present == TRUE) )
    fprintf(fDiffHTML,'  <li><a href="ElevCNRL2.png">GPS: Base-Station L2 C/No versus Elevation</a></li>\n');
    fprintf(fDiffHTML,'  <li><a href="RoverElevCNRL2.png">GPS: Rover-Station L2 C/No versus Elevation</a></li>\n');
  end
  if ( (QuickMode < 3) && GLONASSInstalled && GL2_P_present )
    fprintf(fDiffHTML,'  <li><a href="ElevGLNCNRL2.png">GLN: Base-Station L2 C/No versus Elevation</a></li>\n');
    fprintf(fDiffHTML,'  <li><a href="RoverElevGLNCNRL2.png">GLN: Rover-Station L2 C/No versus Elevation</a></li>\n');
  end

  if (L2Present == TRUE)
    fprintf(fDiffHTML,'  <li><a href="ElevCNRL2Diff.png">GPS: Base minus Rover L2 C/No versus Elevation</a></li>\n');
    fprintf(fDiffHTML,'  <li><a href="CNRL2Diff.png">GPS: Base minus Rover L2 C/No versus time</a></li>\n');
  end
  if(GLONASSInstalled && GL2_P_present)
      fprintf(fDiffHTML,'  <li><a href="ElevGLNCNRL2Diff.png">GLN: Base minus Rover L2 C/No versus Elevation</a></li>\n');
      fprintf(fDiffHTML,'  <li><a href="GLNCNRL2Diff.png">GLN: Base minus Rover L2 C/No versus time</a></li>\n');
  end
    
  if(L5Present == TRUE)
    if(QuickMode < 3)
      fprintf(fDiffHTML,'  <li><a href="ElevCNRL5.png">GPS: Base-Station L5 C/No versus Elevation</a></li>\n');
      fprintf(fDiffHTML,'  <li><a href="RoverElevCNRL5.png">GPS: Rover-Station L5 C/No versus Elevation</a></li>\n');
    end
    fprintf(fDiffHTML,'  <li><a href="ElevCNRL5Diff.png">GPS: Base minus Rover L5 C/No versus Elevation</a></li>\n');
    fprintf(fDiffHTML,'  <li><a href="CNRL5Diff.png">GPS: Base minus Rover L5 C/No versus time</a></li>\n');
  end

  fprintf(fDiffHTML,'</ul><br><br>\n');
  fprintf(fDiffHTML,'<center>');
  fprintf(fDiffHTML,'<table border="1">');
  fprintf(fDiffHTML,'<tr>\n');
  
  % Write out the table header
  fprintf(fDiffHTML,'<th>SV #</th>\n');
  fprintf(fDiffHTML,'<th>System</th>\n');
  fprintf(fDiffHTML,'<TH># L1 Epochs</TH>\n');
  
  if(L2Present == TRUE)
    fprintf(fDiffHTML,'<TH># L2 Epochs</TH>\n');
  end
  
  if(L5Present == TRUE)
    fprintf(fDiffHTML,'<TH># L5 Epochs</TH>\n');
  end

  fprintf(fDiffHTML,'<TH>L1 Mean Difference [dB]</TH>\n');
  if(L2Present == TRUE)
    fprintf(fDiffHTML,'<TH>L2 Mean Difference [dB]</TH>\n');
  end
 
  if(L5Present == TRUE)
    fprintf(fDiffHTML,'<TH>L5 Mean Difference [dB]</TH>\n');
  end
  
  fprintf(fDiffHTML,'<TH>L1 Sigma [dB]</TH>\n');
  
  if(L2Present == TRUE)
    fprintf(fDiffHTML,'<TH>L2 Sigma [dB]</TH>\n');
  end
  
  if(L5Present == TRUE)
    fprintf(fDiffHTML,'<TH>L5 Sigma [dB]</TH>\n');
  end
  
  fprintf(fDiffHTML,'</tr>\n');


  % First process GPS
  for SV=1:32
    if(isfinite(CNRDiffL1(1,SV)))
      % Output the data to a space delimited text file

      fprintf(fDiffHTML,'<tr>\n');
      if( (L2Present == TRUE) & (L5Present == TRUE) )
        fprintf(fDiff,'%d 0 %d %d %d %.3f %.3f %.3f %.3f %.3f %.3f\n', ...
                   SV, ...
                   CNRDiffL1(2,SV), ...
                   CNRDiffL2(2,SV), ...
                   CNRDiffL5(2,SV), ...
                   CNRDiffL1(1,SV), ...
                   CNRDiffL2(1,SV), ...
                   CNRDiffL5(1,SV), ...
                   CNRDiffL1(3,SV), ...
                   CNRDiffL2(3,SV), ...
                   CNRDiffL5(3,SV));
        fprintf(fDiffHTML,'<td>%d</td><td>GPS</td>\n <td>%d</td>\n <td>%d</td>\n <td>%d</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n  <td>%.3f</td>\n', ...
                   SV, ...
                   CNRDiffL1(2,SV), ...
                   CNRDiffL2(2,SV), ...
                   CNRDiffL5(2,SV), ...
                   CNRDiffL1(1,SV), ...
                   CNRDiffL2(1,SV), ...
                   CNRDiffL5(1,SV), ...
                   CNRDiffL1(3,SV), ...
                   CNRDiffL2(3,SV), ...
                   CNRDiffL5(3,SV));
      elseif(L2Present == TRUE)
        fprintf(fDiff,'%d 0 %d %d %.3f %.3f %.3f %.3f\n', ...
                   SV, ...
                   CNRDiffL1(2,SV), ...
                   CNRDiffL2(2,SV), ...
                   CNRDiffL1(1,SV), ...
                   CNRDiffL2(1,SV), ...
                   CNRDiffL1(3,SV), ...
                   CNRDiffL2(3,SV));
        fprintf(fDiffHTML,'<td>%d</td><td>GPS</td>\n <td>%d</td>\n <td>%d</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n  <td>%.3f</td>\n', ...
                   SV, ...
                   CNRDiffL1(2,SV), ...
                   CNRDiffL2(2,SV), ...
                   CNRDiffL1(1,SV), ...
                   CNRDiffL2(1,SV), ...
                   CNRDiffL1(3,SV), ...
                   CNRDiffL2(3,SV));
      elseif(L5Present == TRUE)
        fprintf(fDiff,'%d 0 %d %d %.3f %.3f %.3f %.3f\n', ...
                   SV, ...
                   CNRDiffL1(2,SV), ...
                   CNRDiffL5(2,SV), ...
                   CNRDiffL1(1,SV), ...
                   CNRDiffL5(1,SV), ...
                   CNRDiffL1(3,SV), ...
                   CNRDiffL5(3,SV));
        fprintf(fDiffHTML,'<td>%d</td><td>GPS</td>\n <td>%d</td>\n <td>%d</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n  <td>%.3f</td>\n', ...
                   SV, ...
                   CNRDiffL1(2,SV), ...
                   CNRDiffL5(2,SV), ...
                   CNRDiffL1(1,SV), ...
                   CNRDiffL5(1,SV), ...
                   CNRDiffL1(3,SV), ...
                   CNRDiffL5(3,SV));
      else
        fprintf(fDiff,'%d 0 %d %.3f %.3f\n', ...
                   SV, ...
                   CNRDiffL1(2,SV), ...
                   CNRDiffL1(1,SV), ...
                   CNRDiffL1(3,SV));
        fprintf(fDiffHTML,'<td>%d</td><td>GPS</td>\n <td>%d</td>\n <td>%.3f</td>\n <td>%.3f</td>\n', ...
                   SV, ...
                   CNRDiffL1(2,SV), ...
                   CNRDiffL1(1,SV), ...
                   CNRDiffL1(3,SV));
      end
      fprintf(fDiffHTML,'</tr>\n');
    end
  end



  %********************************************************************
  %      Now form the weighted mean L1 and L2 differences 
  %********************************************************************

  i = find(isfinite(CNRDiffL1(2,:)) & (CNRDiffL1(3,:) < CNRTestSigmaMask ) );
  NumL1Epochs = sum(CNRDiffL1(2,i)); 
  MeanL1Diff  = sum(CNRDiffL1(2,i).*CNRDiffL1(1,i)) / NumL1Epochs;
  MeanL1Sigma = mean(CNRDiffL1(3,i));

  if(L2Present == TRUE)
    i = find(isfinite(CNRDiffL2(2,:)) & (CNRDiffL2(3,:) < CNRTestSigmaMask ) );
    NumL2Epochs = sum(CNRDiffL2(2,i));
    MeanL2Diff  = sum(CNRDiffL2(2,i).*CNRDiffL2(1,i)) / NumL2Epochs;
    MeanL2Sigma = mean(CNRDiffL2(3,i));
  end

  if(L5Present == TRUE)
    i = find(isfinite(CNRDiffL5(2,:)) & (CNRDiffL5(3,:) < CNRTestSigmaMask ) );
    NumL5Epochs = sum(CNRDiffL5(2,i));
    MeanL5Diff  = sum(CNRDiffL5(2,i).*CNRDiffL5(1,i)) / NumL5Epochs;
    MeanL5Sigma = mean(CNRDiffL5(3,i));
  end

  if(L2Present == TRUE)  & (L5Present == TRUE)
    fprintf(fDiffHTML,'<th>Total</th><th>GPS</th>\n <th>%d</th>\n <th>%d</th>\n <th>%d</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n', ...
                   NumL1Epochs, ...
                   NumL2Epochs, ...
                   NumL5Epochs, ...
                   MeanL1Diff, ...
                   MeanL2Diff, ...
                   MeanL5Diff, ...
                   MeanL1Sigma, ...
                   MeanL2Sigma, ...
                   MeanL5Sigma);
  elseif(L2Present == TRUE)
    fprintf(fDiffHTML,'<th>Total</th><th>GPS</th>\n <th>%d</th>\n <th>%d</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n', ...
                   NumL1Epochs, ...
                   NumL2Epochs, ...
                   MeanL1Diff, ...
                   MeanL2Diff, ...
                   MeanL1Sigma, ...
                   MeanL2Sigma);
  elseif(L5Present == TRUE)
    fprintf(fDiffHTML,'<th>Total</th><th>GPS</th>\n <th>%d</th>\n <th>%d</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n', ...
                   NumL1Epochs, ...
                   NumL5Epochs, ...
                   MeanL1Diff, ...
                   MeanL5Diff, ...
                   MeanL1Sigma, ...
                   MeanL5Sigma);
  else
    fprintf(fDiffHTML,'<th>Total</th><th>GPS</th>\n <th>%d</th>\n <th>%.3f</th>\n <th>%.3f</th>\n', ...
                   NumL1Epochs, ...
                   MeanL1Diff, ...
                   MeanL1Sigma);
  end

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

  
  %
  % GLN
  %
  if(GLONASSInstalled)
    fprintf(fDiffHTML,'<br/><br/><center>');
    fprintf(fDiffHTML,'<table border="1">');
    fprintf(fDiffHTML,'<tr>\n');
  
    % Write out the table header
    fprintf(fDiffHTML,'<th>SV #</th>\n');
    fprintf(fDiffHTML,'<th>System</th>\n');
    fprintf(fDiffHTML,'<TH># L1 Epochs</TH>\n');
  
    if(GL2_P_present == TRUE)
      fprintf(fDiffHTML,'<TH># L2 Epochs</TH>\n');
    end
  
    if(GL1_P_present == TRUE)
      fprintf(fDiffHTML,'<TH># L1-P Epochs</TH>\n');
    end

    fprintf(fDiffHTML,'<TH>L1 Mean Difference [dB]</TH>\n');
    
    if(GL2_P_present == TRUE)
      fprintf(fDiffHTML,'<TH>L2 Mean Difference [dB]</TH>\n');
    end
    
    if(GL1_P_present == TRUE)
      fprintf(fDiffHTML,'<TH>L1-P Mean Difference [dB]</TH>\n');
    end
    
    fprintf(fDiffHTML,'<TH>L1 Sigma [dB]</TH>\n');
    
    if(GL2_P_present == TRUE)
      fprintf(fDiffHTML,'<TH>L2 Sigma [dB]</TH>\n');
    end
    
    if(GL1_P_present == TRUE)
      fprintf(fDiffHTML,'<TH>L1-P Sigma [dB]</TH>\n');
    end
    fprintf(fDiffHTML,'</tr>\n');

    % Now process GLONASS
    for SV=1:24
      if(isfinite(GLNCNRDiffL1CA(1,SV)))
        % Output the data to a space delimited text file

        fprintf(fDiffHTML,'<tr>\n');
        if( (GL2_P_present == TRUE) & (GL1_P_present == TRUE) )
          fprintf(fDiff,'%d 2 %d %d %d %.3f %.3f %.3f %.3f %.3f %.3f\n', ...
                     SV, ...
                     GLNCNRDiffL1CA(2,SV), ...
                     GLNCNRDiffL2(2,SV), ...
                     GLNCNRDiffL1P(2,SV), ...
                     GLNCNRDiffL1CA(1,SV), ...
                     GLNCNRDiffL2(1,SV), ...
                     GLNCNRDiffL1P(1,SV), ...
                     GLNCNRDiffL1CA(3,SV), ...
                     GLNCNRDiffL2(3,SV), ...
                     GLNCNRDiffL1P(3,SV));
          fprintf(fDiffHTML,'<td>%d</td><td>GLN</td>\n <td>%d</td>\n <td>%d</td>\n <td>%d</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n  <td>%.3f</td>\n', ...
                     SV, ...
                     GLNCNRDiffL1CA(2,SV), ...
                     GLNCNRDiffL2(2,SV), ...
                     GLNCNRDiffL1P(2,SV), ...
                     GLNCNRDiffL1CA(1,SV), ...
                     GLNCNRDiffL2(1,SV), ...
                     GLNCNRDiffL1P(1,SV), ...
                     GLNCNRDiffL1CA(3,SV), ...
                     GLNCNRDiffL2(3,SV), ...
                     GLNCNRDiffL1P(3,SV));
        elseif(GL2_P_present == TRUE)
          fprintf(fDiff,'%d 2 %d %d %.3f %.3f %.3f %.3f\n', ...
                     SV, ...
                     GLNCNRDiffL1CA(2,SV), ...
                     GLNCNRDiffL2(2,SV), ...
                     GLNCNRDiffL1CA(1,SV), ...
                     GLNCNRDiffL2(1,SV), ...
                     GLNCNRDiffL1CA(3,SV), ...
                     GLNCNRDiffL2(3,SV));
          fprintf(fDiffHTML,'<td>%d</td><td>GLN</td>\n <td>%d</td>\n <td>%d</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n  <td>%.3f</td>\n', ...
                     SV, ...
                     GLNCNRDiffL1CA(2,SV), ...
                     GLNCNRDiffL2(2,SV), ...
                     GLNCNRDiffL1CA(1,SV), ...
                     GLNCNRDiffL2(1,SV), ...
                     GLNCNRDiffL1CA(3,SV), ...
                     GLNCNRDiffL2(3,SV));
        elseif(GL1_P_present == TRUE)
          fprintf(fDiff,'%d 2 %d %d %.3f %.3f %.3f %.3f\n', ...
                     SV, ...
                     GLNCNRDiffL1CA(2,SV), ...
                     GLNCNRDiffL1P(2,SV), ...
                     GLNCNRDiffL1CA(1,SV), ...
                     GLNCNRDiffL1P(1,SV), ...
                     GLNCNRDiffL1CA(3,SV), ...
                     GLNCNRDiffL1P(3,SV));
          fprintf(fDiffHTML,'<td>%d</td><td>GLN</td>\n <td>%d</td>\n <td>%d</td>\n <td>%.3f</td>\n <td>%.3f</td>\n <td>%.3f</td>\n  <td>%.3f</td>\n', ...
                     SV, ...
                     GLNCNRDiffL1CA(2,SV), ...
                     GLNCNRDiffL1P(2,SV), ...
                     GLNCNRDiffL1CA(1,SV), ...
                     GLNCNRDiffL1P(1,SV), ...
                     GLNCNRDiffL1CA(3,SV), ...
                     GLNCNRDiffL1P(3,SV));
        else
          fprintf(fDiff,'%d 2 %d %.3f %.3f\n', ...
                     SV, ...
                     GLNCNRDiffL1CA(2,SV), ...
                     GLNCNRDiffL1CA(1,SV), ...
                     GLNCNRDiffL1CA(3,SV));
          fprintf(fDiffHTML,'<td>%d</td><td>GLN</td>\n <td>%d</td>\n <td>%.3f</td>\n <td>%.3f</td>\n', ...
                     SV, ...
                     GLNCNRDiffL1CA(2,SV), ...
                     GLNCNRDiffL1CA(1,SV), ...
                     GLNCNRDiffL1CA(3,SV));
        end
        fprintf(fDiffHTML,'</tr>\n');
      end
    end

    i = find(isfinite(GLNCNRDiffL1CA(2,:)) & (GLNCNRDiffL1CA(3,:) < CNRTestSigmaMask ) );
    NumL1CAEpochs = sum(GLNCNRDiffL1CA(2,i)); 
    MeanL1CADiff  = sum(GLNCNRDiffL1CA(2,i).*GLNCNRDiffL1CA(1,i)) / NumL1CAEpochs;
    MeanL1CASigma = mean(GLNCNRDiffL1CA(3,i));

    if(GL2_P_present == TRUE)
      i = find(isfinite(GLNCNRDiffL2(2,:)) & (GLNCNRDiffL2(3,:) < CNRTestSigmaMask ) );
      NumL2Epochs = sum(GLNCNRDiffL2(2,i));
      MeanL2Diff  = sum(GLNCNRDiffL2(2,i).*GLNCNRDiffL2(1,i)) / NumL2Epochs;
      MeanL2Sigma = mean(GLNCNRDiffL2(3,i));
    end

    if(GL1_P_present == TRUE)
      i = find(isfinite(GLNCNRDiffL1P(2,:)) & (GLNCNRDiffL1P(3,:) < CNRTestSigmaMask ) );
      NumL1PEpochs = sum(GLNCNRDiffL1P(2,i));
      MeanL1PDiff  = sum(GLNCNRDiffL1P(2,i).*GLNCNRDiffL1P(1,i)) / NumL1PEpochs;
      MeanL1PSigma = mean(GLNCNRDiffL1P(3,i));
    end

    fprintf(fDiffHTML,'</tr><tr>');
    if(GL2_P_present == TRUE)  & (GL1_P_present == TRUE)
      fprintf(fDiffHTML,'<th>Total</th><th>GLN</th>\n <th>%d</th>\n <th>%d</th>\n <th>%d</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n', ...
                     NumL1CAEpochs, ...
                     NumL2Epochs, ...
                     NumL1PEpochs, ...
                     MeanL1CADiff, ...
                     MeanL2Diff, ...
                     MeanL1PDiff, ...
                     MeanL1CASigma, ...
                     MeanL2Sigma, ...
                     MeanL1PSigma);
    elseif(GL2_P_present == TRUE)
      fprintf(fDiffHTML,'<th>Total</th><th>GLN</th>\n <th>%d</th>\n <th>%d</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n', ...
                     NumL1CAEpochs, ...
                     NumL2Epochs, ...
                     MeanL1CADiff, ...
                     MeanL2Diff, ...
                     MeanL1CASigma, ...
                     MeanL2Sigma);
    elseif(GL1_P_present == TRUE)
      fprintf(fDiffHTML,'<th>Total</th><th>GLN</th>\n <th>%d</th>\n <th>%d</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n <th>%.3f</th>\n', ...
                     NumL1CAEpochs, ...
                     NumL1PEpochs, ...
                     MeanL1CADiff, ...
                     MeanL1PDiff, ...
                     MeanL1CASigma, ...
                     MeanL1PSigma);
    else
      fprintf(fDiffHTML,'<th>Total</th><th>GLN</th>\n <th>%d</th>\n <th>%.3f</th>\n <th>%.3f</th>\n', ...
                     NumL1CAEpochs, ...
                     MeanL1CADiff, ...
                     MeanL1CASigma);
    end

    fprintf(fDiffHTML,'</table>\n');
    fprintf(fDiffHTML,'</center>\n');
  end % if(GLONASSInstalled)
  
  fprintf(fDiffHTML,'The averages above are weighted according to the number of epochs of data present for each satellite.<br>\n');
  fprintf(fDiffHTML,'Satellites with C/No difference Sigmas >= %.2fdB are ignored in the averages above.<BR>\n',CNRTestSigmaMask); 
  if(CNRTestElMask > 0)
    fprintf(fDiffHTML,'Data from epochs where the Elevation angle was less than %.2f Degrees were not included in the above statistics<br>\n',CNRTestElMask);
  end
  fprintf(fDiffHTML,'<hr>\n');
  fprintf(fDiffHTML, '%s', timeStr);
  fprintf(fDiffHTML,'</body>\n');
  fprintf(fDiffHTML,'</html>\n');

  fclose(fDiff);
  fclose(fDiffHTML);
end  % End if(SkipCNRTest == 0)


if(QuickMode < 2)
  %*******************************************************************
  % Compute the trends (DD-sigmas versus elevation and C/No)
  % and histograms of the DD observables.
  %*******************************************************************]

  %
  % First process the GPS data
  %
  
  if (L2Present)
    pddL2 = phaseDDiffL2(:,1:33);
    rddL2 = rangeDDiffL2(:,1:33);
  else
    pddL2 = NaN;
    rddL2 = NaN;
  end
  
  if (L5Present)
    pddL5 = phaseDDiffL5(:,1:33);
    rddL5 = rangeDDiffL5(:,1:33);
  else
    pddL5 = NaN;
    rddL5 = NaN;
  end
    
  PlotTrend( phaseDDiffL1(:,1:33), ...
             pddL2, ...
             pddL5, ...
             rangeDDiffL1(:,1:33), ...
             rddL2, ...
             rddL5, ...
             Elev, ...
             CNRL1, ...
             CNRL2, ...
             CNRL5, ...
             timeStr, ...
             PlotPNG, ...
             CarrMax, ...
             CodeMax, ...
             L2Present, ...
             L5Present, ...
             0);
  
  clear pddL2 pddL5 rddL2 rddL5;
           
  %
  % Now process the GLONASS data
  %
  if(GLONASSInstalled)
    
    if (GL2_P_present)
      pddGL2 = phaseDDiffGL2(:,1:25);
      rddGL2 = rangeDDiffGL2(:,1:25);
    else
      pddGL2 = NaN;
      rddGL2 = NaN;
    end

    if (GL1_P_present)
      pddGL1_P = phaseDDiffGL1_P(:,1:25);
      rddGL1_P = rangeDDiffGL1_P(:,1:25);
    else
      pddGL1_P = NaN;
      rddGL1_P = NaN;
    end
        
    PlotTrend(phaseDDiffGL1_CA(:,1:25), ...
              pddGL2, ...
              pddGL1_P, ...
              rangeDDiffGL1_CA(:,1:25), ...
              rddGL2, ...
              rddGL1_P, ...
              NaN, ...
              NaN, ...
              NaN, ...
              NaN, ...
              timeStr, ...
              PlotPNG, ...
              CarrMax/5, ...
              CodeMax, ...
              L2Present, ...
              GL1_P_present, ...
              2);
            
    clear pddGL2 rddGL2 pddGL1_P rddGL1_P;
            
  end % if(GLONASSInstalled)
  
end
% Clean-up
close all

% No longer required.
clear rangeDDiffL1;
clear rangeDDiffL2;
clear Azi;
clear CNRL1;
clear CNRL2;

if(QuickMode < 2)
  ElevL1      = [phaseDDiffL1(:,1) Elev]; % pre-pend the time
  clear Elev;
  
  if(GLONASSInstalled)
    ElevGLNL1CA = [phaseDDiffGL1_CA(:,1) ElevL1CA]; % pre-pend the time
    clear ElevL1CA;
  end

  if(L2Present == TRUE)
    ElevL2 = [phaseDDiffL2(:,1) ElevL2];
  else
    ElevL2 = [];
  end
  
  if(GLONASSInstalled)
    if(GL2_P_present == TRUE)
      ElevGLNL2 = [phaseDDiffGL2(:,1) ElevL2G];
    else
      ElevGLNL2 = [];
    end
  end

  clear phaseDDiffL2;
  if(GLONASSInstalled)
    clear phaseDDiffGL2;
  end

  %*******************************************************************
  % Create a plot of the satellites tracked versus time.
  %*******************************************************************
  plotSVTrack(ElevL1,ElevL2,phaseDDiffL1,PlotPNG,L2Present,timeStr,0);
  clear ElevL1;
  
  if(GLONASSInstalled && ~isempty(phaseDDiffGL1_CA))
    plotSVTrack(ElevGLNL1CA,ElevGLNL2,phaseDDiffGL1_CA,PlotPNG,L2Present,timeStr,2);
    clear ElevGLNL1CA;
  end
else  
  clear Elev;
  clear phaseDDiffL2;
end

% Clean up the unused variables
clear ElevL2;
clear phaseDDiffL1;

%*********************************************************************
% Perform the ionospheric processing.
%*********************************************************************
if( (QuickMode == 0)  & (L2Present == TRUE))
  close all
  plotIono(ionoDDiff,PlotPNG,timeStr);
end

if(QuickMode < 2)
  fprintf(fid, '%s', timeStr);
  fprintf(fid,'<br>\n');
    
  % Save the Version Numbers into the HTML so we know which version created it.
  if strcmp(piProVer.buildType,'Test')
    fprintf(fid,'Created using a %s build of PiPro derived from %s, which uses %s.\n', ...
      piProVer.buildType, piProVer.piproVersion, piProVer.pilineVersion);
  else
    fprintf(fid,'Created using a %s build of PiPro with CVS Tag %s, which uses %s.\n', ...
      piProVer.buildType, piProVer.piproVersion, piProVer.pilineVersion);
  end

  % Terminate all the HTML tags and close the document
  fprintf(fid,'</font></BODY></HTML>');
  fclose(fid);
end

% 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.
if( Compiled == TRUE )
  close all
end

fclose('all');

disp('Goodbye...');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% $Log: load_pi.m,v $
% Revision 1.64  2008/02/25 16:35:30  stuart
% o Update FDMA numbers
%
% Revision 1.63  2008/01/04 20:50:32  vnadkar
% Various upgrades and fixes to piPro.
%
% M load_pi.m
% o Support 32 GPS satellites.
% o Allow user to save output of PiLineMex.dll into a MAT-file.
% o Allow user to load output of MAT-file instead of re-running PiLineMex.dll.
%
% Revision 1.62  2008/01/04 01:13:06  stuart
% o Update to the latest known FDMA channel numbers
%
% Revision 1.61  2007/10/10 20:23:35  vnadkar
% Add support for Galileo Data spoofed as GPS data.
%
% M load_pi.m
% o Specifically, support L5 when L2 is absent.
%
% Revision 1.60  2007/10/09 20:40:52  vnadkar
% Add support for Chogori data
%
% Chogori data does not contain GLONASS L1P data. Previously piPro would
% crash when these data were missing. This fix allows for piPro to run when
% GLONASS L1P data is absent from the data files
%
% M piPro/load_pi.m
% o Allow for missing GLONASS L1P data.
% o The L2Present flag was being used for both GPS and GLONASS
% data. Now create separate flags GL2_P_present and GL2_CA_present
% for GLONASS L2 data. Note that GL2_CA is not currently supported
% in piPro. Replace L2Present with GL2_P_present in the relevant
% locations.
%
% M piPro/loadPiHeaderHTM.m
% o Use the GL2_P_present flag for GLONASS instead of sharing
% the L2Present flag with GPS.
%
% Revision 1.59  2007/03/13 01:02:25  vnadkar
% o Print the PiPro version number before any processing is done. This
% allows the user to know what the version of PiPro is even if processing
% fails.
%
% Revision 1.58  2007/02/01 00:56:40  vnadkar
% M load_pi.m
% o Update canned frequency channel mappings.
% o Add comments, specifying where to get the canned mappings from.
%
% Revision 1.57  2007/02/01 00:43:20  vnadkar
% M load_pi.m
% o Keep figures hidden when plotting residuals for GLONASS.
%
% M suptitle.m
% o Don't make a figure visible when adding a top level title using suptitle.
%
% Revision 1.56  2007/01/29 19:36:53  vnadkar
% load_pi.m
% makePiPro.m
% pipro_version.m
% o Support version numbers for PiPro.
%
% README_pipro.txt
% o Installation instructions for PiPro executable (load_pi.exe).
%
% Revision 1.55  2006/04/27 00:22:03  vnadkar
% o Updates to allow load_pi to run with missing data types (e.g. missing GLONASS etc).
%
% Revision 1.54  2005/09/20 20:45:10  acartmell
% o Fix for GPS only processing.
%   Only align the GLONASS data arrays if there is some GLONASS data present.
%
% Revision 1.53  2005/08/17 02:10:46  acartmell
% o Fix the GPS all satellite, user scaled phase double difference plot error (I think).
%   May have been Windows changing to Current plot status as it worked in debug.
%   Better closing of plots seems to fix it so far.
%
% o Fix for data sets where the first SV has L2 only data.
%   Moved the increment of the colour index to outside the L1 data loop.
%
% Revision 1.52  2005/07/19 13:02:28  timo
% Fix errors in fprintf() when strings contain file names containing
% paths (which use \ which is a special char). Must use %s for
% the string and not use the string directly.
%
% Revision 1.51  2005/06/23 03:26:16  acartmell
% o Update for detection of outlier data points
%   - Generates two versions of the 'all satellites' plots.
%     One applies the user defined Y-axis limits, the other autoscales the Y-axis.
%     The X-axis is limited in both cases by the data time span, as normal.
%     Added new link on the web page to access the autoscaled plots.
%   - For the per satellite DDiff plots (both code and carrier):
%       if the max(abs(DDiff)) > (2*User Supplied Axis Limit) then autoscale the Y-axis to view outliers
%       otherwise apply the user supplied limits
%     The autoscaling selection is done separately for each frequency.
%     The X-axis is limited by the data time span, as normal.
%     The autoscaled plot will replace the standard plot.
%
% Revision 1.50  2005/06/01 18:46:28  acartmell
% o Fix to handle PiLine returning arrays of different sizes (and non-aligned timtags) for L1-C/A, L1-P and L2 data.
% o Generate temporary arrays for elevation storage to avoid crashes when running with CNR plots off.
%
% Revision 1.49  2005/05/26 22:20:36  acartmell
% o Update to produce, save and add web page link to plots that display the DD phase DD range for all GLONASS satellites.
%
% Revision 1.48  2005/05/25 15:20:26  stuart
% o Revert load_pi.m back to a function
% o Add Timo's latest PiLine DLL
%
% Revision 1.47  2005/05/18 14:07:15  stuart
% o Update phase processing
% - Add GLONASS support
%
% Revision 1.46  2005/04/18 22:02:35  stuart
% o Add two satellite system support
%
% Revision 1.45  2005/04/18 19:46:46  vnadkar
% Make plots invisible. This allows the user to work on other things while
% load_pi runs in the background. The outputs from the plots are printed
% to files.
%
% Revision 1.44  2004/08/26 16:24:48  stuart
% o Correct commandline parsing. If the decimate interval was added it 
% forced you to also add the start and stop times, removed this dependency
%
% Revision 1.43  2004/03/16 12:21:46  timo
% Fix colorbar and polar axes handle assignments.
%
% Revision 1.42  2004/01/30 16:23:45  stuart
% o Fix for an axes handling problem from Timo.
%
% Revision 1.41  2004/01/28 01:51:56  stuart
% o Improve L5 processing by adding correlation calculation,
%
% Revision 1.38  2004/01/10 00:05:13  stuart
% o Update for L5
%
% Revision 1.37  2003/10/02 02:00:59  stuart
% o Solve a problem when the number of rows (epochs of data) was less
% than the number of columns. The length command was being used which 
% always returns the longest dimension.
%
% Revision 1.36  2003/09/30 17:12:57  stuart
% o Corrected command line parsing.
%
% Revision 1.34  2003/01/07 21:53:08  stuart
% o Add support for decimation to piPro
%
% Revision 1.33  2002/06/05 00:48:00  stuart
% o Added the ability to not print a legend on some of the C/No plots. 
%   This was a request made from the RF group, with lots of satellites 
%   the legend can sometimes obscure important events.
%
% Revision 1.32  2002/05/31 00:44:51  stuart
% o Corrected a bug in the ACF code
% o Added a new mode which allows the C/No tests to be turned off.
%
% Revision 1.31  2002/02/21 17:23:39  stuart
% o Updates to use Timo's latest labelaxes function.
%
% Revision 1.30  2002/02/13 01:49:18  stuart
% o Solve a bug where the case of a variable was wrong.
%
% Revision 1.29  2002/01/05 01:42:03  stuart
% o Added an elevation and C/No sigma filter threshold for use in the 
%   C/No analysis for Pogo testing
%
% Revision 1.28  2001/12/14 15:47:25  stuart
% o Updates for use in Pogo testing
% o Fixed a few minor bugs
% o Some restructuring to make the code more maintainable.
%
% Revision 1.27  2001/11/27 19:29:48  stuart
% o Updated to account for the new order in PiLineMex.dll of the base 
%   and rover C/No data.
%
% Revision 1.26  2001/11/15 23:57:52  stuart
% o Add processing of the Rover C/No
%
% Revision 1.25  2001/11/06 23:56:07  stuart
% o Clean up on some of the labels on varous graphs
% o Added a "QuickMode" so only a sub-set of the graphs get created
%
% Revision 1.24  2001/10/26 17:45:35  stuart
% o Timo corrected a problem in runPiLine which meant that a 
%   temporary fix in load_pi could be removed.
%
% Revision 1.23  2001/10/25 21:58:16  stuart
% o Change the logo to a .png file (its smaller)
% o Remove some debug code that prevented input options from being 
%   accepted.
%
% Revision 1.22  2001/10/25 21:36:46  stuart
% o Bug fixes for when there is no L2 data
% o Updated Logo data
%
% Revision 1.21  2001/09/19 18:19:26  stuart
% o Resolved a couple of minor bugs
% o Created a seperate file for the DD trend & distribution graphs and 
%   analysis
% o Added DD versus C/No trend analysis
% o Added a second string to suptitle
% o Matlab creates a .gif logo that is referenced in the HTML results
%
% Revision 1.20  2001/09/07 19:39:49  stuart
% o Resolved a bug with the input option
% o Minor code tidy up
%
% Revision 1.19  2001/08/29 15:07:53  stuart
% o Solved a problem in the compiled version
% o Added text to the graphs that gives the length of the session and 
%   files used
% o Improved the HTML output
%
% Revision 1.18  2001/08/20 19:26:23  stuart
% o Updates to add observable histograms
%
% Revision 1.17  2001/06/28 19:48:16  stuart
% o Call MEX/DLL piline directly in load_pi
% o Produce a CSV file of stats in pos
% o Update Perl scripts to generate the 'pos' webpages.
%
% Revision 1.16  2001/05/07 16:55:38  stuart
% o Add:
%   Satellite Tracking Plot
%   Ionospheric Double Difference
%   L1 CNR against Elevation
%   L2 CNR against Elevation.
%
% Revision 1.15  2001/05/03 22:28:56  stuart
% o Improved CNR versus elevation/azimuth plots with code provided by 
%   Timo.
%
% Revision 1.14  2001/05/01 20:00:53  stuart
% o Add extra comments
% o Solve a bug where if there was a single reference SV for the whole 
%   file the elevation/azimuth data for that satellite would not be 
%   plotted
% o Add text status messages to the command window so the user knows 
%   what is happening (ie so that on a large data set we know it 
%   hasn't crashed).
%
% Revision 1.13  2001/05/01 15:34:54  stuart
% o By not removing the NaN can simulate a pen-up/pen-down, ie when 
%   there is missing data the plot does not draw a line between the two 
%   surrounding points.
% o Add first attempt at providing a polar plot (el/az) with colour 
%   representing the CNR in dBHz.
%
% Revision 1.12  2001/04/27 00:24:01  stuart
% o Updates to support elevation/azimuth plots
%
% Revision 1.11  2001/04/25 15:58:50  stuart
% o Provide a command line on/off for the FFT plots
% o Add a polar plot of elevation/azimuth.
%
% Revision 1.10  2001/04/23 15:33:38  stuart
% o Bug fixes and new scripts to plot elevation/azimuth.
%
% Revision 1.9  2001/04/03 22:10:06  stuart
% o Updates to get the .m files to compile and take advantage of the 
% fact that the latest Matlab compiler will allow .png files to be 
% written directly.
%
% Revision 1.7  2001/04/02 22:04:51  stuart
% o Updated with latest QA scripts
%
% Revision 1.6  2000/11/28 11:58:43  timo
% Update to allow piLine residuals file to be input directly (when
% piLine option used which does not output elevation/CNR).
% Also ignores first 12 lines of comments in piLine files.
%
% Revision 1.6  2000/11/23 16:34:29  timo
% Updates to enable direct reading of piLine residuals files (when
% piLine options set to only output residuals without elevations and
% CNRs). Also ignores first 12 lines of piLine files (comment section
% prior to residuals).
%
% Revision 1.5  2000/09/24 19:48:57  timo
% Remove highlight markers, >>>>> (changes between this and
% last rev). No functional changes.
%
% Revision 1.4  2000/09/22 18:49:36  timo
% Update to allow Matlab compilation (StandAlone application)
%
% Revision 1.3  2000/08/21 13:57:33  stuart
% o Correct $Id
% o Add $Log
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
