function [h1,h2] = plotDdResids (residTypeStr, ymax, chooseSv)
% function plotDdResids (residTypeStr) 
% Plot double difference residuals of code or carrier phase 
% in various ways. 
% 
% residTypeStr - Residual Type string is one of
%                'code', 'phase', with the default being
%                'phase'
%
% ymax         - an optional input specifying the maximum limit
%                on the y axis. The y axis is set to have the 
%                limits [-ymax, ymax]
%
% Vivek Nadkarni
% July 7, 2004
% Trimble

% Plotting flags
printPlots = 1; % save plots in PNG format
plotFigC   = 1; % combined figure 1
plotFigC2  = 1; % combined figure 2
plotFigS2  = 1; % separate figures by SV
plotFigH   = 1; % histogram
printStats = 1; % print out statistics page

% Kludge
useTimeOffset = 0; % Set time offset to start with startTime. 
if exist('chooseSv','var')
  if ~isempty(chooseSv)
    printPlots = 0;
    plotFigC = 0;
    plotFigC2 = 0;
    plotFigH = 0;
    printStats = 0; 
    
    % This won't work for versions older than Matlab 7. Set
    % enableLinkAxes = 0 for older versions. 
    enableLinkAxes = 1; 
  end
end
  

% Legends in plotFigS are not correct. This is a legacy figure
% that is no longer used. 
plotFigS = 0; 

% Select plotting of code or phase
PHASE = 1;
CODE = 2; 
if (nargin == 0)
  residTypeStr = 'phase';
end
if strcmp(residTypeStr, 'phase')
  residType = PHASE;
  residTypeName = 'Phase';
elseif strcmp(residTypeStr, 'code')
  residType = CODE; 
  residTypeName = 'Code';
else
  error('Invalid residual type'); 
end

if (nargin < 2)
  % Use default residual ranges for code and phase, if
  % no value is entered. 
  ymax = 0;
end

%% Load appropriate data sets
close all;
loadCfg;
if exist([directoryname 'gielog_switch.txt'], 'file')
  switchFileExists = 1;
  switchLoop = load([directoryname 'gielog_switch.txt']);
  % switchLoop contains
  % time in ms, svid, state
  % state = 16 <=> in Loop B switching to Loop L
  % state = 17 <=> in Loop L switching to Loop B
  % state = 255 <=> Loss of Lock in Loop L
  loopBtoL = 16; 
  loopLtoB = 17; 
  loopLlost = 255;
  
  % adjust Time to seconds, and round to closest Ts.
  switchLoop(:,1) = switchLoop(:,1)/1000;
  switchLoopTime = Ts*round(switchLoop(:,1)/Ts);
  
  % kludge, to account for repeated timestamps
  tmpidx = 1; % set to nonempty to perform while at least once
  while(~isempty(tmpidx))
    diffSwLT = diff(switchLoopTime);
    tmpidx = find(diffSwLT == 0);
    switchLoopTime(tmpidx+1) = switchLoopTime(tmpidx+1) + Ts;
  end
  
else
  switchFileExists = 0; 
end

if exist([directoryname 'gielog_freq.txt'], 'file')
  freqFileExists = 1;
  freqErrs = load([directoryname 'gielog_freq.txt']);
  % switchLoop contains
  % time in ms, svid, freq err, other stuff
  
  % adjust Time to seconds
  freqErrs(:,1) = freqErrs(:,1)/1000;
else
  freqFileExists = 0; 
end


if (residType == PHASE)
  load([directoryname 'combinedDdprData']);
  ddResGie = ddprGie;
  ddResRov = ddprRov;
else
  load([directoryname 'combinedDdcrData']);
  ddResGie = ddcrGie;
  ddResRov = ddcrRov;
end

%% Disable warning
warnstruct = warning('query', 'Matlab:divideByZero');
disp('disable divide by zero warning'); 
warning('off', 'Matlab:divideByZero');


% pause on; 
pause off; 

if (useTimeOffset)
  timeOffset = startTime; 
else
  timeOffset = 0; 
end

startTime = startTime - timeOffset;
endTime = endTime - timeOffset;

% Time duration
hrs = floor(duration/3600);
mins = floor((duration - hrs*3600)/60);
xlimits = [startTime - duration*0.02, endTime + duration*0.02];

