function PlotTrend(phaseDDiffL1, ...
                   phaseDDiffL2, ...
                   phaseDDiffL5_L1P, ...
                   rangeDDiffL1, ...
                   rangeDDiffL2, ...
                   rangeDDiffL5_L1P, ...
                   Elev, ...
                   CNRL1, ...
                   CNRL2, ...
                   CNRL5_L1P, ...
                   timeStr, ...
                   PlotPNG, ...
                   CarrMax, ...
                   CodeMax, ...
                   L2Present, ...
                   L5_L1P_Present, ...
                   system);

%*********************************************************************
% $Id: PlotTrend.m,v 1.10 2006/04/27 00:22:03 vnadkar Exp $
%
% Computes and plots the following:
% * DD residual sigma versus elevation angle
% * DD residual sigma versus C/No
% * Histogram of the DD residual errors
%
%*********************************************************************
 
TRUE  = 1;
FALSE = 0;

%
% Plot Histograms of the data
%

if(system == 0)
  svs = 32;
else
  svs = 24;
end

% Convert to millicycles (GPS) / millimetres (GLONASS)
obsPhaseL1 = phaseDDiffL1(:,2:(svs+1))./10;  
obsRangeL1 = rangeDDiffL1(:,2:(svs+1))./100; % Convert to centimetres

if(L2Present == TRUE)
  % Convert to millicycles (GPS) / millimetres (GLONASS)
  obsPhaseL2 = phaseDDiffL2(:,2:(svs+1))./10;
  obsRangeL2 = rangeDDiffL2(:,2:(svs+1))./100; % Convert to centimetres
end

if(L5_L1P_Present == TRUE)
  % Convert to millicycles (GPS) / millimetres (GLONASS)
  obsPhaseL5_L1P = phaseDDiffL5_L1P(:,2:(svs+1))./10;
  obsRangeL5_L1P = rangeDDiffL5_L1P(:,2:(svs+1))./100; % Convert to centimetres
end

