% GPSCEP.m
% Author: Jim Connor (jimconnor99@hotmail.com)


% ********************************************************************************
% Code to calculate the CEP. Basic idea is to find the radius of a circle
% where half of the data points lie. The algorithm starts by calculating the
% distances of all the points. Those distances are then compared to a given radius
% and if half of the points are inside the radius, then that radius is the CEP.
% CEP is defined as the circle where 50% of the points lie
% ********************************************************************************

for i=1:size(north)  
  x(i)= sqrt(north(i)^2 + east(i)^2);   % Calculate distances
end


for r=1:2000     % increasing radius, starts at 1 meter
    
    if R==1
        r=r/100;        % Allows accuracy of two decimal points for RTCM CEPs
    else
        r=r/1000;       % Allows accuracy of three decimal points for RT20 CEPs
    end
    
    inside =0;    % resets the 'number of points inside the radius' count
             
       for j=1:size(north)
           if (x(j) <= r ), inside=inside + 1; end  % if the distance is smaller than the radius, count it
           if inside >= (size(north)/2) break; end  % if half of the points are counted stop
       end
               
  if inside >= (size(north)/2) break; end   % used to break the outer loop when finished

end