% xlimits = [5.202e5, 5.210e5];

if (residType == PHASE)
  % Set the default range of values that double differenced phase 
  % residuals could take.
  residRange = [-1.2,1.2]; % cycles
else
  % Set the default range of values that double differenced code 
  % residuals could take. 
  residRange = [-15,15]; % meters
end
if (ymax > 0)
  ylimits = [-ymax,ymax];
  residRange = [-ymax, ymax]; 
else
  ylimits = residRange;
end
  
% Starting location for text below legend. 
xText = xlimits(2)+duration*0.025;

% Create several figures - some for plotting everything together
% and some for plotting each SV separately.
labelFileNames = strrep(fileNames, '_', '\_');

% Underscores are a command string indicating
% subscript. Replace with real underscores for the labels.
if (plotFigC)
  % combined figure
  figC = newfigure;
  hold on;
  grid minor; 
  suptitlebar(['Double Diff ' residTypeName ' Residuals by SV']);
end
if (plotFigC2)
  % second combined figure
  figC2 = newfigure;
  subplot(2,1,1);
  grid on; 
  hold on;
  subplot(2,1,2);
  grid on;
  hold on;
end
if (plotFigS)
  % separate figure
  figS = newfigure;
end
if (plotFigS2)
  % separate figure
  figS2 = newfigure;
  legStrS2 = {'common' 'GIE only' 'Rover only'};
end
if (printStats)
  fid = fopen([directoryname residTypeStr 'Stats.txt'],'w');

  fprintf(fid, '%% %s\r\n', [directoryname residTypeStr 'Stats.txt']);
  fprintf(fid, '%% svid, epochs, mu, sigma, epochs, mu, sigma, type %% type string \r\n'); 
  fprintf(fid, '%%                           type: all = 0, only = 1, common = 2 \r\n');
  fprintf(fid, '%%\r\n');
  fprintf(fid, '%% \t GIE \t\t\t\t ROVER \r\n');

end

% Create a legend cellstring
legStr = {'Rover' 'Gie' 'only Rov' 'only Gie'};

% Set up color list
colour_LUT = setColourCodes; 
legidx = 0;

% Specify bins for histogram and initialize 
% histogram storage
binWidth = (residRange(2) - residRange(1))/40; 
bins = residRange(1)+binWidth/2:binWidth:residRange(2)-binWidth/2;

histGonly = zeros(size(bins));
histGcommon = histGonly; 
histRonly = histGonly;
histRcommon = histGonly; 

% Initialize statistics
MAX_SVS = 32; 
epochsG = zeros(MAX_SVS,1);
meanG = zeros(MAX_SVS,1);
sigmaG = zeros(MAX_SVS,1);
epochsGonly = zeros(MAX_SVS,1);
meanGonly = zeros(MAX_SVS,1);
sigmaGonly = zeros(MAX_SVS,1);
epochsGcommon = zeros(MAX_SVS,1);
meanGcommon = zeros(MAX_SVS,1);
sigmaGcommon = zeros(MAX_SVS,1);

