function Processed = do_Xcorr(phaseDDiffL1, phaseDDiffL2, ...
                              rangeDDiffL1, rangeDDiffL2, ...
                              SavePlot, MinDataSize);

% *********************************************************************
% $Id: do_Xcorr.m,v 1.9 2005/04/18 19:46:46 vnadkar Exp $
%
% do_Xcorr
%
% Inputs:
%   phaseDDiffL1 - L1 Double Difference Phase residuals
%   phaseDDiffL2 - L2 Double Difference Phase residuals
%   rangeDDiffL1 - L1 Double Difference Range residuals
%   rangeDDiffL2 - L2 Double Difference range residuals
%   SavePlot     - 1 = Save .png files
%   MinDataSize  - Minimum contiguous data span before an ACF  is
%                  calculated
%
% Usage: 
% 
% Takes the output from piline  and looks for the longest period of
% double difference which is near continuous. Because of bugs in TRS 
% the occasional epoch of data is lost, this script attempts to 
% account for this problem. The longest set of continuous data must 
% be over "MinDatSize" seconds. If it is an  auto-correlation of the 
% data will be formed and saved to a file. An array is created and 
% returned to the calling function, the array has a "1" for each PRN
% that a ACF was created.
% 
% Requires:
%  xcorr.m
%
% 8/31/00 - modified to work as part of "load_pi"
%
% *********************************************************************

% nargin = 4;

FALSE = 0;
TRUE  = 1;
if(nargin < 4)
  error('do_Xcorr requires at least 4 input parameters');
elseif(nargin == 4)
  SavePlot = TRUE;
  MinDataSize = 1000;
elseif(nargin == 5)
  MinDataSize = 1000;
end

% The ACF gets inaccurate the bigger the time shift as less data is
% available to compute the function. I've tried, and failed, to get a
% good reference that discusses the accuracy of the ACF as a function of
% shift. So I'm limiting it to 10% of the minimum contiguous data set we
% are allowing to use used fo ran ACF computation.
MaxShift    = floor(MinDataSize / 10);

%
% Open files where we will store ACF results.
%
fidL1Phase  = fopen('ACFL1Phase.txt','w');
fidL2Phase  = fopen('ACFL2Phase.txt','w');
fidL1Range  = fopen('ACFL1Range.txt','w');
fidL2Range  = fopen('ACFL2Range.txt','w');

numPlottedRange = 0;
numPlottedPhase = 0;

k = [];

