% 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