function [final, Total_V] = rotatevor3d(M);
% myvoronoi3d(M) ----  Calculates Voronoi volumes for an array of points
% in 3D.  Analogous to the 'myvoronoi' function, with an extra column in 
% the input for z-coordinates.  Here also, the points-of-interest are 
% bounded far outside the field-of-interest to avoid having infinite volumes.
% The vertices of a dodecahedron form the boundary.
%
% The input matrix 'M' should have individual numbers in the 1st column, then
% x y and z coordinates in the 2nd, 3rd, and 4th columns respectively.
%
% Output are the relative volumes (in a vector) of the polyhedra for each 
% point, i.e. the percent-area each occupies.


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

% Next, the boundary points will be placed, far outside of the space of interest.  
% The vertices of a dodecahedron are used, 20 points total, and are placed in the 
% 'BOUND' matrix.

% find the max radius (the longest straight-line distance between any two
% points in the distribution

currentMax = 0;
j=0;

for i=1:(size(M,1))
    for j=1:(size(M,1))
        if dist3d(M(i,2), M(i,3), M(i,4), M(j,2), M(j,3), M(j,4))>currentMax
            currentMax = dist3d(M(i,2), M(i,3), M(i,4), M(j,2), M(j,3), M(j,4));
        end
    end
end

radius = currentMax;

% find the center of mass by summing the x, y, and z values of each point,
% then dividing by the number of points

currentX = 0;
currentY = 0;
currentZ = 0;

for i=1:(size(M,1))
    currentX = (M(i,2))+currentX;
    currentY = (M(i,3))+currentY;
    currentZ = (M(i,4))+currentZ;
end

centerX = currentX/(size(M,1));
centerY = currentY/(size(M,2));
centerZ = currentZ/(size(M,3));
center = [centerX, centerY, centerZ];

radius = currentMax;

% define a regular dodecahedron, centered at the origin with radius of one.

BOUND = [[.3568220896, -.9341723582, 0.]; [-.3568220896, -.9341723582, 0.]; [-.5773502686, -.5773502686, .5773502686]; [0., -.3568220896, .9341723586]; [.5773502686, -.5773502686, .5773502686]; [0., .3568220896, -.9341723582]; [-.5773502686, .5773502686, -.5773502686]; [-.3568220896, .9341723586, 0.]; [.3568220896, .9341723586, 0.]; [.5773502686, .5773502686, -.5773502686]; [0., -.3568220896, -.9341723582]; [-.5773502686, -.5773502686, -.5773502686]; [-.9341723582, 0., -.3568220896]; [-.9341723582, 0., .3568220896]; [-.5773502686, .5773502686, .5773502686]; [0., .3568220896, .9341723586]; [.5773502686, .5773502686, .5773502686]; [.9341723586, 0., .3568220896]; [.9341723586, 0., -.3568220896]; [.5773502686, -.5773502686, -.5773502686]];

% shift the center to the center of mass, as determined above.

for i = 1:size(BOUND,1)
    BOUND(i,1)=BOUND(i,1)+centerX;
    BOUND(i,2)=BOUND(i,2)+centerY;
    BOUND(i,3)=BOUND(i,3)+centerZ;
end

% scale the polyhedron to 1.5 times the largest radius within the group

BOUND = BOUND*1.5*radius; 

% generate a random rotation matrix using the rotate3d.m script

k = rotate3d;

% generate the final bounding matrix

BOUND = BOUND*k;

% A fourth column is added to 'BOUND', so that the input matrix M may be
% attached correctly.  
BOUND = [zeros(20,1) BOUND];



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


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

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

VOL = [];                             % 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.
    for j=1:size(X,1)
        for k=1:size(X,2)
    
    [K A] = convhulln(X);               % Volume is computed for the polyhedron.
    
    VOL = [VOL; A];                 % Save the volumes into the matrix.
end

% Volumes are added up, as TOTAL_V,
final = [VOL];
TOTAL_V = sum(VOL)