function plotDdResVsIndepVar (residTypeStr, indepVarStr)
% function plotDdResVsIndepVar (residTypeStr, indepVarStr)
% Plot the statistics of double difference residuals versus
% an independent variable such as CN0, Elevation or Azimuth.
% 
% residTypeStr - Residual Type string is one of
%                'code', 'phase', with the default being
%                'phase'
% indepVarStr - Independent Variable String is one of 
%             'cn0', 'elv', 'az' with the default being 'cn0'
% 
% Vivek Nadkarni
% Trimble
% July 14, 2004


% Plotting Flags
printPlots = 1; % save plots to PNG
plotBar    = 0; % create a bar chart instead of a plot

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

% Select plotting vs CN0 or Elevation or Azimuth
CN0 = 1;
ELV = 2;
AZ = 3; 
if strcmp(indepVarStr, 'cn0')
  plotType = CN0; 
elseif strcmp(indepVarStr, 'elv')
  plotType = ELV;
elseif strcmp(indepVarStr, 'az')
  plotType = AZ; 
end  

% Load appropriate files
close all;
loadCfg;
if (residType == PHASE)
  load([directoryname 'combinedDdprData']);
  ddResGie = ddprGie; 
  ddResRov = ddprRov; 
else % (residType == CODE)
  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');

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

% Time duration
hrs = floor(duration/3600);
mins = floor((duration - hrs*3600)/60);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initialize Based on Plot Type
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (plotType == CN0)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Compute effective CN0, and then clean up
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  v1=10.^(-ddResRov(:,rkey.BaseRefSvSNR)*0.1);
  v2=10.^(-ddResRov(:,rkey.BaseTestSvSNR)*0.1);
  v3=10.^(-ddResRov(:,rkey.RovRefSvSNR)*0.1);
  v4=10.^(-ddResRov(:,rkey.RovTestSvSNR)*0.1);
  v=v1+v2+v3+v4;
  effCN0R = -10*log10(v);

  v1=10.^(-ddResGie(:,rkey.BaseRefSvSNR)*0.1);
  v2=10.^(-ddResGie(:,rkey.BaseTestSvSNR)*0.1);
  v3=10.^(-ddResGie(:,rkey.RovRefSvSNR)*0.1);
  v4=10.^(-ddResGie(:,rkey.RovTestSvSNR)*0.1);
  v=v1+v2+v3+v4;
  effCN0G = -10*log10(v);

  clear v1 v2 v3 v4 v;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % indepVar stands for independent variable against
  % which the double difference residual statistics
  % will be plotted. 
  indepVarG = effCN0G;
  indepVarR = effCN0R; 

  % Choose max and min CN0 in dB
  cn0Min = 0;
  cn0Max = 60;
  % Choose bin width in dB
  binWidth = 2; % dB of CN0

  %   nBins = (cn0Max-cn0Min)/binWidth;
  %   bins = (cn0Min + binWidth/2):binWidth:(cn0Max-binWidth/2);

  indepMin = cn0Min;
  indepMax = cn0Max; 
  
  indepVarName = 'Composite C/N_0'; 
  xLabelName = 'Composite C/N_0 in dB-Hz';
end

if (plotType == ELV)  
  % Use the elevation of the test SV
  % Can also use the average elevation of the test and
  % reference sv.
  elvR = ddResRov(:,rkey.TestSvElv);
  elvG = ddResGie(:,rkey.TestSvElv);
  indepVarG = elvG;
  indepVarR = elvR; 

  % Choose Max and Min elevations in degrees
  elvMin = 0;  % degrees
  elvMax = 90; % degrees
  % Choose bin width in degrees
  binWidth = 5; % degrees
  
  indepMin = elvMin;
  indepMax = elvMax;
  
  indepVarName = 'Test SV Elevation';
  xLabelName = 'Elevation in degrees';
end

if (plotType == AZ)
  % Use the azimuth of the test SV
  azR = ddResRov(:,rkey.TestSvAz);
  azG = ddResGie(:,rkey.TestSvAz);
  
  indepVarG = azG;
  indepVarR = azR;

  azMin = 0;  % degrees
  azMax = 360; % degrees
  % Choose bin width in degrees
  binWidth = 20; % degrees

  indepMin = azMin; 
  indepMax = azMax; 
  
  indepVarName = 'Test SV Azimuth';
  xLabelName = 'Azimuth in degrees';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% End of Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set up Bins