%
% GPS only processing
%
if(system == 0)
  ftrend  = fopen('trend.txt','w');
  dThres = 50; % Need at least this many points in each bin
  k = 1;
  for j=2.5:5:87.5
    i = find( (Elev >= j) &  (Elev < (j+5)) );
    phaseEl(k,1)  = j+2.5; % Centre of the bin
    rangeEl(k,1)  = j+2.5; % Centre of the bin

    if(length(i) >  dThres)
      % L1 Phase
      m = find( isfinite(obsPhaseL1(i)) == 1);

      % Remove any bad outliers +/- 1000 millicycles.
      n = find( abs(obsPhaseL1(i(m))) < 1000);
      if(length(n) > dThres)
        phaseEl(k,2)  = std(obsPhaseL1(i(m(n))));
      else
        phaseEl(k,2)  = NaN;
      end

      if(L2Present == TRUE)
        %L2 Phase
        m = find( isfinite(obsPhaseL2(i)) == 1);
    
        % Remove any bad outliers +/- 1000 millicycles.
        n = find( abs(obsPhaseL2(i(m))) < 1000);
        if(length(n) > dThres)
          phaseEl(k,3)  = std(obsPhaseL2(i(m(n))));
        else
          phaseEl(k,3)  = NaN;
        end
      else
        phaseEl(k,3)  = NaN;
      end

      % L1 Range
      m = find( isfinite(obsRangeL1(i)) == 1);
    
      % Remove any bad outliers +/- 10 m
      n = find( abs(obsRangeL1(i(m))) < 1000);
    
      if(length(n) > dThres)
        rangeEl(k,2)  = std(obsRangeL1(i(m(n))));
      else
        rangeEl(k,2)  = NaN;
      end

      if(L2Present == TRUE)
        %L2 Range
        m = find( isfinite(obsRangeL2(i)) == 1);
    
        % Remove any bad outliers +/- 10 m
        n = find( abs(obsRangeL2(i(m))) < 1000);
        if(length(n) > dThres)
          rangeEl(k,3)  = std(obsRangeL2(i(m(n))));
        else
          rangeEl(k,3)  = NaN;
        end
      else
        rangeEl(k,3)  = NaN;
      end
    else
      phaseEl(k,2)  = NaN;
      phaseEl(k,3)  = NaN;
      rangeEl(k,2)  = NaN;
      rangeEl(k,3)  = NaN;
    end

    fprintf(ftrend,'%.1f %.3f %.3f %.3f %.3f\n', ...
           phaseEl(k,1), ...
           phaseEl(k,2), ...
           phaseEl(k,3), ...
           rangeEl(k,2), ...
           rangeEl(k,3));
  
    k = k + 1;
  end
  fclose(ftrend);

  ftrend  = fopen('trendCno.txt','w');
  dThres = 50; % Need at least this many points in each bin

  k = 1;
  CNOBinSize = 2;

  if(length(CNRL1) > 0)
    for j=1:CNOBinSize:60
      i = find( (CNRL1 >= j) &  (CNRL1 < (j+CNOBinSize)) );

      phaseCNo(k,1)  = j+CNOBinSize/2; % Centre of the bin
      rangeCNo(k,1)  = j+CNOBinSize/2; % Centre of the bin

      if(length(i) >  dThres)
        % L1 Phase
        m = find( isfinite(obsPhaseL1(i)) == 1);

        % Remove any bad outliers +/- 1000 millicycles.
        n = find( abs(obsPhaseL1(i(m))) < 1000);
        if(length(n) > dThres)
          phaseCNo(k,2)  = std(obsPhaseL1(i(m(n))));
        else
          phaseCNo(k,2)  = NaN;
        end
    
        % L1 Range
        m = find( isfinite(obsRangeL1(i)) == 1);
    
        % Remove any bad outliers +/- 10 m
        n = find( abs(obsRangeL1(i(m))) < 1000);
    
        if(length(n) > dThres)
          rangeCNo(k,2)  = std(obsRangeL1(i(m(n))));
        else
          rangeCNo(k,2)  = NaN;
        end
      else
        phaseCNo(k,2)  = NaN;
        rangeCNo(k,2)  = NaN;
      end

      if(L2Present == TRUE)
  
        i = find( (CNRL2 >= j) &  (CNRL2 < (j+CNOBinSize)) );

        if(length(i) >  dThres)
          %L2 Phase
          m = find( isfinite(obsPhaseL2(i)) == 1);
    
          % Remove any bad outliers +/- 1000 millicycles.
          n = find( abs(obsPhaseL2(i(m))) < 1000);
          if(length(n) > dThres)
            phaseCNo(k,3)  = std(obsPhaseL2(i(m(n))));
          else
            phaseCNo(k,3)  = NaN;
          end

          %L2 Range
          m = find( isfinite(obsRangeL2(i)) == 1);
    
          % Remove any bad outliers +/- 10 m
          n = find( abs(obsRangeL2(i(m))) < 1000);
          if(length(n) > dThres)
            rangeCNo(k,3)  = std(obsRangeL2(i(m(n))));
          else
            rangeCNo(k,3)  = NaN;
          end
        else
          phaseCNo(k,3)  = NaN;
          rangeCNo(k,3)  = NaN;
        end
      else
        phaseCNo(k,3)  = NaN;
        rangeCNo(k,3)  = NaN;
      end

      if(L5_L1P_Present == TRUE)
  
        i = find( (CNRL5_L1P >= j) &  (CNRL5_L1P < (j+CNOBinSize)) );

        if(length(i) >  dThres)
          %L5 Phase
          m = find( isfinite(obsPhaseL5_L1P(i)) == 1);
    
          % Remove any bad outliers +/- 1000 millicycles.
          n = find( abs(obsPhaseL5_L1P(i(m))) < 1000);
          if(length(n) > dThres)
            phaseCNo(k,4)  = std(obsPhaseL5_L1P(i(m(n))));
          else
            phaseCNo(k,4)  = NaN;
          end

          %L5 Range
          m = find( isfinite(obsRangeL5_L1P(i)) == 1);
    
          % Remove any bad outliers +/- 10 m
          n = find( abs(obsRangeL5_L1P(i(m))) < 1000);
          if(length(n) > dThres)
            rangeCNo(k,4)  = std(obsRangeL5_L1P(i(m(n))));
          else
            rangeCNo(k,4)  = NaN;
          end
        else
          phaseCNo(k,4)  = NaN;
          rangeCNo(k,4)  = NaN;
        end
      else
        phaseCNo(k,4)  = NaN;
        rangeCNo(k,4)  = NaN;
      end


      fprintf(ftrend,'%.1f %.3f %.3f %.3f %.3f\n', ...
           phaseCNo(k,1), ...
           phaseCNo(k,2), ...
           phaseCNo(k,3), ...
           rangeCNo(k,2), ...
           rangeCNo(k,3));
  
      k = k + 1;
    end
  end 
  fclose(ftrend);
