function [REL_AREAS] = vorareas(M);
% vorareas(M) ----
% This function will take standardized point configurations, and make voronoi 
% tessellations of them.  The edges are bounded artificially, by putting an octagon of
% points around the field, far outside of it.  These will absord the infinited edge
% cells, leaving the points of interest with finite areas.  This will not 
% come out perfectly or consistently shaped.  So, RELATIVE "domains of danger" will 
% be found, since the tessellation will have different total areas for each case.  At 
% least the areas are finite.  
%
% 
if nargin == 1
    plot_req = 0; 
elseif strcmp(p,'y') ==1
    plot_req = 1;
end

% First, some parameters of the field.  The 'width' is the actual half-length of the 
% field dimensions; and L will be the distance out that the artificial
%  edges occur. 
W = sqrt(9*100*pi);
L = W*4;
OB = (1/sqrt(2)) * L;

% The boundary is made by placing an octogon around the field, far outside of it.
% BOUND is a matrix containing the coordinates of the octogon vertices.
BOUND = [0 0 -L 0; 0 0 L 0; 0 -L 0 0; 0 L 0 0; 0 -OB -OB 0; 0 OB OB 0; 0 OB -OB 0;...
         0 -OB OB 0];

% The edge points are added to the input array M, then the x and y coordinates 
% extracted from columns 2 and 3, since this 2-column format is needed
% for the 'voronoin' function.
N0 = [M; BOUND];
N = [N0(:,2) N0(:,3)];


% Now, the voronoi calculations can be made.  The output V will be a matrix containing
% the coordinates of all vertices for the polygons; C is a cell, each element 
% (rowwise) being a matrix containing the indices into V of the necessary vertices 
% for the polygon around that point. 
[V C]  = voronoin(N);
    

% Then, for each point i, the polygon is constructed from the voronoi data in C and
% V, and the area computed by convhulln.  K will just be a repeat of the vertices, 
% but A will be the area and therefore will be saved.  By basing the no. of reps on 
% size(M), it will ignore the tessellations around the edge points (i.e., the 
% infinite throwaways).

AREAS = [];                             % Initialize matrix for area storage.
for i = 1:size(M,1)                     % Loop over each point.
    
    X = V(C{i},:);                      % All the necessary vertices are put into 
                                        %  matrix X.
    
    [K A] = convhulln(X);               % Area is computed for the polygon.
    
    AREAS = [AREAS; A];                 % Save the areas into the matrix.
end

% Areas are added up, and the 'RELative AREAS' matrix initialized, for storage of, 
% you guessed it, the relative areas of each polygon, for each point.
TOTAL_A = sum(AREAS);


REL_AREAS = [];

% Calculate relative areas, and store them.
for j = 1:size(AREAS,1)                 % Loop through each polygon.
    
    relarea = AREAS(j,1) / TOTAL_A;     % Calculate rel a.
    
    REL_AREAS = [REL_AREAS; relarea];   % Store each in the matrix.
end

% Values are multiplied by 100 (i.e. put in % form) to make them easier 
% to read.
REL_AREAS = REL_AREAS * 100;