%                             gdopv.m
%  Scope:   This MATLAB macro computes an approximate value of the geometric
%           dilution of precision (GDOP) when four line-of-sight unit vectors
%           are specified, by using an approximate geometric method [1].
%  Usage:   x = gdopv(g)
%  Description of parameters:
%           g  - input, line-of-sight unit vectors, each vector is
%                stored in another row (it should be 3 columns)
%           x  - output, approximate value of the geometric dilution of 
%                precision 
%  Reference:
%           [1] Massatt, P., Rudnick, K., Geometric formulas for dilution 
%               of precision calculations. Navigation, Journal of the Insti-
%               tute of Navigation, Vol. 37, No. 4, 1990-91, pp. 379-391. 
%  Last update:  06/15/00
%  Copyright (C) 1996-00 by LL Consulting. All Rights Reserved.

function  x = gdopv(g)

[m,n] = size(g);
if n ~= 3
   error('Error 1 -  GDOPV: see length of the unit line of sight vectors');
end
if m ~= 4
   error('Error 2 -  GDOPV: see number of the unit line of sight vectors');
end

eps1 = 1.e-10;
gdopmax = 1.e+10;

c(1) = 1. + g(1,1) * g(2,1) + g(1,2) * g(2,2) + g(1,3) * g(2,3);
c(2) = 1. + g(1,1) * g(3,1) + g(1,2) * g(3,2) + g(1,3) * g(3,3);
c(3) = 1. + g(1,1) * g(4,1) + g(1,2) * g(4,2) + g(1,3) * g(4,3);
c(4) = 1. + g(2,1) * g(3,1) + g(2,2) * g(3,2) + g(2,3) * g(3,3);
c(5) = 1. + g(2,1) * g(4,1) + g(2,2) * g(4,2) + g(2,3) * g(4,3);
c(6) = 1. + g(3,1) * g(4,1) + g(3,2) * g(4,2) + g(3,3) * g(4,3);

temp(1) = c(1) * c(6);
temp(2) = c(3) * c(4);
temp(3) = c(2) * c(5);

square = 4. * (c(1) * c(1) + c(2) * c(2) + c(3) * c(3) + c(4) * c(4)...
         + c(5) * c(5) + c(6) * c(6));

triple = 2. * (c(1) * (c(2) * c(4) + c(3) * c(5)) ...
              + c(6) * (c(2) * c(3) + c(4) * c(5)));

quad = temp(1) * (temp(1) - 2. * temp(2)) + ...
       temp(2) * (temp(2) - 2. * temp(3)) + ...
       temp(3) * (temp(3) - 2. * temp(1));

numer = 32. - square + triple;
denom = 16. - square + 2. * triple + quad;
if (numer < 0.) | (denom < eps1)
   x = gdopmax;
else
   x = sqrt(numer/denom);
end