The script is to generate the grain size distribution of the read particles

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear all

tic


NoOfParticles=dir('ClosedSurface\*.csv');
numberOfcsvs=numel(NoOfParticles);

for i=1:numberOfcsvs
csvName=strcat('ClosedSurface\','surfacevoxel-',num2str(i,'%02i'),'.csv');
I{i}=xlsread(csvName);
x{i} = I{i}(:,1);
y{i} = I{i}(:,2);
z{i} = I{i}(:,3);
% Mass Center
XXX(i) = sum(x{i})/numel(x{i}); %non zero numbers of the surface.
YYY(i) = sum(y{i})/numel(y{i});
ZZZ(i) = sum(z{i})/numel(z{i});

% Transfer to around origin

x{i} = x{i} - XXX(i);
y{i} = y{i} - YYY(i);
z{i} = z{i} - ZZZ(i);

% Ellipsoid fitting
% fit ellipsoid in the form Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz = 1
D{i} = [ x{i} .* x{i}, y{i} .* y{i}, z{i} .* z{i}, 2 * x{i} .* y{i}, ...
2 * x{i} .* z{i}, 2 * y{i} .* z{i}, 2 * x{i}, ...
2 * y{i}, ...
2 * z{i} ];

% solve the normal system of equations
v{i} = ( D{i}' * D{i} ) \ ( D{i}' * ones( size( x{i}, 1 ), 1 ) );
% Algebraic form of the ellipsoid
E{i} = [ v{i}(1) v{i}(4) v{i}(5) v{i}(7); ...
v{i}(4) v{i}(2) v{i}(6) v{i}(8); ...
v{i}(5) v{i}(6) v{i}(3) v{i}(9); ...
v{i}(7) v{i}(8) v{i}(9) -1 ];
% find the center of the ellipsoid
center{i} = -E{i}( 1:3, 1:3 ) \ [ v{i}(7); v{i}(8); v{i}(9) ];
% form the corresponding translation matrix
T = eye( 4 );
T( 4, 1:3 ) = center{i}';
% translate to the center
R{i} = T * E{i} * T';
% solve the eigenproblem
% eigenvectors corresponding to the semi axis directions
[ evecs{i}, evals{i} ] = eig( R{i}( 1:3, 1:3 ) / -R{i}( 4, 4 ) );
% Three semi-axis
radii{i} = sqrt( 1 ./ diag( evals{i} ) );

% Particle orientations of all the particles
ParticleOrientation(i,:)= evecs{i}(1,:);
% Particle Semi-Axes of all the particles
ParticleSemiAxes(i,:)= radii{i}';
end

% Set Major Semi-axis as the approximation of particle sizes
ParticleSize = ParticleSemiAxes(:,2);
[ParticleSize,Index] = sort(ParticleSize);
ParticleVolume = 4*pi*ParticleSemiAxes(:,1).*ParticleSemiAxes(:,2).*ParticleSemiAxes(:,2)/3;

for i = 1: numel(Index)
ParticleVolume1(i,:) = ParticleVolume(Index(i));
end

%-------------------------------------**
% For CDF calculation

%Real particle Volume
ParticleVolumeR = ParticleVolume1/(795.6/7)^3;

% Particle Size measured by pixels
% xlswrite('ParticleSizePixel.xlsx',ParticleSize);

% Particle Size measured by millimeters
% xlswrite('ParticleSizeMM.xlsx',sort(ParticleSemiAxes(:,2))*7*2/795.6);

% Particle size Normalization as sieve openings

ParticleSize = sort(ParticleSemiAxes(:,2))*7*2/795.6;

SieveSize = ParticleSize/(0.5+1/1.414);

%binranges
binranges = 0.4:0.02:0.7;

bincounts = histc(SieveSize, binranges);

%Frequecy of each bin
SievePro = bincounts/numel(SieveSize);

% Cummulative Probability
SieveProb = cumsum(SievePro);

figure
ylim([0,1]);
hold all
p = polyfit(binranges',SieveProb,6);
f = polyval(p,binranges);

semilogx(binranges,SieveProb,'o',binranges,f,'-');
title('Grain Size Distribution of Ottawa 20-30 Sand','fontsize',13)
xlabel('Particle Size(mm)','fontsize',13);
ylabel('Percent Finer by Volume(Weight)(%)','fontsize',13);
set(gca,'fontsize',12);
set(gca,'xscale','log');
xscale log;
xlim([0.1,1]);
grid on
legend('Separated Particles','polynomial Fitting of the separated particles')

toc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     Grain size distribution of Ottawa 20-30 sands