% function plot_obs(Compiled,CodeMax,CarrMax,Title,PNGPlot);

Compiled = 0;
CodeMax  = 2;
CarrMax  = 250;
Title    = 'LSQA';
PNGPlot  = 0;
nargin   = 5;

%
% $Id: plot_obs.m,v 1.7 2005/04/18 19:46:46 vnadkar Exp $
%
% Script to analyze the statistical observable data
% provided by the automatic test system
%
% To use this script....
%
% Place the following files from a TRS test system in the 
% current directory (standalone run PiPro on a pair of Survey Data 
% files):
%
%    stats_l1.ph
%    stats_l1.rg
%    stats_l2.ph
%    stats_l2.rg
%    
%    el_l1.ph
%    el_l1.rg
%    el_l2.ph
%    el_l2.rg
%
% plot_obs(Compiled,CodeMax,CarrMax,Title,PNGPlot);
%   
%   Compiled - (=0) Called by Matlab; (=1) StandAlone application.
%   CodeMax  - Y-scaling for Code [Metres]
%   CarrMax  - Y-scaling for Carr [Metres]
%   Title    - Supertitle place on each graph
%   PlotPNG  - (=0) Don't Create .png files; (=1) Create .png files

% 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.
%
% 5/7/00 - Stuart Riley X8696
%
% 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.

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

% Check for input parameters, and convert from strings to numbers.
% plot_obs(Compiled,CodeMax,CarrMax,Title,PNGPlot);

% 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
end

if( nargin == 0 )
  Compiled  = FALSE ;
  Title     = '';
  AutoScale = TRUE;
  PNGPlot   = TRUE;
elseif(nargin <= 2 )
  Title     = '';
  AutoScale = TRUE;
  PNGPlot   = TRUE;
elseif(nargin == 3)
  Title     = '';
  if( Compiled == TRUE )
    CodeMax = str2num( CodeMax ) ;
    CarrMax = str2num( CarrMax ) ;
  % else leave CodeMax and CarrMax as passed (numbers).
  end
  AutoScale = FALSE;
  PNGPlot   = TRUE;
%!!!TA The above line was missing.
elseif(nargin == 4)
  if( Compiled == TRUE )
    CodeMax = str2num( CodeMax ) ;
    CarrMax = str2num( CarrMax ) ;
  % else leave CodeMax and CarrMax as passed (numbers).
  end
  AutoScale = FALSE;
  PNGPlot   = TRUE;
else
  if( Compiled == TRUE )
    CodeMax = str2num( CodeMax ) ;
    CarrMax = str2num( CarrMax ) ;
    PNGPlot = str2num( PNGPlot) ;
  % else leave CodeMax, CarrMax, and PNGPlot as passed (numbers).
  end
  AutoScale = FALSE;
end

if( Compiled == FALSE )
  % Load the stats which are one hour segments
  % grouped into 5degree elevation bins.
  load el_l1.ph
  elev_phase_l1 = el_l1;
  load el_l1.rg
  elev_range_l1 = el_l1;
  load el_l2.ph
  elev_phase_l2 = el_l2;
  load el_l2.rg
  elev_range_l2 = el_l2;

  % Load the stats which are one hour segments
  % grouped into 3dB CNR bins.
  load cnr_l1.ph
  cnr_phase_l1 = cnr_l1;
  load cnr_l1.rg
  cnr_range_l1 = cnr_l1;
  load cnr_l2.ph
  cnr_phase_l2 = cnr_l2;
  load cnr_l2.rg
  cnr_range_l2 = cnr_l2;

  % Load the stats which are one hour segments
  % analyzed on a per satellite basis.

  load stats_l1.ph
  stats_phase_l1 = stats_l1;
  load stats_l2.ph
  stats_phase_l2 = stats_l2;
  load stats_l1.rg
  stats_range_l1 = stats_l1;
  load stats_l2.rg
  stats_range_l2 = stats_l2;