end % system == 0

%
% Now plot data as a histogram
%

% L1 Phase
figure('Visible','off');

if(system == 0)
  disp('GPS: Creating L1 Phase distribution');
  myString = 'GPS: L1 Phase Double Difference: \mu=';
  myUnits = 'mcycle';
else
  disp('GLN: Creating L1-C/A Phase distribution');
  myString = 'GLN: L1-C/A Phase Double Difference: \mu=';
  myUnits = 'mm';
end

i = find( (obsPhaseL1 < CarrMax/2) & (obsPhaseL1 > -CarrMax/2));
obsPhaseL1 = obsPhaseL1(i);

if ~isempty(obsPhaseL1)
  [Chi, dof, Prob, muL1Phase, sigL1Phase, Num] =  ...
    plotHist(obsPhaseL1,10,-CarrMax/2,CarrMax/2,FALSE);

  title([myString ,num2str(muL1Phase,'%.1f'),...
    myUnits, '; \sigma=', num2str(sigL1Phase,'%.1f'), ...
    myUnits, '; Num. Epochs=',num2str(Num)]);

  if(system == 0)
    xlabel(['Error [millicycle]    ',timeStr],'Interpreter','none');
    if(PlotPNG == TRUE)
      print('-dpng',['L1PhaseDist.png']);
      close
    end
  else
    xlabel(['Error [mm]    ',timeStr],'Interpreter','none');
    if(PlotPNG == TRUE)
      print('-dpng',['GLNL1CAPhaseDist.png']);
      close
    end
  end
else
  if(PlotPNG == TRUE)
    figure;
    if (system == 0)
      title('No L1 Phase Observables');
      print('-dpng',['L1PhaseDist.png']);
    else
      title('No GLN L1 Phase Observables');
      print('-dpng',['GLNL1CAPhaseDist.png']);
    end
    close;
  end  
end

if(L2Present == TRUE)
  % L2 Phase
  figure('Visible','off');

  if(system == 0)
    disp('GPS: Creating L2 Phase distribution');
    myString = 'GPS: L2 Phase Double Difference: \mu=';
  else
    disp('GLN: Creating L2 Phase distribution');
    myString = 'GLN: L2 Phase Double Difference: \mu=';
  end

  i = find( (obsPhaseL2 < CarrMax/2) & (obsPhaseL2 > -CarrMax/2));
  obsPhaseL2 = obsPhaseL2(i);

  if ~isempty(obsPhaseL2)
    [Chi, dof, Prob, muL2Phase, sigL2Phase, Num] = ...
      plotHist(obsPhaseL2,10,-CarrMax/2,CarrMax/2,FALSE);

    title([myString,num2str(muL2Phase,'%.1f'),...
      myUnits, '; \sigma=', num2str(sigL2Phase,'%.1f'), ...
      myUnits, '; Num. Epochs=',num2str(Num)]);

    if(system == 0)
      xlabel(['Error [millicycle]    ',timeStr],'Interpreter','none');

      if(PlotPNG == TRUE)
        print('-dpng',['L2PhaseDist.png']);
        close
      end
    else
      xlabel(['Error [mm]    ',timeStr],'Interpreter','none');

      if(PlotPNG == TRUE)
        print('-dpng',['GLNL2PhaseDist.png']);
        close
      end
    end
  else
    if(PlotPNG == TRUE)
      figure;
      if (system == 0)
        title('No L2 Phase Observables');
        print('-dpng',['L2PhaseDist.png']);
      else
        title('No GLN L2 Phase Observables');
        print('-dpng',['GLNL2CAPhaseDist.png']);
      end
      close;
    end
  end