nBins = (indepMax-indepMin)/binWidth;
bins = (indepMin + binWidth/2):binWidth:(indepMax-binWidth/2);

%% Sift through bins to generate histogram
for bin = 1:nBins
  binLo = indepMin + (bin-1)*binWidth;
  binHi = indepMin + bin*binWidth;
  
  seg = find((binLo < indepVarG) & (indepVarG <= binHi));
  resG = ddResGie(seg,rkey.RESID);
  
  seg = find((binLo < indepVarR) & (indepVarR <= binHi));
  resR = ddResRov(seg,rkey.RESID);
  
  muG(bin) = mean(resG);
  sigG(bin) = std(resG);  
  nG(bin) = length(resG);

  muR(bin) = mean(resR);
  sigR(bin) = std(resR);
  nR(bin) = length(resR);  
end

%% Clean up data
seg = find(nG<=1);
muG(seg) = NaN;
sigG(seg) = NaN;

seg = find(nR<=1);
muR(seg) = NaN;
sigR(seg) = NaN;

goodBins = bins(...
  ~isnan(muG) | ~isnan(sigG) | ...
  ~isnan(muR) | ~isnan(sigR) );
xlimits = [goodBins(1) - binWidth/2, goodBins(end) + binWidth/2]; 

ylimitsMu = [ min([muG, muR]) , max([muG, muR]) ];
muAdjust = (ylimitsMu(2) - ylimitsMu(1))*0.1;
ylimitsMu = [ylimitsMu(1) - muAdjust, ylimitsMu(2) + muAdjust];
ylimitsSigma = [ 0 , max([sigG,sigR])*1.1 ];


%% Create plots

newfigure;

%%%%%%%%%%%%%%
%% Subplot 1
%%%%%%%%%%%%%%
subplot(2,2,1);
if (plotBar)
  bar(bins, muG);
else
  plot(bins,muG,'-x');
end
grid on;
xlim(xlimits);
ylim(ylimitsMu);
titlestr = sprintf('\\mu GIE (%d epochs)', sum(nG));
title(titlestr);
if (residType == PHASE)
  ylabel('L1 cycles');
else % (residType == CODE)
  ylabel('meters');
end

%%%%%%%%%%%%%%
%% Subplot 2
%%%%%%%%%%%%%%
subplot(2,2,2);
if (plotBar)
  bar(bins, muR);
else
  plot(bins,muR,'-x');
end
grid on;
xlim(xlimits);
ylim(ylimitsMu);
titlestr = sprintf('\\mu Rover (%d epochs)', sum(nR));
title(titlestr);
if (residType == PHASE)
  ylabel('L1 cycles');
else % (residType == CODE)
  ylabel('meters');
end

%%%%%%%%%%%%%%
%% Subplot 3
%%%%%%%%%%%%%%
subplot(2,2,3);
if (plotBar)
  bar(bins, sigG);
else
  plot(bins,sigG,'-x');
end
grid on;
xlim(xlimits);
ylim(ylimitsSigma);
titlestr = sprintf('\\sigma GIE - %s', labelFileNames{GIE});
title(titlestr);
% xlabel('Composite C/N_0 in dB-Hz');
xlabel(xLabelName); 
if (residType == PHASE)
  ylabel('L1 cycles');
else % (residType == CODE)
  ylabel('meters');
end

%%%%%%%%%%%%%%
%% Subplot 4
%%%%%%%%%%%%%%
subplot(2,2,4);
if (plotBar)
  bar(bins, sigR);
else
  plot(bins,sigR,'-x');
end
grid on;
xlim(xlimits);
ylim(ylimitsSigma);
titlestr = sprintf('\\sigma Rover - %s', labelFileNames{ROVER});
title(titlestr);
xlabelstr = sprintf( ...
  'Base - %s.  %02dhrs:%02dmins', ...
  labelFileNames{BASE},  hrs, mins);
xlabel(xlabelstr);
if (residType == PHASE)
  ylabel('L1 cycles');
else % (residType == CODE)
  ylabel('meters');
end

suptitlestr = sprintf('Double Diff %s Residual Stats vs. %s', ...
  residTypeName, indepVarName);
suptitlebar(suptitlestr); 

% Print to PNG file
if (printPlots)
  pngFileName = ...
    [resultsdir 'dd' residTypeName 'Res-' indepVarStr];
  print( '-dpng', pngFileName);
end

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

return; 