for obs=1:2
  if(obs == 1)
    dataL1 = rangeDDiffL1;
    dataL2 = rangeDDiffL2;
    TYPE = 'Pseudorange';
  else
    dataL1 = phaseDDiffL1;
    dataL2 = phaseDDiffL2;
    TYPE = 'Phase';
  end

  for SV=1:32
    MaxValidL1 = 0;
    MaxValidL2 = 0;

    % Find any NaNs that indicate there is no
    % data available
    i = find(isnan( dataL1(:,(SV+1))) == 1);
    
    % Keep track of which satellites we have do for. This array gets
    % returned to the calling function and is used to setup the HTML
    % index page.
    %
    % Default to zero, gets set to 1 later if a valid plot is created.
    Processed(SV) = 0; 

    if(length(i) ~= 0)
      %
      % There are some NaN's so there is missing data. Look for the
      % longest contiguous section of data.
      %
	  clear k;
      k(1) = 0;
      k(2:(length(i)+1)) = i;
      k((length(i)+2))   = length(dataL1)+1;
    
      % Correct for small gaps in the data - TRS is notorious for
      % droping data packets. We do not want to prevent the ACF
      % calculation due to a TRS data collection problem.
      %
      % First find all data segments that are longer than 10 seconds.
      % "k" contains the indexs of all the NaNs (that indicate there is
      % no data available), so the difference between the indicies is
      % the time we do have good data.
      % 
      z = find(diff(k) > 10);

      % If we have at least two segments of data
      if(length(z) >= 2)

        % If there are any missing epochs we are going to set them to
        % the value that has least impact. Therefore set to the mean of
        % the data set.
        j = find(isfinite( dataL1(:,(SV+1))) == 1);
        mu = mean(dataL1(j,(SV+1)));
        
        % Check each segment of data. If there is good data either epoch
        % beside a NaN assume TRS dropped a packet and so set this epoch
        % to zero so it can be used in the ACF.
        for xx = 2:length(z)
          if(    isnan( dataL1(k(z(xx))-1, SV+1))==0 & ... 
                 isnan( dataL1(k(z(xx))+1, SV+1))==0)
            dataL1(k(z(xx)),(SV+1)) = mu;
          end
        end

        % Recompute the index of missing data epochs now we have
        % corrected for single epochs.
        i = find(isnan(dataL1(:,(SV+1))) == 1);
		clear k;
        k(1) = 0;
        k(2:(length(i)+1)) = i;
        k((length(i)+2))   = length(dataL1)+1;
      end

      % Find the maximum contigous data segment.
      MaxValidL1 = max(diff(k));
      index = find(diff(k) == MaxValidL1);
      % incase there are two segments with the same timespan use the
      % first.
      index = index(1) + 1;
      % Grab the data we want to process further.
      ValidData = dataL1((k(index - 1)) + 1:(k(index)-1),(SV+1));
      MaxValidL1 = MaxValidL1 - 1;
    else
      % There is data for the entire observation period so use it all
      % for the ACF calculation.
      MaxValidL1 = length(dataL1);
      ValidData  = dataL1(:,SV+1);
    end

    if(MaxValidL1 > MinDataSize)
      X1 = xcorr(ValidData - mean(ValidData),'coeff');
      figure('Visible','off');
      plot(0:MaxShift,X1( MaxValidL1: (MaxValidL1+MaxShift)),'b');
      hold
      
      %
      % We have some L1 data - now process the L2 data. See the L1 data
      % processing for comments.
      % 
      i = find(isnan(dataL2(:,(SV+1))) == 1);
      if(length(i) ~= 0)
		clear k;
        k(1) = 0;
        k(2:(length(i)+1)) = i;
        k((length(i)+2))   = length(dataL2)+1;
      
        % Correct for small gaps in the data
        z = find(diff(k) > 10);
        if(length(z) >= 2)
          % If there are any missing epochs we are going to set them to
          % the value that has least impact. Therefore set to the mean 
          % of the data set.
          j = find(isfinite( dataL2(:,(SV+1))) == 1);
          mu = mean(dataL2(j,(SV+1)));
          
          for xx = 2:length(z)
            if(    isnan( dataL2(k(z(xx))-1, SV+1))==0 & ... 
                   isnan( dataL2(k(z(xx))+1, SV+1))==0)
              dataL2(k(z(xx)),(SV+1)) = mu;
            end
          end
          i = find(isnan(dataL2(:,(SV+1))) == 1);
		  clear k;
          k(1) = 1;
          k(2:(length(i)+1)) = i;
          k((length(i)+2))   = length(dataL2)+1;
        end

  
        MaxValidL2 = max(diff(k));
        index = find(diff(k) == MaxValidL2);
        % incase there are two segments with the same timespan use the
        % first.
        index = index(1) + 1;
        % Grab the data we want to process further.
        ValidData = dataL2((k(index - 1)) + 1:(k(index)-1),(SV+1));
        MaxValidL2 = MaxValidL2 - 1;
      else 
        % There is data for the entire observation period so use it all
        % for the ACF calculation.
        MaxValidL2 = length(dataL2);
        ValidData  = dataL2(:,SV+1);
      end

      if(MaxValidL2 > MinDataSize)
        X2 = xcorr(ValidData - mean(ValidData),'coeff');
        plot(0:MaxShift,X2( MaxValidL2: (MaxValidL2+MaxShift)),'r');
        plot(0:MaxShift,X1( MaxValidL1: (MaxValidL1+MaxShift)),'bx');
        plot(0:MaxShift,X2( MaxValidL2: (MaxValidL2+MaxShift)),'rx');
        legend('L1','L2');

        % Save for later
        if(obs == 1)
          % Dump to a text file:
          %
          % Format:
          % SV #Epoch e^-1 X(0) X(1) X(1).... X(299)
          %
          j = find( X1( MaxValidL1: length(X1)) > exp(-1));
          if(length(j) == 0)
            timeConst = NaN;
          else
            timeConst = length(j);
          end
          fprintf(fidL1Range,'%d %d %d ',SV,MaxValidL1,timeConst);
          
          j = find( X2( MaxValidL2: length(X2)) > exp(-1));
          if(length(j) == 0)
            timeConst = NaN;
          else
            timeConst = length(j);
          end
          fprintf(fidL2Range,'%d %d %d ',SV,MaxValidL2,timeConst);

          for j=0:299
            fprintf(fidL1Range,'%.3f ',X1(MaxValidL1 + j));
            fprintf(fidL2Range,'%.3f ',X2(MaxValidL2 + j));
          end
          fprintf(fidL1Range,'\n');
          fprintf(fidL2Range,'\n');

          L2ACFDataRange(numPlottedRange+1,:) = ...
                            X2( MaxValidL2: (MaxValidL2+MaxShift))';
          L1ACFDataRange(numPlottedRange+1,:) = ...
                            X1( MaxValidL1: (MaxValidL1+MaxShift))';
        else
          % Dump to a text file:
          %
          % Format:
          % SV #Epoch e^-1 X(0) X(1) X(1).... X(299)
          %
          j = find( X1( MaxValidL1: length(X1)) > exp(-1));
          if(length(j) == 0)
            timeConst = NaN;
          else
            timeConst = length(j);
          end
          fprintf(fidL1Phase,'%d %d %d ',SV,MaxValidL1,timeConst);
          
          j = find( X1( MaxValidL2: length(X2)) > exp(-1));
          if(length(j) == 0)
            timeConst = NaN;
          else
            timeConst = length(j);
          end
          fprintf(fidL2Phase,'%d %d %d ',SV,MaxValidL2,timeConst);

          for j=0:299
            fprintf(fidL1Phase,'%.3f ',X1(MaxValidL1 + j));
            fprintf(fidL2Phase,'%.3f ',X2(MaxValidL2 + j));
          end
          fprintf(fidL1Phase,'\n');
          fprintf(fidL2Phase,'\n');

          L2ACFDataPhase(numPlottedPhase+1,:) = ...
                            X2( MaxValidL2: (MaxValidL2+MaxShift))';
          L1ACFDataPhase(numPlottedPhase+1,:) = ...
                            X1( MaxValidL1: (MaxValidL1+MaxShift))';
        end
      else
        plot(0:MaxShift,X1( MaxValidL1: (MaxValidL1+MaxShift)),'b');
        legend('L1');
        % Save for later
        if(obs == 1)
          L2ACFDataRange(numPlottedRange+1,:) = ones(1,(MaxShift+1)).*NaN;
          L1ACFDataRange(numPlottedRange+1,:) =  ...
                            X1( MaxValidL1: MaxValidL1+MaxShift)';
        else
          L2ACFDataPhase(numPlottedPhase+1,:) = ones(1,(MaxShift+1)).*NaN;
          L1ACFDataPhase(numPlottedPhase+1,:) = ...
                            X1( MaxValidL1: MaxValidL1+MaxShift)';
        end
      end
      
	  A=axis;
      A(1) = -5;
      A(2) = MaxShift;
      axis(A);
      grid
      title([TYPE,' Auto-correlation PRN #',num2str(SV)]);
      xlabel('Time Shift [Sec]');
      ylabel('ACF');

      if(SavePlot == TRUE)
        print('-dpng',['ACF_',TYPE,'_',num2str(SV),'.PNG']);
        Processed(SV) = 1; % Valid plot has been created.
      end
      close
      %
      % Form an array which gives a list of the satellites we have
      % created plots for.
      %
      if(obs == 1)
        SVPlotted(numPlottedRange + 1) = SV;
        numPlottedRange = numPlottedRange + 1;
      else
        SVPlotted(numPlottedPhase + 1) = SV;
        numPlottedPhase = numPlottedPhase + 1;
      end
    end
  end 