end

if(L5_L1P_Present == TRUE)
  % L5 Phase
  figure('Visible','off');
  
  if(system == 0)
    disp('GPS: Creating L5 Phase distribution');
    myString = 'GPS: L5 Phase Double Difference: \mu=';
  else
    disp('GLN: Creating L1-P Phase distribution');
    myString = 'GLN: L1-P Phase Double Difference: \mu=';
  end

  i = find( (obsPhaseL5_L1P < CarrMax/2) & (obsPhaseL5_L1P > -CarrMax/2));
  obsPhaseL5_L1P = obsPhaseL5_L1P(i);

  if ~isempty(obsPhaseL5_L1P)
    [Chi, dof, Prob, muL5Phase, sigL5Phase, Num] = ...
      plotHist(obsPhaseL5_L1P,10,-CarrMax/2,CarrMax/2,FALSE);

    title([myString,num2str(muL5Phase,'%.1f'),...
      myUnits, '; \sigma=', num2str(sigL5Phase,'%.1f'), ...
      myUnits, '; Num. Epochs=',num2str(Num)]);

    if(system == 0)
      xlabel(['Error [millicycle]    ',timeStr],'Interpreter','none');

      if(PlotPNG == TRUE)
        print('-dpng',['L5PhaseDist.png']);
        close
      end
    else
      xlabel(['Error [mm]    ',timeStr],'Interpreter','none');

      if(PlotPNG == TRUE)
        print('-dpng',['GLNL1PPhaseDist.png']);
        close
      end
    end
  else
    if(PlotPNG == TRUE)
      figure;
      if (system == 0)
        title('No L5 Phase Observables');
        print('-dpng',['L5PhaseDist.png']);
      else
        title('No GLN L1 P Phase Observables');
        print('-dpng',['GLNL1PPhaseDist.png']);
      end
      close;
    end
  end

end

% L1 Range
figure('Visible','off');

if(system == 0)
  disp('GPS: Creating L1 Range distribution');
  myString = 'GPS: L1 Range Double Difference: \mu=';
else
  disp('GLN: Creating L1-C/A Range distribution');
  myString = 'GLN: L1-C/A Range Double Difference: \mu=';
end

i = find( (obsRangeL1 < 100*CodeMax) & (obsRangeL1 > -100*CodeMax));
obsRangeL1 = obsRangeL1(i);

if ~isempty(obsRangeL1)
  [Chi, dof, Prob, muL1Range, sigL1Range, Num] = ...
    plotHist(obsRangeL1,10,-100*CodeMax,100*CodeMax,FALSE);

  title([myString, ...
    num2str(muL1Range/100,'%.2f'), ...
    'm; \sigma=', num2str(sigL1Range/100,'%.2f'), ...
    'm; Num. Epochs=',num2str(Num)]);

  xlabel(['Error [cm]    ',timeStr],'Interpreter','none');

  if(PlotPNG == TRUE)
    if(system == 0)
      print('-dpng',['L1RangeDist.png']);
    else
      print('-dpng',['GLNL1CARangeDist.png']);
    end
    close
  end