else
  % Load the stats which are one hour segments
  % grouped into 5degree elevation bins.

  % Using "load" does not compile with the Matlab compiler, and is
  % replaced by the alternative technique below for all data files.

  fid = fopen( 'el_l1.ph', 'rt' ) ;
  el_l1 = fscanf( fid, '%f', [6, inf] ) ;  % Start with a vector.
  el_l1 = el_l1' ;  % Transpose as fscanf() fills by column.
  fclose( fid ) ;
  elev_phase_l1 = el_l1;

  % Load remaining files using same technique as above.
  fid = fopen( 'el_l1.rg', 'rt' ) ;
  el_l1 = fscanf( fid, '%f', [6, inf] ) ;
  el_l1 = el_l1' ;
  fclose( fid ) ;
  elev_range_l1 = el_l1;

  fid = fopen( 'el_l2.ph', 'rt' ) ;
  el_l2 = fscanf( fid, '%f', [6, inf] ) ;
  el_l2 = el_l2' ;
  fclose( fid ) ;
  elev_phase_l2 = el_l2;

  fid = fopen( 'el_l2.rg', 'rt' ) ;
  el_l2 = fscanf( fid, '%f', [6, inf] ) ;
  el_l2 = el_l2' ;
  fclose( fid )
  elev_range_l2 = el_l2;

  % Load the stats which are one hour segments
  % grouped into 3dB CNR bins.

  % See above for same technique

  fid = fopen( 'cnr_l1.ph', 'rt' ) ;
  cnr_l1 = fscanf( fid, '%f', [6, inf] ) ;
  cnr_l1 = cnr_l1' ;
  fclose( fid )
  cnr_phase_l1 = cnr_l1;

  fid = fopen( 'cnr_l1.rg', 'rt' ) ;
  cnr_l1 = fscanf( fid, '%f', [6, inf] ) ;
  cnr_l1 = cnr_l1' ;
  fclose( fid )
  cnr_range_l1 = cnr_l1;

  fid = fopen( 'cnr_l2.ph', 'rt' ) ;
  cnr_l2 = fscanf( fid, '%f', [6, inf] ) ;
  cnr_l2 = cnr_l2' ;
  fclose( fid )
  cnr_phase_l2 = cnr_l2;

  fid = fopen( 'cnr_l2.rg', 'rt' ) ;
  cnr_l2 = fscanf( fid, '%f', [6, inf] ) ;
  cnr_l2 = cnr_l2' ;
  fclose( fid )
  cnr_range_l2 = cnr_l2;

  % Load the stats which are one hour segments
  % analyzed on a per satellite basis.

  fid = fopen( 'stats_l1.ph', 'rt' ) ;
  stats_l1 = fscanf( fid, '%f', [7, inf] ) ;
  stats_l1 = stats_l1' ;
  fclose( fid )
  stats_phase_l1 = stats_l1;

  fid = fopen( 'stats_l1.rg', 'rt' ) ;
  stats_l1 = fscanf( fid, '%f', [7, inf] ) ;
  stats_l1 = stats_l1' ;
  fclose( fid )
  stats_range_l1 = stats_l1;

  fid = fopen( 'stats_l2.ph', 'rt' ) ;
  stats_l2 = fscanf( fid, '%f', [7, inf] ) ;
  stats_l2 = stats_l2' ;
  fclose( fid )
  stats_phase_l2 = stats_l2;

  fid = fopen( 'stats_l2.rg', 'rt' ) ;
  stats_l2 = fscanf( fid, '%f', [7, inf] ) ;
  stats_l2 = stats_l2' ;
  fclose( fid )
  stats_range_l2 = stats_l2;
end


% Number of epochs in each statistic which must be
% exceeded before the data is plotted/analyzed further.
ACCEPTABLE_SAMPLE_SET = 50;
fig(10) = figure('Visible','off');
fig(1) = figure('Visible','off');