epochsR = zeros(MAX_SVS,1);
meanR = zeros(MAX_SVS,1);
sigmaR = zeros(MAX_SVS,1);
epochsRonly = zeros(MAX_SVS,1);
meanRonly = zeros(MAX_SVS,1);
sigmaRonly = zeros(MAX_SVS,1);
epochsRcommon = zeros(MAX_SVS,1);
meanRcommon = zeros(MAX_SVS,1);
sigmaRcommon = zeros(MAX_SVS,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup Complete
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start generating data for combined plots and 
% actual separate plots. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Loop through the satellites
if (exist('chooseSv', 'var') == 1)
  % If we have chosen an SV, only plot that SV
  startSv = chooseSv; 
  endSv = chooseSv;
else
  startSv = 1; 
  endSv = 32;
end

for sv = startSv:endSv
        
  % Set up Rover data
  seg = find(ddResRov(:,rkey.SV) == sv);
  timeR = ddResRov(seg,rkey.TIME) - timeOffset*1000;
  resR = ddResRov(seg,rkey.RESID);
  slipR = ddResRov(seg,rkey.SLIP);

  % Set up GIE data
  seg = find(ddResGie(:,rkey.SV) == sv);
  timeG = ddResGie(seg,rkey.TIME) - timeOffset*1000;
  resG = ddResGie(seg,rkey.RESID);
  slipG = ddResGie(seg,rkey.SLIP);

  % Classify epochs into common and uncommon segments.
  [segRonly, segGonly, segRcommon, ...
    segGcommon, skipRov, skipGie] = ...
    classifyEpochs (timeR, timeG, startTime, endTime, Ts);

  if (skipGie && skipRov)
    continue; 
  end
  
  %% Set up full time duration vectors with data in appropriate spots.
  fullTime = startTime:Ts:endTime;

  if (~skipGie)
    % Initialize
    fullDataGcommon = fullTime*NaN;
    fullDataGonly = fullDataGcommon; 
    
    tIdx = (timeG - startTime*1000)/(Ts*1000);
    % Correct for matlab's 1 based numbering
    tIdx = tIdx + 1;
    
    fullDataGcommon(tIdx(segGcommon)) = resG(segGcommon);
    fullDataGonly(tIdx(segGonly)) = resG(segGonly); 
    
    fullDataG = fullDataGcommon; 
    fullDataG(tIdx(segGonly)) = resG(segGonly);    
  end

  if (~skipRov)
    % Initialize
    fullDataRcommon = fullTime*NaN;
    fullDataRonly = fullDataRcommon; 
    
    tIdx = (timeR - startTime*1000)/(Ts*1000);
    % Correct for matlab's 1 based numbering
    tIdx = tIdx + 1;
    
    fullDataRcommon(tIdx(segRcommon)) = resR(segRcommon);
    fullDataRonly(tIdx(segRonly)) = resR(segRonly); 
    
    fullDataR = fullDataRcommon; 
    fullDataR(tIdx(segRonly)) = resR(segRonly);
  end

  % Set up colors
  colorGcommon = 'b.';
  colorRcommon = 'r*';
  colorGonly = 'g+';
  colorRonly = 'm+';

  if ~(skipGie && skipRov)
    color = colour_LUT(sv,:);
    legidx = legidx + 1;
  end

  %%%%%%%%%%%%%%%%%%%%%%%
  %% Generate Statistics
  %%%%%%%%%%%%%%%%%%%%%%%
  
  if (plotFigS2 || plotFigH || printStats)
    epochsG(sv) = length(resG);
    meanG(sv) = mean(resG);
    sigmaG(sv) = std(resG);
    epochsGonly(sv) = length(segGonly);
    meanGonly(sv) = mean(resG(segGonly));
    sigmaGonly(sv) = std(resG(segGonly));
    epochsGcommon(sv) = length(segGcommon);
    meanGcommon(sv) = mean(resG(segGcommon));
    sigmaGcommon(sv) = std(resG(segGcommon));

    epochsR(sv) = length(resR);
    meanR(sv) = mean(resR);
    sigmaR(sv) = std(resR);
    epochsRonly(sv) = length(segRonly);
    meanRonly(sv) = mean(resR(segRonly));
    sigmaRonly(sv) = std(resR(segRonly));
    epochsRcommon(sv) = length(segRcommon);
    meanRcommon(sv) = mean(resR(segRcommon));
    sigmaRcommon(sv) = std(resR(segRcommon));
  end
  
  if (plotFigH)
    % Accumulate histogram information
    histGonly = histGonly + hist(resG(segGonly), bins);
    histGcommon = histGcommon + hist(resG(segGcommon), bins);

    histRonly = histRonly + hist(resR(segRonly), bins);
    histRcommon = histRcommon + hist(resR(segRcommon), bins);
  end
  
  %%%%%%%%%%%%%%%%
  %% CombinedFig
  %%%%%%%%%%%%%%%%
  if (plotFigC)
    figure(figC);
    
    if (residType == PHASE)
      residScale = 1.0;
    else % (residType == CODE)
      residScale = 0.1;
    end
    % Plot the measurement residuals
    if (~isempty(segRcommon))
      handleC(1) = plot(timeR(segRcommon)/1000,...
        sv+resR(segRcommon)*residScale,colorRcommon);
    end
    if (~isempty(segGcommon))
      handleC(2) = plot(timeG(segGcommon)/1000,...
        sv+resG(segGcommon)*residScale,colorGcommon);
    end
    if (~isempty(segRonly))
      handleC(3) = plot(timeR(segRonly)/1000,...
        sv+resR(segRonly)*residScale,colorRonly);
    end
    if (~isempty(segGonly))
      handleC(4) = plot(timeG(segGonly)/1000,...
        sv+resG(segGonly)*residScale,colorGonly);
    end
  end % (plotfigC)

  %%%%%%%%%%%%%%%%
  %% CombinedFig2
  %%%%%%%%%%%%%%%%
  % Plot two subplots -- one of all the GIE epochs, and
  % one of all the Rover epochs, of all SVs.

  if (plotFigC2)
    figure(figC2);

    subplot(2,1,1);
    hold on;
    if (~skipGie)
      handleC2G(legidx) = plot(fullTime, fullDataG, 'Color', color);
    elseif (~skipRov)
      handleC2G(legidx) = plot(startTime, NaN, 'Color', color);
    end

    subplot(2,1,2);
    hold on;
    if (~skipRov)
      handleC2R(legidx) = plot(fullTime, fullDataR, 'Color', color);
    elseif (~skipGie)
      handleC2R(legidx) = plot(startTime, NaN, 'Color', color);
    end

    legStrC2{legidx} = sprintf('SV %d', sv);
  end % (plotfigC2)

  %%%%%%%%%%%%%%%%
  %% SeparateFig
  %%%%%%%%%%%%%%%%

  if (plotFigS)
    % Set up separate figure
    figure(figS);
    clf(figS);
    ylabel([residTypeName ' Residuals']);
    xlabel('Time');
    xlim(xlimits);
    ylim(ylimits);
    hold on;

    % Plot the measurement residuals
    if (~isempty(segRcommon))
      handleS(1) = plot(timeR(segRcommon)/1000,...
        resR(segRcommon),colorRcommon);
    end
    if (~isempty(segGcommon))
      handleS(2) = plot(timeG(segGcommon)/1000,...
        resG(segGcommon),colorGcommon);
    end
    if (~isempty(segRonly))
      handleS(3) = plot(timeR(segRonly)/1000,...
        resR(segRonly),colorRonly);
    end
    if (~isempty(segGonly))
      handleS(4) = plot(timeG(segGonly)/1000,...
        resG(segGonly),colorGonly);
    end
    titlestr = sprintf('Double Diff %s Residuals for SV %d', ...
      residTypeName, sv);
    suptitlebar(titlestr);
    legend(handleS, legStr, -1);

    if (skipGie && skipRov)
      % Print a message saying no data.
      hold off;
      plot(1,1, 'w');
      %% Set up the plots to be cleared
      text(0.8,1,'No Data', 'FontSize', 16);
    end

    %% Print the separate png file
    % Append a leading '0' if the SV id is less than 10.
    if (printPlots && ~(skipGie && skipRov))
      if (sv < 10)
        pngFileName = [resultsdir 'dd' residTypeName 'Res-old-SV0' num2str(sv)];
      else
        pngFileName = [resultsdir 'dd' residTypeName 'Res-old-SV' num2str(sv)];
      end
      print( '-dpng', pngFileName);
    end

  end % (plotFigS)

  
  %%%%%%%%%%%%%%%%
  %% SeparateFig2
  %%%%%%%%%%%%%%%%
  % Plot two subplots -- one of all the GIE epochs, and
  % one of all the Rover epochs, of individual SVs.

  if (plotFigS2)
    figure(figS2);
    set(figS2, 'Position', [1 1 1280 957]);

    clf; 

    h1 = subplot(2,1,1);
    hold on;
    if (~skipGie)
      handleS2(1) = plot(fullTime, fullDataGcommon, 'g');
      handleS2(2) = plot(fullTime, fullDataGonly, 'b');
      handleS2(3) = plot(startTime, NaN, 'r');
    elseif (~skipRov)
      handleS2(1) = plot(startTime, NaN, 'g');
      handleS2(2) = plot(startTime, NaN, 'b');      
      handleS2(3) = plot(startTime, NaN, 'r');
    end
    xlim(xlimits);
    ylim(ylimits);
    grid minor; 
    titlestr = ...
      sprintf('GIE (%s) All: \\mu=%.3f \\sigma=%.3f epochs=%d', ...
      labelFileNames{GIE}, meanG(sv), sigmaG(sv), epochsG(sv));
    title(titlestr);


    h2 = subplot(2,1,2);
    hold on;
    if (~skipRov)
      plot(fullTime, fullDataRcommon, 'g');
      plot(fullTime, fullDataRonly, 'r');
    elseif (~skipGie)
      plot(startTime, NaN, 'g');
      plot(startTime, NaN, 'r');      
    end

    if (~(skipRov && skipGie))
      titlestr = sprintf('Double Diff %s Residuals for SV %d', ...
        residTypeName, sv);
      suptitle(titlestr);
      xlim(xlimits);
      ylim(ylimits)
      grid minor;
      titlestr = ...
        sprintf('Rover (%s) All: \\mu=%.3f \\sigma=%.3f epochs=%d', ...
        labelFileNames{ROVER}, meanR(sv), sigmaR(sv), epochsR(sv));
      title(titlestr);
      xlabelstr = sprintf( ...
        'Time[GPS Sec] %02dhrs:%02dmins - Base (%s)', ...
        hrs, mins, labelFileNames{BASE});
      xlabel(xlabelstr);
      if (residType == PHASE)
        ylabel('Residuals - cycles');
      else
        ylabel('Residuals - meters'); 
      end

      legend(handleS2, legStrS2, -1);
      p1 = get(h1, 'Position');
      p2 = get(h2, 'Position');
      p2(3) = p1(3);
      set(h2, 'Position', p2);
      
      %% Add text statistics
      % Base text on coordinates of h1.
      subplot(h1);
      
      % Gie Only Epochs
      txtstr = 'GIE only'; 
      text(xText, -0.3*ylimits(2), txtstr);
      
      txtstr = ...
        sprintf('\\mu=%.3f. \\sigma=%.3f.', ...
        meanGonly(sv), sigmaGonly(sv));
      text(xText, -0.6*ylimits(2), txtstr);
      
      txtstr = ...
        sprintf('epochs=%d', epochsGonly(sv));
      text(xText, -0.9*ylimits(2), txtstr);

      % Gie Common Epochs
      txtstr = 'GIE common'; 
      text(xText, -1.3*ylimits(2), txtstr);

      txtstr = ...
        sprintf('\\mu=%.3f. \\sigma=%.3f.', ...
        meanGcommon(sv), sigmaGcommon(sv));
      text(xText, -1.6*ylimits(2), txtstr);

      txtstr = ...
        sprintf('epochs=%d', epochsGcommon(sv));
      text(xText, -1.9*ylimits(2), txtstr); 

      %       % Gie Only Epochs
      %       txtstr = 'GIE only';
      %       text(xText, 2.3*ylimits(2), txtstr);
      %
      %       txtstr = ...
      %         sprintf('\\mu=%.3f. \\sigma=%.3f.', ...
      %         meanGonly(sv), sigmaGonly(sv));
      %       text(xText, 2.0*ylimits(2), txtstr);
      %
      %       txtstr = ...
      %         sprintf('epochs=%d', epochsGonly(sv));
      %       text(xText, 1.7*ylimits(2), txtstr);
      %
      %       % Gie Common Epochs
      %       txtstr = 'GIE common';
      %       text(xText, 1.3*ylimits(2), txtstr);
      %
      %       txtstr = ...
      %         sprintf('\\mu=%.3f. \\sigma=%.3f.', ...
      %         meanGcommon(sv), sigmaGcommon(sv));
      %       text(xText, 1*ylimits(2), txtstr);
      %
      %       txtstr = ...
      %         sprintf('epochs=%d', epochsGcommon(sv));
      %       text(xText, 0.7*ylimits(2), txtstr);
      %

      % Base text on coordinates of h2.
      subplot(h2);

      % Rover Only Epochs
      txtstr = 'Rover only'; 
      text(xText, 0.3*ylimits(2), txtstr);
      
      txtstr = ...
        sprintf('\\mu=%.3f. \\sigma=%.3f.', ...
        meanRonly(sv), sigmaRonly(sv));
      text(xText, 0*ylimits(2), txtstr);

      txtstr = ...
        sprintf('epochs=%d', epochsRonly(sv));
      text(xText, -0.3*ylimits(2), txtstr); 

      % Rover Common Epochs      
      txtstr = 'Rover common'; 
      text(xText, -0.7*ylimits(2), txtstr);

      txtstr = ...
        sprintf('\\mu=%.3f. \\sigma=%.3f.', ...
        meanRcommon(sv), sigmaRcommon(sv));
      text(xText, -1*ylimits(2), txtstr);

      txtstr = ...
        sprintf('epochs=%d', epochsRcommon(sv));
      text(xText, -1.3*ylimits(2), txtstr); 
      
    end

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Vivek's additions for Loop L
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    subplot(h1);
    for j = 1:2
      if (j==2)
        subplot(h2);
      end

      if (switchFileExists)
      
        oldmethod = 0;
        
        if (oldmethod)
          seg = find(switchLoop(:,2) == sv);

          segidx = find(switchLoop(seg,3) == loopBtoL);
          stem(switchLoop(seg(segidx),1),-ymax*ones(size(segidx)),'k','filled');

          segidx = find(switchLoop(seg,3) == loopLtoB);
          stem(switchLoop(seg(segidx),1),+ymax*ones(size(segidx)),'k','filled');

          segidx = find(switchLoop(seg,3) == loopLlost);
          stem(switchLoop(seg(segidx),1),+ymax*ones(size(segidx)),'k','filled');
        else % test
          % Test out method to plot state of Loop L.
          
          loopState = zeros(size(fullTime));
          % First select all the switches that occurred for this SV
          svSelect = find(switchLoop(:,2) == sv);
          % Now select out the time stamps for these switches and put
          % them in an array switchTime
          switchTime = switchLoopTime(svSelect);
          
          % Find the indices into fullTime which correspond to the
          % timestamps at which switches occurred. 
          [comT, idxSw, idxFull] = intersect(switchTime, fullTime);
          
          
          % Select the various types of switches, and indicate them
          % at the appropriate times within loopState. 
          seg = find(switchLoop(svSelect,3) == loopBtoL);
          loopState(idxFull(seg)) = 1;
          
          seg = find(switchLoop(svSelect,3) == loopLtoB);
          loopState(idxFull(seg)) = -1;

          seg = find(switchLoop(svSelect,3) == loopLlost);
          loopState(idxFull(seg)) = -1;
          
          loopState = cumsum(loopState);
          
          plot(fullTime, loopState*ymax/20, 'k');
          
        end
          
      else % (switchFileExists)
        for k = 1:length(tStartLoopL)
          if (sv == svLoopL(k))
            stem(tStartLoopL(k), -ymax, 'k', 'filled');
            stem(tEndLoopL(k), +ymax, 'k', 'filled');
          end
        end

      end % (switchFileExists)

      if (freqFileExists)

        seg = find(freqErrs(:,2) == sv);
        % stem(freqErrs(seg,1),+ymax*ones(size(seg)),'m.');
        plot(freqErrs(seg,1),+ymax/2*ones(size(seg)),'m.');        
        
      else % (freqFileExists)
      
        for k = 1:length(tEvent)
          if (sv == svEvent(k))
            stem(tEvent(k), ymax, 'rd', 'filled');
            text(tEvent(k), ymax, descEvent{k});
          end
        end

      end % (freqFileExists)
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % End of Vivek's additions for Loop L
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % If only plotting a few figures, link the axes.
    if (enableLinkAxes)
      hAll = [h1, h2];
      linkaxes(hAll, 'x');
      zoom on; 
    end
    
    %% Print the separate png file
    % Append a leading '0' if the SV id is less than 10.
    if (printPlots && ~(skipGie && skipRov))
      if (sv < 10)
        pngFileName = [resultsdir 'dd' residTypeName 'Res-SV0' num2str(sv)];
      else
        pngFileName = [resultsdir 'dd' residTypeName 'Res-SV' num2str(sv)];
      end
      saveas(gcf, [pngFileName], 'bmp');
    end   
    
  end % (plotFigS2)
  
  if (printStats)
    %%       GIE                ROVER
    %% svid, epochs, mu, sigma, epochs, mu, sigma, type %% type 
    %% (all = 0, only = 1, common = 2)

    fprintf(fid, '%2d % 8d % 2.3f % 2.3f % 8d % 2.3f % 2.3f %1d %% %s\r\n', ...
      sv, ...
      epochsG(sv), meanG(sv), sigmaG(sv), ...
      epochsR(sv), meanR(sv), sigmaR(sv), ...
      0, 'all');
      
    fprintf(fid, '%2d % 8d % 2.3f % 2.3f % 8d % 2.3f % 2.3f %1d %% %s\r\n', ...
      sv, ...
      epochsGonly(sv), meanGonly(sv), sigmaGonly(sv), ...
      epochsRonly(sv), meanRonly(sv), sigmaRonly(sv), ...
      1, 'only');
    
    fprintf(fid, '%2d % 8d % 2.3f % 2.3f % 8d % 2.3f % 2.3f %1d %% %s\r\n', ...
      sv, ...
      epochsGcommon(sv), meanGcommon(sv), sigmaGcommon(sv), ...
      epochsRcommon(sv), meanRcommon(sv), sigmaRcommon(sv), ...
      1, 'common');

  end

  if ((plotFigS || plotFigS2) && ~printPlots)
    pause;
  end

end % sv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of loop through the satellites
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Perform combined figure specific tasks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (plotFigC)
  figure(figC);
  if (residType == PHASE)
    ylabel('Phase Residuals - svid + cycles');
  else
    ylabel('Code Residuals - svid + meters/10');
  end
  xlabelstr = sprintf( ...
    'Time[GPS Sec] %02dhrs:%02dmins', hrs, mins);
  xlabel(xlabelstr);
  xlim(xlimits);
  ylim([0, 33]);
    
  txtstr = sprintf('GIE:%s', labelFileNames{GIE});   
  text(xText,20,txtstr);

  txtstr = sprintf('Rover:%s', labelFileNames{ROVER});
  text(xText,16,txtstr);

  txtstr = sprintf('Base:%s', labelFileNames{BASE});
  text(xText,12,txtstr);

  
  legend(handleC, legStr, -1);
  suptitlebar(['Double Diff ' residTypeName ' Residuals by SV']);
  %% Print the combined figure
  if (printPlots)
    pngFileName = [resultsdir 'dd' residTypeName 'Res-old'];
    print( '-dpng', pngFileName);
  end
  
end % plotFigC

if (plotFigC2)
  figure(figC2);

  h1 = subplot(2,1,1);
  legend(handleC2G, legStrC2, -1);
  titlestr = sprintf('GIE (%s) All Epochs',labelFileNames{GIE});
  title(titlestr);
  xlim(xlimits);
  ylim(ylimits);

  h2 = subplot(2,1,2);
  titlestr = sprintf('Rover (%s) All Epochs',labelFileNames{ROVER});
  title(titlestr);
  xlabelstr = sprintf( ...
    'Time[GPS Sec] %02dhrs:%02dmins - Base (%s)', ...
    hrs, mins, labelFileNames{BASE});
  xlabel(xlabelstr);
  if (residType == PHASE)
    ylabel('Phase Residuals - cycles');
  else
    ylabel('Code Residuals - meters'); 
  end
  xlim(xlimits);
  ylim(ylimits);
  
  p1 = get(h1, 'Position');
  p2 = get(h2, 'Position');
  p2(3) = p1(3);
  set(h2, 'Position', p2);
  %% Print the combined figure

  suptitlebar(['Double Diff ' residTypeName ' Residuals for All SVs']);
  % Suptitlebar changed the position of h1. Therefore 
  % change it back. 
  p1 = get(h1, 'Position');
  p1(3) = p2(3);
  set(h1, 'Position', p1);
  

  if (printPlots)
    pngFileName = [resultsdir 'dd' residTypeName 'Res'];
    print( '-dpng', pngFileName);
  end

end % plotFigC2


% Compute composite(total) statistics
% Each variable is labelled with a trailing 'T' which stands for
% "Total" statistics
[meanGT, sigmaGT, epochsGT] = ...
  computeCompositeStats(meanG, sigmaG, epochsG);

[meanGonlyT, sigmaGonlyT, epochsGonlyT] = ...
  computeCompositeStats(meanGonly, sigmaGonly, epochsGonly);

[meanGcommonT, sigmaGcommonT, epochsGcommonT] = ...
  computeCompositeStats(meanGcommon, sigmaGcommon, epochsGcommon);

[meanRT, sigmaRT, epochsRT] = ...
  computeCompositeStats(meanR, sigmaR, epochsR);

[meanRonlyT, sigmaRonlyT, epochsRonlyT] = ...
  computeCompositeStats(meanRonly, sigmaRonly, epochsRonly);

[meanRcommonT, sigmaRcommonT, epochsRcommonT] = ...
  computeCompositeStats(meanRcommon, sigmaRcommon, epochsRcommon);
  

% Compute total histogram values
histG = histGonly+histGcommon;
histR = histRonly+histRcommon;


% Normalize histogram information
histG = histG/epochsGT;
histR = histR/epochsRT;
histGonly = histGonly/epochsGonlyT;
histRonly = histRonly/epochsRonlyT;
histGcommon = histGcommon/epochsGcommonT;
histRcommon = histRcommon/epochsRcommonT;
histMax = max([histG, histR, histGonly, histRonly, ...
  histGcommon, histRcommon ]);


if (plotFigH)  
  % histogram
  figH = newfigure;
  
  h1 = histPlot([3,2,1], bins, histG, residRange, 'GIE all', ...
    meanGT, sigmaGT, epochsGT, histMax);
  
  h2 = histPlot([3,2,2], bins, histR, residRange, 'Rover all', ...
    meanRT, sigmaRT, epochsRT, histMax);

  h3 = histPlot([3,2,3], bins, histGcommon, residRange, 'GIE common', ...
    meanGcommonT, sigmaGcommonT, epochsGcommonT, histMax);
  
  h4 = histPlot([3,2,4], bins, histRcommon, residRange, 'Rover common', ...
    meanRcommonT, sigmaRcommonT, epochsRcommonT, histMax);

  h5 = histPlot([3,2,5], bins, histGonly, residRange, 'GIE only', ...
    meanGonlyT, sigmaGonlyT, epochsGonlyT, histMax);
  if (residType == PHASE)
    xlabelstr = sprintf('L1 cycles. GIE - %s. %02dhrs:%02dmins', ...
      labelFileNames{GIE}, hrs, mins);
  else % (residType == CODE)
    xlabelstr = sprintf('meters. GIE - %s. %02dhrs:%02dmins', ...
      labelFileNames{GIE}, hrs, mins);
  end
  xlabel(xlabelstr);

  % This order was used because suptitlebar deletes one of 
  % the non-standard axes. 
  h6 = subplot(3,2,6);
  suptitlebar(['Double Difference ' residTypeName ' Residual Histograms']);

  histPlot([3,2,6, h6], bins, histRonly, residRange, 'Rover only', ...
    meanRonlyT, sigmaRonlyT, epochsRonlyT, histMax);
  xlabelstr = sprintf('Rover - %s. Base - %s.', ...
    labelFileNames{ROVER}, labelFileNames{BASE});
  xlabel(xlabelstr);

  if (printPlots)
    pngFileName = [resultsdir 'dd' residTypeName 'Res-hist'];
    print( '-dpng', pngFileName);
  end
  
end % plotFigH

if (printStats)
  %%       GIE                ROVER
  %% svid, epochs, mu, sigma, epochs, mu, sigma, type %% type
  %% (all = 0, only = 1, common = 2)

  fprintf(fid, '%2d % 8d % 2.3f % 2.3f % 8d % 2.3f % 2.3f %1d %% %s\r\n', ...
    0, ...
    epochsGT, meanGT, sigmaGT, ...
    epochsRT, meanRT, sigmaRT, ...
    0, 'all, all SVs');

  fprintf(fid, '%2d % 8d % 2.3f % 2.3f % 8d % 2.3f % 2.3f %1d %% %s\r\n', ...
    0, ...
    epochsGonlyT, meanGonlyT, sigmaGonlyT, ...
    epochsRonlyT, meanRonlyT, sigmaRonlyT, ...
    1, 'only, all SVs');

  fprintf(fid, '%2d % 8d % 2.3f % 2.3f % 8d % 2.3f % 2.3f %1d %% %s\r\n', ...
    0, ...
    epochsGcommonT, meanGcommonT, sigmaGcommonT, ...
    epochsRcommonT, meanRcommonT, sigmaRcommonT, ...
    1, 'common, all SVs');

  fclose(fid);
end

%% Reset warning
warning(warnstruct);
disp('reset divide by zero warning'); 

return;