else
  if(PlotPNG == TRUE)
    figure;
    if(system == 0)
      title('No L1 Range Observables');
      print('-dpng',['L1RangeDist.png']);
    else
      title('No GLN L1 Range Observables');
      print('-dpng',['GLNL1CARangeDist.png']);
    end
    close
  end
end

if(L2Present)
  % L2 Range
  figure('Visible','off');
  
  if(system == 0)
    disp('GPS: Creating L2 Range distribution');
    myString = 'GPS: L2 Range Double Difference: \mu=';
  else
    disp('GLN: Creating L2 Range distribution');
    myString = 'GLN: L2 Range Double Difference: \mu=';
  end

  i = find( (obsRangeL2 < 100*CodeMax) & (obsRangeL2 > -100*CodeMax));
  obsRangeL2 = obsRangeL2(i);

  if ~isempty(obsRangeL2)
    [Chi, dof, Prob, muL2Range, sigL2Range, Num] = ...
      plotHist(obsRangeL2,10,-100*CodeMax,100*CodeMax,FALSE);

    title([myString, ...
      num2str(muL2Range/100,'%.2f'), ...
      'm; \sigma=', num2str(sigL2Range/100,'%.2f'), ...
      'm; Num. Epochs=',num2str(Num)]);

    xlabel(['Error [cm]    ',timeStr],'Interpreter','none');

    if(PlotPNG == TRUE)
      if(system == 0)
        print('-dpng',['L2RangeDist.png']);
      else
        print('-dpng',['GLNL2RangeDist.png']);
      end
      close
    end
  else
    if(PlotPNG == TRUE)
      figure;
      if(system == 0)
        title('No L2 Range Observables');
        print('-dpng',['L2RangeDist.png']);
      else
        title('No GLN L2CA Range Observables');
        print('-dpng',['GLNL2CARangeDist.png']);
      end
      close
    end
  end
end

if(L5_L1P_Present)
  % L5 Range
  figure('Visible','off');

  if(system == 0)
    disp('GPS: Creating L5 Range distribution');
    myString = 'GPS: L5 Range Double Difference: \mu=';
  else
    disp('GLN: Creating L1-P Range distribution');
    myString = 'GLN: L1-P Range Double Difference: \mu=';
  end
 
  i = find( (obsRangeL5_L1P < 100*CodeMax) & (obsRangeL5_L1P > -100*CodeMax));
  obsRangeL5_L1P = obsRangeL5_L1P(i);

  if ~isempty(obsRangeL5_L1P)
    [Chi, dof, Prob, muL5Range, sigL5Range, Num] = ...
      plotHist(obsRangeL5_L1P,10,-100*CodeMax,100*CodeMax,FALSE);

    title([myString, ...
      num2str(muL5Range/100,'%.2f'), ...
      'm; \sigma=', num2str(sigL5Range/100,'%.2f'), ...
      'm; Num. Epochs=',num2str(Num)]);

    xlabel(['Error [cm]    ',timeStr],'Interpreter','none');

    if(PlotPNG == TRUE)
      if(system == 0)
        print('-dpng',['L5RangeDist.png']);
      else
        print('-dpng',['GLNL1PRangeDist.png']);
      end
      close
    end
  else
    if(PlotPNG == TRUE)
      figure;
      if(system == 0)
        title('No L5 Range Observables');
        print('-dpng',['L5RangeDist.png']);
      else
        title('No GLN L1P Range Observables');
        print('-dpng',['GLNL1PRangeDist.png']);
      end
      close
    end
  end    
end