end


%
% Now plot all the ACFs on the same graph
%

colour_LUT = setColourCodes;

figure('Visible', 'off');
subplot(2,2,1);
for SV = 1:numPlottedRange
  plot(0:MaxShift,L1ACFDataRange(SV,:),'color',colour_LUT(SV,:));
  if(SV == 1)
    hold;
  end
end
grid
title('L1 Pseudorange');
ylabel('ACF');
axis([-5 MaxShift -.2 1.1])

subplot(2,2,2);
for SV = 1:numPlottedRange
  plot(0:MaxShift,L2ACFDataRange(SV,:),'color',colour_LUT(SV,:));
  if(SV == 1)
    hold;
  end
end
grid
title('L2 Pseudorange');
ylabel('ACF');
axis([-5 MaxShift -.2 1.1])


subplot(2,2,3);
for SV = 1:numPlottedPhase
  plot(0:MaxShift,L1ACFDataPhase(SV,:),'color',colour_LUT(SV,:));
  if(SV == 1)
    hold;
  end
end
grid
title('L1 Carrier Phase');
ylabel('ACF');
xlabel('Time Shift [Sec]');
axis([-5 MaxShift -.2 1.1])


subplot(2,2,4);
for SV = 1:numPlottedPhase
  plot(0:MaxShift,L2ACFDataPhase(SV,:),'color',colour_LUT(SV,:));
  if(SV == 1)
    hold;
  end
end
grid
title('L2 Carrier Phase');
ylabel('ACF');
xlabel('Time Shift [Sec]');
axis([-5 MaxShift -.2 1.1])

if(SavePlot == TRUE)
  print -dpng ALL_ACF.png
end


fclose(fidL1Phase);
fclose(fidL2Phase);
fclose(fidL1Range);
fclose(fidL2Range);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% $Log: do_Xcorr.m,v $
% Revision 1.9  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.8  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.7  2001/09/19 18:19:25  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.6  2001/09/07 19:39:49  stuart
% o Resolved a bug with the input option
% o Minor code tidy up
%
% Revision 1.5  2000/11/27 22:01:56  stuart
% o Corrected bug
%
% Revision 1.4  2000/09/24 19:49:55  timo
% Remove highlight markers, >>>>> (changes between this and
% last rev). No functional changes.
%
% Revision 1.3  2000/09/22 18:51:13  timo
% Update to allow Matlab compilation (StandAlone application).
% Updates highlighed with >>>>> (to be removed).
%
% Revision 1.2  2000/08/21 13:57:33  stuart
% o Correct $Id
% o Add $Log
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

