% radius_scaling_factor = calc_rad_scal_factor(packing_type, porosity)
% calculates scaling factor for sphere radius (or diameter) to
% obtain required porosity. As input it takes packing type (BCC, FCC, SC)
% and porosity. The function is based on calculation of unknown radius
% assuming that a unit cell is formed by spheres of unit diameter. If
% there is no overlap between spheres, radius calculation is
% straightforward. In the case of overlaps, it is assumed that only pair
% overlaps exist (i.e., there is no space points which are simultaneously
% in three spheres). Then the volume _not_ occupied by spheres is formed by
% volume of each individual sphere reduced by the volume of spherical caps
% due to overlap of a given sphere with its neighbours. The final solid
% volume is constructed in described way and then cubical equation is solved
% numericaly to find unknown radius given input porosity.
function [r_scal_factor] = calc_rad_scal_factor(pack_type, porosity)
density = 1 - porosity;
d = 1.0; % intersphere distance for touching spheres
switch pack_type
case 1
amount_of_contacts = 8;
Vcell = (2.0 / sqrt(3))^3;
Nspheres_per_cell = 2;
touching_density = pi*sqrt(3)/8;
case 2
amount_of_contacts = 12;
Vcell = sqrt(2)^3;
Nspheres_per_cell = 4;
touching_density = pi/sqrt(18);
case 3
amount_of_contacts = 6;
Vcell = 1.0;
Nspheres_per_cell = 1;
touching_density = pi/6;
end
if density < (touching_density + eps*1e2)
r_scal_factor = ((1 - porosity) / touching_density) ^ (1/3);
else
Vspheres = Vcell * density;
sphere_caps_to_remove = amount_of_contacts*Nspheres_per_cell;
% from http://en.wikipedia.org/wiki/Cubic_function#General_formula_for_roots
% V = 1/12 * pi * (4R + d) * (2R - d)^2
% where V is the volume of two overlapping spheres, R their equal radius
% (= 0.5 in the case of touching spheres), and d is the separation
% distance between sphere centers; d = 1.0 for SC packing
% let us construct the polynomial in the form of a*x^3 + b*x^2 + c*x + d = 0
cubic_poly = [Nspheres_per_cell*4/sphere_caps_to_remove - 2; ...
2*d - d/2; ...
0; ...
-(d/2)^3 - 3*Vspheres/(sphere_caps_to_remove*pi) ];
r = roots(cubic_poly);
% we are intrested only in the smallest positive root: after
% intersection, the joint shape (formed for example by two spheres) will
% have equal volume for two radius values, the one closer to initial
% sphere radius of 0.5 and much larger, closer to unity
r_scal_factor = min( r(r > 0) ) * 2; % 2 due to using radius above
end