%
% GPS Only
% 
if(system == 0)
  %
  % Now plot the standard deviation as a function of elevation
  % 

  figure('Visible','off');
  x(1) = 0;
  x(2) = 90;

  if(L2Present == TRUE)
    subplot(2,2,1)
  else
    subplot(2,1,1)
  end
  plot(phaseEl(:,1), phaseEl(:,2),'r-x');
  y(1) = sigL1Phase;
  y(2) = sigL1Phase;
  H = line(x,y);
  set(H,'color','b');
  set(H,'LineWidth',2);
  A = axis;
  A(1) = 0;
  A(2) = 90;
  axis(A)
  grid;
  legend('\sigma - Versus Elev','\sigma - All Data');
  ylabel('L1 Phase \sigma [mcycle]');
  xlabel('Elevation Angle [degrees]');
  title(['L1 Phase: \sigma=',num2str(sigL1Phase,'%.1f'), ...
         'mcycle; \mu=',num2str(muL1Phase,'%.1f'),'mcycle'])

  if(L2Present == TRUE)       
    subplot(2,2,2)
    plot(phaseEl(:,1), phaseEl(:,3),'r-x');
    y(1) = sigL2Phase;
    y(2) = sigL2Phase;
    H = line(x,y);
    set(H,'color','b');
    set(H,'LineWidth',2);
    B = axis;
    B(1) = 0;
    B(2) = 90;

    if(B(3) < A(3))
      ymin = B(3);
    else
      ymin = A(3);
    end
    B(3) = ymin;

    if(B(4) > A(4))
      ymax = B(4);
    else
      ymax = A(4);
    end
    B(4) = ymax;

    axis(B)
    grid;
    ylabel('L2 Phase \sigma [mcycle]');
    xlabel('Elevation Angle [degrees]');
    title(['L2 Phase: \sigma=',num2str(sigL2Phase,'%.1f'), ...
           'mcycle; \mu=',num2str(muL2Phase,'%.1f'),'mcycle'])
    subplot(2,2,1)
    A(3) = ymin;
    A(4) = ymax;
    axis(A)
  end

  if(L2Present == TRUE)
    subplot(2,2,3)
  else
    subplot(2,1,2)
  end
  plot(rangeEl(:,1), rangeEl(:,2)./100,'r-x');
  y(1) = sigL1Range/100;
  y(2) = sigL1Range/100;
  H = line(x,y);
  set(H,'color','b');
  set(H,'LineWidth',2);
  A = axis;
  A(1) = 0;
  A(2) = 90;
  axis(A)
  grid;
  ylabel('L1 Range \sigma [m]');
  xlabel('Elevation Angle [degrees]');
  title(['L1 Range: \sigma=',num2str(sigL1Range/100,'%.2f'), ...
         'm; \mu=',num2str(muL1Range/100,'%.2f'),'m'])


  if(L2Present == TRUE)       
    subplot(2,2,4)
    plot(rangeEl(:,1), rangeEl(:,3)./100,'r-x');
    y(1) = sigL2Range/100;
    y(2) = sigL2Range/100;
    H = line(x,y);
    set(H,'color','b');
    set(H,'LineWidth',2);
    B = axis;
    B(1) = 0;
    B(2) = 90;

    if(B(3) < A(3))
      ymin = B(3);
    else
      ymin = A(3);
    end
    B(3) = ymin;

    if(B(4) > A(4))
      ymax = B(4);
    else
      ymax = A(4);
    end
    B(4) = ymax;

    axis(B)
    grid;
    ylabel('L2 Range \sigma [m]');
    xlabel('Elevation Angle [degrees]');
    title(['L2 Range: \sigma=',num2str(sigL2Range/100,'%.2f'), ...
           'm; \mu=',num2str(muL2Range/100,'%.2f'),'m'])

    subplot(2,2,3)
    A(3) = B(3);
    A(4) = B(4);
    axis(A)
  end

  suptitle('Double Difference Residuals Versus Elevation Angle: ',timeStr); 

  if(PlotPNG == TRUE)
    print('-dpng',['DDVersusElev.png']);
    close
  end


  if(length(CNRL1) > 0)
    figure('Visible','off');
    x(1) = 0;
    x(2) = 60;
    L1CNoMin = 20;
    L1CNoMax = 60;
    L2CNoMin = 0;
    L2CNoMax = 50;

    if(L2Present == TRUE)
      subplot(2,2,1)
    else
      subplot(2,1,1)
    end
    plot(phaseCNo(:,1), phaseCNo(:,2),'r-x');
    y(1) = sigL1Phase;
    y(2) = sigL1Phase;
    H = line(x,y);
    set(H,'color','b');
    set(H,'LineWidth',2);
    A = axis;
    A(1) = L1CNoMin;
    A(2) = L1CNoMax;
    axis(A)
    grid;
    legend('\sigma - Versus C/No','\sigma - All Data');
    ylabel('L1 Phase \sigma [mcycle]');
    xlabel('C/No [dBHz]');
    title(['L1 Phase: \sigma=',num2str(sigL1Phase,'%.1f'), ...
           'mcycle; \mu=',num2str(muL1Phase,'%.1f'),'mcycle'])

    if(L2Present == TRUE)
      subplot(2,2,2)
      plot(phaseCNo(:,1), phaseCNo(:,3),'r-x');
      y(1) = sigL2Phase;
      y(2) = sigL2Phase;
      H = line(x,y);
      set(H,'color','b');
      set(H,'LineWidth',2);
      B = axis;
      B(1) = L2CNoMin;
      B(2) = L2CNoMax;

      if(B(3) < A(3))
        ymin = B(3);
      else
        ymin = A(3);
      end
      B(3) = ymin;

      if(B(4) > A(4))
        ymax = B(4);
      else
        ymax = A(4);
      end
      B(4) = ymax;

      axis(B)
      grid;
      ylabel('L2 Phase \sigma [mcycle]');
      xlabel('C/No [dBHz]');
      title(['L2 Phase: \sigma=',num2str(sigL2Phase,'%.1f'), ...
             'mcycle; \mu=',num2str(muL2Phase,'%.1f'),'mcycle'])

      subplot(2,2,1)
      A(3) = ymin;
      A(4) = ymax;
      axis(A)
    end

    if(L2Present == TRUE)
      subplot(2,2,3)
    else
      subplot(2,1,2)
    end
    plot(rangeCNo(:,1), rangeCNo(:,2)./100,'r-x');
    y(1) = sigL1Range/100;
    y(2) = sigL1Range/100;
    H = line(x,y);
    set(H,'color','b');
    set(H,'LineWidth',2);
    A = axis;
    A(1) = L1CNoMin;
    A(2) = L1CNoMax;
    axis(A)
    grid;
    ylabel('L1 Range \sigma [m]');
    xlabel('C/No [dBHz]');
    title(['L1 Range: \sigma=',num2str(sigL1Range/100,'%.2f'), ...
           'm; \mu=',num2str(muL1Range/100,'%.2f'),'m'])

    if(L2Present == TRUE)       
      subplot(2,2,4)
      plot(rangeCNo(:,1), rangeCNo(:,3)./100,'r-x');
      y(1) = sigL2Range/100;
      y(2) = sigL2Range/100;
      H = line(x,y);
      set(H,'color','b');
      set(H,'LineWidth',2);
      B = axis;
      B(1) = L2CNoMin;
      B(2) = L2CNoMax;

      if(B(3) < A(3))
        ymin = B(3);
      else
        ymin = A(3);
      end
      B(3) = ymin;

      if(B(4) > A(4))
        ymax = B(4);
      else
        ymax = A(4);
      end
      B(4) = ymax;

      axis(B)
      grid;
      ylabel('L2 Range \sigma [m]');
      xlabel('C/No [dBHz]');
      title(['L2 Range: \sigma=',num2str(sigL2Range/100,'%.2f'), ...
             'm; \mu=',num2str(muL2Range/100,'%.2f'),'m'])
      subplot(2,2,3)
      A(3) = B(3);
      A(4) = B(4);
      axis(A)
    end

    suptitle('Double Difference Residuals Versus Rover C/No',timeStr);

    if(PlotPNG == TRUE)
      print('-dpng',['DDVersusCNo.png']);
      close
    end
  end

end % if(system == 0)