% L1 phase Sigma
subplot(2,2,1)
j = find(elev_phase_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(elev_phase_l1(j,2),elev_phase_l1(j,5),'x');
grid
xlabel('Elev. Angle [Deg]')
ylabel('Sigma [mcycle]')
title('L1 Phase DD Res');

max_l1 = max(elev_phase_l1(j,5));

% Form the weighed average
weighted_phase_l1 = weight(0:5:90, elev_phase_l1, j, 5);
hold
plot(0:5:90,weighted_phase_l1(:,2),'r')

% L2 phase Sigma
subplot(2,2,2)
j = find(elev_phase_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(elev_phase_l2(j,2),elev_phase_l2(j,5),'x');
grid
xlabel('Elev. Angle [Deg]')
ylabel('Sigma [mcycle]')
title('L2 Phase DD Res');

max_l2 = max(elev_phase_l2(j,5));

% Form the weighed average
weighted_phase_l2 = weight(0:5:90, elev_phase_l2, j, 5);
hold
plot(0:5:90,weighted_phase_l2(:,2),'r')

% figure(fig(10))
set(0,'CurrentFigure',fig(10));
subplot(2,1,1);

H = bar(0:5:90, [weighted_phase_l1(:,2), weighted_phase_l2(:,2)],'grouped');

% L1 Data
XData = get(H(1),'XData');
YData = get(H(1),'YData');
X = mean([XData(2,:)' , XData(3,:)']');
for i=1:19
  kk = find(elev_phase_l1(:,2) == (i-1)*5);
  if(kk > 0)
    TextArray(i) = {num2str(sum(elev_phase_l1(kk,3)))};
  else
    TextArray(i) = {''};
  end
end
text(X,YData(2,:)+2,TextArray,'Rotation',90,'FontSize',7);

% L2 Data
XData = get(H(2),'XData');
YData = get(H(2),'YData');
X = mean([XData(2,:)' , XData(3,:)']');
for i=1:19
  kk = find(elev_phase_l2(:,2) == (i-1)*5);
  if(kk > 0)
    TextArray(i) = {num2str(sum(elev_phase_l2(kk,3)))};
  else
    TextArray(i) = {''};
  end
end
text(X,YData(2,:)+2,TextArray,'Rotation',90,'FontSize',7);

legend('L1','L2');
xlabel('Elevation [Degrees]');
ylabel('Sigma [mcycle]');
title('Phase DD');
grid on

% figure(fig(1))
set(0,'CurrentFigure',fig(1));


% Now make sure L1 & L2 sigma graphs have the same scale.

m = max([max_l1,max_l2]);
m = ceil(m/50)*50;
if(AutoScale == TRUE)
  % Need to make sure we can create a valid axis
  if(m ~= 0)
    axis([0 100 0 m]);  % L2
    subplot(2,2,1);
    axis([0 100 0 m]);  % L1
  end
else
  axis([0 100 0 CarrMax]);
  subplot(2,2,1);
  axis([0 100 0 CarrMax]);
end

% L1 phase mean
subplot(2,2,3)
j = find(elev_phase_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(elev_phase_l1(j,2),elev_phase_l1(j,4),'x');
grid
xlabel('Elev. Angle [Deg]')
ylabel('Mean [mcycle]')
title('L1 Phase DD Res');

% Form the weighed average
weighted_phase_l1 = weight(0:5:90, elev_phase_l1, j, 4);
hold
plot(0:5:90,weighted_phase_l1(:,2),'r')

max_l1 = max(elev_phase_l1(j,4));
min_l1 = min(elev_phase_l1(j,4));

% L2 phase mean
subplot(2,2,4)
j = find(elev_phase_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(elev_phase_l2(j,2),elev_phase_l2(j,4),'x');
grid
xlabel('Elev. Angle [Deg]')
ylabel('Mean [mcycle]')
title('L2 Phase DD Res');

% Form the weighed average
weighted_phase_l2 = weight(0:5:90, elev_phase_l2, j, 4);
hold
plot(0:5:90,weighted_phase_l2(:,2),'r')

max_l2 = max(elev_phase_l2(j,4));
min_l2 = min(elev_phase_l2(j,4));

max_mean = max([max_l1,max_l2]);
min_mean = min([min_l1,min_l2]);

max_mean = ceil(max_mean/50)*50;
min_mean = floor(min_mean/50)*50;

if(AutoScale == TRUE)
  % Need to make sure we can create a valid axis
  if(min_mean ~= max_mean)
    axis([0 100 min_mean max_mean]);  % L2
    subplot(2,2,3);
    axis([0 100 min_mean max_mean]);  % L1
  end
else
  axis([0 100 -CarrMax CarrMax]);
  subplot(2,2,3);
  axis([0 100 -CarrMax CarrMax]);
end
suptitle(Title);


% figure(fig(10))
set(0,'CurrentFigure',fig(10));
subplot(2,1,2);

H = bar(0:5:90, [weighted_phase_l1(:,2), weighted_phase_l2(:,2)],'grouped');

% L1 Data
XData = get(H(1),'XData');
YData = get(H(1),'YData');
kk = find(YData < 0);
YData(kk) = 0;
X = mean([XData(2,:)' , XData(3,:)']');
for i=1:19
  kk = find(elev_phase_l1(:,2) == (i-1)*5);
  if(kk > 0)
    TextArray(i) = {num2str(sum(elev_phase_l1(kk,3)))};
  else
    TextArray(i) = {''};
  end
end
text(X,YData(2,:)+2,TextArray,'Rotation',90,'FontSize',7);

% L2 Data
XData = get(H(2),'XData');
YData = get(H(2),'YData');
kk = find(YData < 0);
YData(kk) = 0;
X = mean([XData(2,:)' , XData(3,:)']');
for i=1:19
  kk = find(elev_phase_l2(:,2) == (i-1)*5);
  if(kk > 0)
    TextArray(i) = {num2str(sum(elev_phase_l2(kk,3)))};
  else
    TextArray(i) = {''};
  end
end
text(X,YData(2,:)+2,TextArray,'Rotation',90,'FontSize',7);

legend('L1','L2');
xlabel('Elevation [Degrees]');
ylabel('Mean [mcycle]');
title('Phase DD');
grid on

fig(2) = figure('Visible','off');
% L1 range Sigma
subplot(2,2,1)
j = find(elev_range_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(elev_range_l1(j,2),elev_range_l1(j,5)./100,'x');
grid
xlabel('Elev. Angle [Deg]')
ylabel('Sigma [M]')
title('L1 pseudorange DD Res');

max_l1 = max(elev_range_l1(j,5)./100);

% Form the weighed average
weighted_range_l1 = weight(0:5:90, elev_range_l1, j, 5);
hold
plot(0:5:90,weighted_range_l1(:,2)./100,'r')

% L2 range Sigma
subplot(2,2,2)
j = find(elev_range_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(elev_range_l2(j,2),elev_range_l2(j,5)./100,'x');
grid
xlabel('Elev. Angle [Deg]')
ylabel('Sigma [M]')
title('L2 pseudorange DD Res');

max_l2 = max(elev_range_l2(j,5)./100);

% Form the weighed average
weighted_range_l2 = weight(0:5:90, elev_range_l2, j, 5);
hold
plot(0:5:90,weighted_range_l2(:,2)./100,'r')

% Now make sure L1 & L2 sigma graphs have the same scale.

m = max([max_l1,max_l2]);
m = ceil(m);
if(AutoScale == TRUE)
  % Need to make sure we can create a valid axis
  if(m ~= 0)
    axis([0 100 0 m]);  % L2
    subplot(2,2,1);
    axis([0 100 0 m]);  % L1
  end
else
  axis([0 100 0 CodeMax]);
  subplot(2,2,1);
  axis([0 100 0 CodeMax]);
end

% L1 range mean
subplot(2,2,3)
j = find(elev_range_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(elev_range_l1(j,2),elev_range_l1(j,4)./100,'x');
grid
xlabel('Elev. Angle [Deg]')
ylabel('Mean [M]')
title('L1 pseudorange DD Res');

% Form the weighed average
weighted_range_l1 = weight(0:5:90, elev_range_l1, j, 4);
hold
plot(0:5:90,weighted_range_l1(:,2)./100,'r')

max_l1 = max(elev_range_l1(j,4)./100);
min_l1 = min(elev_range_l1(j,4)./100);

% L2 range mean
subplot(2,2,4)
j = find(elev_range_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(elev_range_l2(j,2),elev_range_l2(j,4)./100,'x');
grid
xlabel('Elev. Angle [Deg]')
ylabel('Mean [M]')
title('L2 pseudorange DD Res');

% Form the weighed average
weighted_range_l2 = weight(0:5:90, elev_range_l2, j, 4);
hold
plot(0:5:90,weighted_range_l2(:,2)./100,'r')

max_l2 = max(elev_range_l2(j,4)./100);
min_l2 = min(elev_range_l2(j,4)./100);

max_mean = max([max_l1,max_l2]);
min_mean = min([min_l1,min_l2]);

%
% Scale to closest 0.5m
%
max_mean = ceil(max_mean*2)/2; 
min_mean = floor(min_mean*2)/2;

axis([0 100 min_mean max_mean]);  % L1
if(AutoScale == TRUE)
  % Need to make sure we can create a valid axis
  if(min_mean ~= max_mean)
    axis([0 100 min_mean max_mean]);  % L2
    subplot(2,2,3);
    axis([0 100 min_mean max_mean]);  % L1
  end
else
  axis([0 100 -CodeMax CodeMax]);
  subplot(2,2,3);
  axis([0 100 -CodeMax CodeMax]);
end

suptitle(Title);

%
% Now plot the statistics against CNR
%

fig(3) = figure('Visible','off');

% L1 phase Sigma
subplot(2,2,1)
j = find(cnr_phase_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_phase_l1(j,2),cnr_phase_l1(j,5),'x');
grid
xlabel('CNR [dBHz]')
ylabel('Sigma [mcycle]')
title('L1 Phase DD Res');

max_l1 = max(cnr_phase_l1(j,5));

% Form the weighed average
weighted_phase_l1 = weight(30:3:60, cnr_phase_l1, j, 5);
hold
plot(30:3:60,weighted_phase_l1(:,2),'r')

% L2 phase Sigma
subplot(2,2,2)
j = find(cnr_phase_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_phase_l2(j,2),cnr_phase_l2(j,5),'x');
grid
xlabel('CNR [dBHz]')
ylabel('Sigma [mcycle]')
title('L2 Phase DD Res');

max_l2 = max(cnr_phase_l2(j,5));

% Form the weighed average
weighted_phase_l2 = weight(0:3:60, cnr_phase_l2, j, 5);
hold
plot(0:3:60,weighted_phase_l2(:,2),'r')

% Now make sure L1 & L2 sigma graphs have the same scale.

m = max([max_l1,max_l2]);
m = ceil(m/50)*50;

if(AutoScale == TRUE)
  % Need to make sure we can create a valid axis
  if(m ~= 0)
    axis([10 50 0 m]);  % L2
    subplot(2,2,1);
    axis([30 60 0 m]);  % L1
  end
else
  axis([10 50 0 CarrMax]);
  subplot(2,2,1);
  axis([30 60 0 CarrMax]);
end


% L1 phase mean
subplot(2,2,3)
j = find(cnr_phase_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_phase_l1(j,2),cnr_phase_l1(j,4),'x');
grid
xlabel('CNR [dBHz]');
ylabel('Mean [mcycle]')
title('L1 Phase DD Res');

% Form the weighed average
weighted_phase_l1 = weight(30:3:60, cnr_phase_l1, j, 4);
hold
plot(30:3:60,weighted_phase_l1(:,2),'r')

max_l1 = max(cnr_phase_l1(j,4));
min_l1 = min(cnr_phase_l1(j,4));

% L2 phase mean
subplot(2,2,4)
j = find(cnr_phase_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_phase_l2(j,2),cnr_phase_l2(j,4),'x');
grid
xlabel('CNR [dBHz]');
ylabel('Mean [mcycle]')
title('L2 Phase DD Res');

% Form the weighed average
weighted_phase_l2 = weight(0:3:60, cnr_phase_l2, j, 4);
hold
plot(0:3:60,weighted_phase_l2(:,2),'r')

max_l2 = max(cnr_phase_l2(j,4));
min_l2 = min(cnr_phase_l2(j,4));

max_mean = max([max_l1,max_l2]);
min_mean = min([min_l1,min_l2]);

max_mean = ceil(max_mean/50)*50;
min_mean = floor(min_mean/50)*50;

if(AutoScale == TRUE)
  % Need to make sure we can create a valid axis
  if(min_mean ~= max_mean)
    axis([10 50 min_mean max_mean]);  % L2
    subplot(2,2,3);
    axis([30 60 min_mean max_mean]);  % L1
  end
else
  axis([10 50 -CarrMax CarrMax]);
  subplot(2,2,3);
  axis([30 60 -CarrMax CarrMax]);
end

suptitle(Title);

fig(4) = figure('Visible','off');
% L1 range Sigma
subplot(2,2,1)
j = find(cnr_range_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_range_l1(j,2),cnr_range_l1(j,5)./100,'x');
grid
xlabel('CNR [dBHz]');
ylabel('Sigma [M]')
title('L1 pseudorange DD Res');

max_l1 = max(cnr_range_l1(j,5)./100);

% Form the weighed average
weighted_range_l1 = weight(0:3:60, cnr_range_l1, j, 5);
hold
plot(0:3:60,weighted_range_l1(:,2)./100,'r')

% L2 range Sigma
subplot(2,2,2)
j = find(cnr_range_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_range_l2(j,2),cnr_range_l2(j,5)./100,'x');
grid
xlabel('CNR [dBHz]');
ylabel('Sigma [M]')
title('L2 pseudorange DD Res');

max_l2 = max(cnr_range_l2(j,5)./100);

% Form the weighed average
weighted_range_l2 = weight(0:3:60, cnr_range_l2, j, 5);
hold
plot(0:3:60,weighted_range_l2(:,2)./100,'r')

% Now make sure L1 & L2 sigma graphs have the same scale.

m = max([max_l1,max_l2]);
m = ceil(m);

if(AutoScale == TRUE)
  % Need to make sure we can create a valid axis
  if(m ~= 0)
    axis([10 50 0 m]);  % L2
    subplot(2,2,1);
    axis([30 60 0 m]);  % L1
  end
else
  axis([10 50 0 CodeMax]);
  subplot(2,2,1);
  axis([30 60 0 CodeMax]);
end

% L1 range mean
subplot(2,2,3)
j = find(cnr_range_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_range_l1(j,2),cnr_range_l1(j,4)./100,'x');
grid
xlabel('CNR [dBHz]')
ylabel('Mean [M]')
title('L1 pseudorange DD Res');
% Form the weighed average
weighted_range_l1 = weight(0:3:60, cnr_range_l1, j, 4);
hold
plot(0:3:60,weighted_range_l1(:,2)./100,'r')


max_l1 = max(cnr_range_l1(j,4)./100);
min_l1 = min(cnr_range_l1(j,4)./100);

% L2 range mean
subplot(2,2,4)
j = find(cnr_range_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_range_l2(j,2),cnr_range_l2(j,4)./100,'x');
grid
xlabel('CNR [dBHz]');
ylabel('Mean [M]')
title('L2 pseudorange DD Res');
% Form the weighed average
weighted_range_l2 = weight(0:3:60, cnr_range_l2, j, 4);
hold
plot(0:3:60,weighted_range_l2(:,2)./100,'r')


max_l2 = max(cnr_range_l2(j,4)./100);
min_l2 = min(cnr_range_l2(j,4)./100);

max_mean = max([max_l1,max_l2]);
min_mean = min([min_l1,min_l2]);

%
% Scale to closest 0.5m
%
max_mean = ceil(max_mean*2)/2; 
min_mean = floor(min_mean*2)/2;

if(AutoScale == TRUE)
  % Need to make sure we can create a valid axis
  if(min_mean ~= max_mean)
    axis([10 50 min_mean max_mean]);  % L2
    subplot(2,2,3);
    axis([30 60 min_mean max_mean]);  % L1
  end
else
  axis([10 50 -CodeMax CodeMax]);
  subplot(2,2,3);
  axis([30 60 -CodeMax CodeMax]);
end

suptitle(Title);

%
% Plot the distribution of the counts.
%

fig(5) = figure('Visible','off');
subplot(2,1,1)
j = find(cnr_range_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_range_l1(j,2),cnr_range_l1(j,3),'x');

% Form the weighed average
k = 1;
for i=0:3:60
  m = find(cnr_range_l1(j,2) == i);
  if(m>0)
    av_sum_l1(k) = sum(cnr_range_l1(j(m),3))/length(m);
  else
    av_sum_l1(k) = NaN;
  end
  k = k + 1;
end
hold
plot(0:3:60,av_sum_l1,'r')

A    = axis;
A(1) = 30;
A(2) = 60;
axis(A);
grid
xlabel('CNR [dBHz]');
ylabel('Number or Epochs');
title('L1 No. of Epochs Versus CNR');

subplot(2,1,2)
j = find(cnr_range_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(cnr_range_l2(j,2),cnr_range_l2(j,3),'x');

% Form the weighed average
k = 1;
for i=0:3:60
  m = find(cnr_range_l2(j,2) == i);
  if(m>0)
    av_sum_l2(k) = sum(cnr_range_l2(j(m),3))/length(m);
  else
    av_sum_l2(k) = NaN;
  end
  k = k + 1;
end
hold
plot(0:3:60,av_sum_l2,'r')

A    = axis;
A(1) = 10;
A(2) = 50;
axis(A);
grid
xlabel('CNR [dBHz]');
ylabel('Number or Epochs');
title('L2 No. of Epochs Versus CNR');
suptitle(Title);

%
% Now plot the statistics on a per satellite basis
%

fig(6) = figure('Visible','off');

% L1 phase Sigma
subplot(2,2,1)
j = find(stats_phase_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(stats_phase_l1(j,2),stats_phase_l1(j,5),'x');
grid
xlabel('Satellite [PRN #]')
ylabel('Sigma [mcycle]')
title('L1 Phase DD Res');


% L2 phase Sigma
subplot(2,2,2)
j = find(stats_phase_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(stats_phase_l2(j,2),stats_phase_l2(j,5),'x');
grid
xlabel('Satellite [PRN #]')
ylabel('Sigma [mcycle]')
title('L2 Phase DD Res');


% L1 phase mean
subplot(2,2,3)
j = find(stats_phase_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(stats_phase_l1(j,2),stats_phase_l1(j,4),'x');
grid
xlabel('Satellite [PRN #]')
ylabel('Mean [mcycle]')
title('L1 Phase DD Res');

% L2 phase mean
subplot(2,2,4)
j = find(stats_phase_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(stats_phase_l2(j,2),stats_phase_l2(j,4),'x');
grid
xlabel('Satellite [PRN #]')
ylabel('Mean [mcycle]')
title('L2 Phase DD Res');
suptitle(Title);

fig(7) = figure('Visible','off');
% L1 range Sigma
subplot(2,2,1)
j = find(stats_range_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(stats_range_l1(j,2),stats_range_l1(j,5)./100,'x');
grid
xlabel('Satellite [PRN #]')
ylabel('Sigma [M]')
title('L1 pseudorange DD Res');

% L2 range Sigma
subplot(2,2,2)
j = find(stats_range_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(stats_range_l2(j,2),stats_range_l2(j,5)./100,'x');
grid
xlabel('Satellite [PRN #]')
ylabel('Sigma [M]')
title('L2 pseudorange DD Res');

% L1 range mean
subplot(2,2,3)
j = find(stats_range_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(stats_range_l1(j,2),stats_range_l1(j,4)./100,'x');
grid
xlabel('Satellite [PRN #]')
ylabel('Mean [M]')
title('L1 pseudorange DD Res');

% L2 range mean
subplot(2,2,4)
j = find(stats_range_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
plot(stats_range_l2(j,2),stats_range_l2(j,4)./100,'x');
grid
xlabel('Satellite [PRN #]')
ylabel('Mean [M]')
title('L2 pseudorange DD Res');
suptitle(Title);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now create a figure with L1 & L2 CNRs against
% elevation angle.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fig(8) = figure('Visible','off');
subplot(2,1,1)
plot(elev_range_l1(:,2),elev_range_l1(:,6),'x');
%grid
%xlabel('Elev. Angle [Deg]');
%ylabel('CNR [dBHz]')
%title('L1 CNR');

hold
plot(elev_range_l2(:,2),elev_range_l2(:,6),'xr');
grid
xlabel('Elev. Angle [Deg]');
ylabel('CNR [dBHz]')
title('L1/L2 CNR');
axis([0 90 0 60])

k = 1;
for i=0:5:90
  m = find(elev_range_l1(:,2) == i);
  if(length(m) > 0) 
    weighted_cnr_l1(k,1) = sum(elev_range_l1(m,3));
    weighted_cnr_l1(k,2) = sum((elev_range_l1(m,3).*  ...
                                  elev_range_l1(m,6))./ ...
                                  weighted_cnr_l1(k,1));
  else
    weighted_cnr_l1(k,1:2) = [NaN,NaN];
  end
  k = k + 1;
end

k = 1;
for i=0:5:90
  m = find(elev_range_l2(:,2) == i);
  if(length(m) > 0) 
    weighted_cnr_l2(k,1) = sum(elev_range_l2(m,3));
    weighted_cnr_l2(k,2) = sum((elev_range_l2(m,3).*  ...
                                  elev_range_l2(m,6))./ ...
                                  weighted_cnr_l2(k,1));
  else
    weighted_cnr_l2(k,1:2) = [NaN,NaN];
  end
  k = k + 1;
end

plot(0:5:90,weighted_cnr_l1(:,2)-weighted_cnr_l2(:,2),'m-.')
plot(0:5:90,weighted_cnr_l1(:,2)-weighted_cnr_l2(:,2),'mx')
plot(0:5:90,weighted_cnr_l1(:,2),'b')
plot(0:5:90,weighted_cnr_l2(:,2),'r')

legend('L1','L2','L1-L2',2)

%
% Now plot elevation against CNR
%

subplot(2,1,2)
plot(cnr_range_l1(:,2),cnr_range_l1(:,6),'x');
hold
plot(cnr_range_l2(:,2),cnr_range_l2(:,6),'xr');
grid
xlabel('CNR [dBHz]');
ylabel('Elevation [Deg]')
title('Elevation(CNR)');
axis([0 60 0 90])

k = 1;
for i=0:3:60
  m = find(cnr_range_l1(:,2) == i);
  if(length(m) > 0) 
    weighted_el_l1(k,1) = sum(cnr_range_l1(m,3));
    weighted_el_l1(k,2) = sum((cnr_range_l1(m,3).*  ...
                                  cnr_range_l1(m,6))./ ...
                                  weighted_el_l1(k,1));
  else
    weighted_el_l1(k,1:2) = [NaN,NaN];
  end
  k = k + 1;
end

k = 1;
for i=0:3:60
  m = find(cnr_range_l2(:,2) == i);
  if(length(m) > 0) 
    weighted_el_l2(k,1) = sum(cnr_range_l2(m,3));
    weighted_el_l2(k,2) = sum((cnr_range_l2(m,3).*  ...
                                  cnr_range_l2(m,6))./ ...
                                  weighted_el_l2(k,1));
  else
    weighted_el_l2(k,1:2) = [NaN,NaN];
  end
  k = k + 1;
end

plot(0:3:60,weighted_el_l1(:,2),'b')
plot(0:3:60,weighted_el_l2(:,2),'r')
legend('L1','L2',2)




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Now plot the outliers....
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

min_time = NaN;
max_time = NaN;
min_pseudo = NaN; 
max_pseudo = NaN;
min_phase = NaN;
max_phase = NaN;
suptitle(Title);

fig(9) = figure('Visible','off');
for j=0:3
  subplot(2,2,j+1)
  if(j==0)
    i = find(stats_range_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
    stats = stats_range_l1(i,:); 
  elseif(j==1)
    i = find(stats_range_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
    stats = stats_range_l2(i,:); 
  elseif(j==2)
    i = find(stats_phase_l1(:,3) > ACCEPTABLE_SAMPLE_SET);
    stats = stats_phase_l1(i,:); 
  elseif(j==3)
    i = find(stats_phase_l2(:,3) > ACCEPTABLE_SAMPLE_SET);
    stats = stats_phase_l2(i,:); 
  end

  weeks = 0;
  time_diff = diff(stats(:,1));
  stats_out(1,:)  = stats(1,:);
  for i=1:(length(stats(:,1))-1)
    if(time_diff(i) < -302400) % Half a week
      weeks = weeks + 1;
    end
    if(time_diff(i) > 302400) % Half a week
      weeks = weeks - 1;
    end
    stats_out(i+1,:) = stats(i+1,:);
    stats_out(i+1,1) = stats_out(i+1,1) + weeks*3600*24*7;
  end

  stats_out = sortrows(stats_out,1);

  if(j<=1)
    plot(stats_out(:,1)./86400,stats_out(:,4)./100);
    hold
    plot(stats_out(:,1)./86400,stats_out(:,6)./100,'r');
    plot(stats_out(:,1)./86400,stats_out(:,7)./100,'r');
    min_pseudo = min([min(stats_out(:,6)),min_pseudo]);
    max_pseudo = max([max(stats_out(:,7)),max_pseudo]);
  else
    plot(stats_out(:,1)./86400,stats_out(:,4));
    hold
    plot(stats_out(:,1)./86400,stats_out(:,6),'r');
    plot(stats_out(:,1)./86400,stats_out(:,7),'r');
    min_phase = min([min(stats_out(:,6)),min_phase]);
    max_phase = max([max(stats_out(:,7)),max_phase]);
  end

  if(j==0)
    title('L1 Pseudorange DD Min/Mean/Max');
    ylabel('L1 Pseudorange DD [M]');
  elseif(j==1)
    title('L2 Pseudorange DD Min/Mean/Max');
    ylabel('L2 Pseudorange DD [M]');
  elseif(j==2)
    title('L1 Phase DD Min/Mean/Max');
    ylabel('L1 Phase DD [mcycle]');
  elseif(j==3)
    title('L2 Phase DD Min/Mean/Max');
    ylabel('L2 Phase DD [mcycle]');
  end
  
  grid
  xlabel('Time of Week [Day]');

  min_time = min([min(stats_out(:,1)),min_time]);
  max_time = max([max(stats_out(:,1)),max_time]);
end

min_time = floor(min_time/86400);
max_time = ceil(max_time/86400);
min_pseudo = floor(min_pseudo/500)*5;
max_pseudo = ceil(max_pseudo/500)*5;
min_phase = floor(min_phase/250)*250;
max_phase = ceil(max_phase/250)*250;

for j=0:3
  subplot(2,2,j+1)
  if(j<=1)
    axis([min_time max_time min_pseudo max_pseudo]);
  else
    axis([min_time max_time min_phase max_phase]);
  end
end

% suptitle('Observable Outliers');
suptitle(Title);

if(PNGPlot == TRUE)
  if( Compiled == FALSE )
    print -dpng -f1 phase1.png
    print -dpng -f2 range1.png
    print -dpng -f3 phase2.png
    print -dpng -f4 range2.png
    print -dpng -f5 cnr1.png
    print -dpng -f6 phase3.png
    print -dpng -f7 range3.png
    print -dpng -f8 cnr2.png
    print -dpng -f9 outlier.png

    % print -dtiff -f1 phase1.tif
    % print -dtiff -f2 range1.tif
    % print -dtiff -f3 phase2.tif
    % print -dtiff -f4 range2.tif
    % print -dtiff -f5 cnr1.tif
    % print -dtiff -f6 phase3.tif
    % print -dtiff -f7 range3.tif
    % print -dtiff -f8 cnr2.tif
    % print -dtiff -f9 outlier.tif
  else
    % png files not supported by the Matlab compiler.
    print -dbitmap -f1 phase1.bmp
    print -dbitmap -f2 range1.bmp
    print -dbitmap -f3 phase2.bmp
    print -dbitmap -f4 range2.bmp
    print -dbitmap -f5 cnr1.bmp
    print -dbitmap -f6 phase3.bmp
    print -dbitmap -f7 range3.bmp
    print -dbitmap -f8 cnr2.bmp
    print -dbitmap -f9 outlier.bmp

    % All bitmap files have been created.
    % Translate .bmp files to .prn files an delete .bmp files.
    dos( 'bmp2png phase*.bmp' ) ;
    dos( 'bmp2png range*.bmp' ) ;
    dos( 'bmp2png cnr*.bmp' ) ;
    dos( 'bmp2png outlier.bmp' ) ;
    dos( 'erase phase*.bmp' ) ;
    dos( 'erase range*.bmp' ) ;
    dos( 'erase cnr*.bmp' ) ;
    dos( 'erase outlier.bmp' ) ;
  end

  %
  % Finally create an HTML file we can use to publish the results to the
  % Intranet.
  %

  fid  = fopen('trend.htm','w');
  % Start the HTML document
  fprintf(fid,'<HTML>\n');
  fprintf(fid,'<BODY BGCOLOR= "#E0E0E0" TEXT="000000" BACKGROUND="backgrnd.gif">\n');
  fprintf(fid,'<HR>\n');
  fprintf(fid,'<a href="phase1.png">Phase Vs Elevation</A><BR>\n');
  fprintf(fid,'<a href="range1.png">range Vs Elevation</A><BR>\n');
  fprintf(fid,'<a href="phase2.png">Phase Vs CNR</A><BR>\n');
  fprintf(fid,'<a href="range2.png">range Vs CNR</A><BR>\n');
  fprintf(fid,'<a href="phase3.png">Phase Vs PRN #</A><BR>\n');
  fprintf(fid,'<a href="range3.png">range Vs PRN #</A><BR>\n');
  fprintf(fid,'<a href="cnr2.png">CNR Vs Elevation</A><BR>\n');
  fprintf(fid,'<a href="outlier.png">Phase/Range Outliers</A><BR>\n');
  fprintf(fid,'<HR>\n');
  fprintf(fid,'Created using Matlab: $Id: plot_obs.m,v 1.7 2005/04/18 19:46:46 vnadkar Exp $\n');
  fprintf(fid,'<HR>\n');
  fprintf(fid,'</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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% $Log: plot_obs.m,v $
% Revision 1.7  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.6  2001/04/02 22:04:53  stuart
% o Updated with latest QA scripts
%
% Revision 1.5  2000/09/24 19:49:27  timo
% Remove highlight markers, >>>>> (changes between this and
% last rev). No functional changes.
%
% Revision 1.4  2000/09/22 18:52:20  timo
% Update to allow Matlab compilation (StandAlone application).
% Updates highlighed with >>>>> (to be removed).
%
% Revision 1.3  2000/08/21 13:57:34  stuart
% o Correct $Id
% o Add $Log
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

