search.m

it is a source code from gps toolbox

function [ncan, cands, disall, ipos] = search(Chic, MaxCan, n, a, D, L)
% SEARCH   finds 'MaxCan' integer vectors whose distances to the real
%          vector 'a' are minimal in the metric of Q = transpose(L) D L.
%          Only integer vectors with a distance less than sqrt(Chic)
%          are regarded.
%
%          The search for gridpoints inside the ambiguity search ellipsoid
%          is a sequential conditional adjustment upon the ambiguities.
%          The search starts by conditioning the last ambiguity a_n to
%          an integer, then the one-but-last a_{n-1} etc., until
%             - the squared norm grows too large (out of the ellipsoid)
%             - an integer for a_1 is found: a full integer vector is
%               encountered (a gridpoint inside the ellipsoid)
%          then it goes back to some previous (towards a_n) ambiguity
%          and considers another integer for it.

% Chic     Chi squared
% MaxCan   number of minimum integer vectors required
% n        dimension of matrix
% a        vector of real valued estimates \hat{a} (float
%                        solution)
% D        diagonal matrix
% L        lower triangular matrix: Q^{-1} = L D L^T
%          although L is lower triangular, in this example
%          program L is stored column-wise in a 2-dimensional
%          array, to avoid the necessity of a dedicated
%          storage scheme.
% lef      vector of length  n
% left     vector of length  n+1
% right    vector of length  n+1
% dist     difference between the integer tried and \hat{a}_i,
%          length n
% endd     vector of length n
% dq       vector of length n
% ncan     number of integer vectors found
% cands    | 2-dimensional array to store the candidates
% disall   | according squared norms \hat{a}-\check{a}
% ipos     column number in 'cands' where the candidate
%          belonging to the minimum distance is stored

% Delft Geodetic Computing Centre/LGR, Paul de Jonge
% copyright by Delft University of Technology, Faculty of Geodesy

% Recoded into MATLAB by Kai Borre 12-04-96

disall = zeros(1,MaxCan);
cands = zeros(n,MaxCan);
tmax = 0;
imax = 0;

if MaxCan < 1
   error('ERROR in SEARCH: number of requested candidates < 1')
else
   if n < 2
      error('ERROR in SEARCH: dimension of system < 2')
   end
   
   ende = 'false';
   right = zeros(n+1,1);
   left = zeros(n+1,1);
   right(n+1,1) = Chic;
   for i = 1:n-1
      dq(i) = D(i+1)/D(i);
   end
   dq(n) = 1/D(n);
   
   ncan = 0;
   i = n+1;
   iold = i;
   
   while ende == 'false'
      i = i-1;                             % go a level deeper
      if iold <= i                         % we were here before
         lef(i) = lef(i) + L(i+1,i);       % only one dist is one bigger
      else
         lef(i) = 0;
         for j = i+1:n
            lef(i) = lef(i) + L(j,i)*dist(j,1);
         end
      end
      iold = i;                     % keep track of level
      
      right(i) = (right(i+1)-left(i+1)) * dq(i);
      reach = sqrt(right(i));      
      % delta=a(i)-reach-lef(i) is the left border
      % ceil (delta) is the left most integer in the interval
      % dist(i,1) is the distance of this left most integer to the a_hat      
      delta = a(i) - reach - lef(i);
      dist(i,1) = ceil(delta) - a(i);      
      if dist(i,1) > reach-lef(i)   % there is nothing at this level
                                    % so, ... back track
         [i, dist, left, ende] = backtrac(n, i, endd, dist, lef, left, ende);
      else
         endd(i) = reach - lef(i) - 1;     % set the right 'border'
         left(i) = (dist(i,1)+lef(i))^2;
      end      
      if i == 1
         [dist, endd, ncan, cands, tmax, imax, disall] = ...
            collects(n, MaxCan, D, lef, left, right, ...
            Chic, dist, endd, ncan, disall, cands, tmax, imax);
         [i, dist, left, ende] = backtrac(n, i, endd, dist, lef, left, ende);
      end
   end
   % the candidate vectors are stored and the index for the best
   % candidate is set (ipos)
   dminn = 1.e20;                 % some big number
   for i = 1:min(ncan, MaxCan)     
      if disall(i) < dminn
         dminn = disall(i);
         ipos = i;
      end
      cands(:,i) = cands(:,i) + a;
   end
end
% cands contains integer transformed reduced ambiguities
%%%%%%%%%%%%%%% end search.m  %%%%%%%%%%%%%%%%%%%%%%%%%%%