there are many bright people who have studied the problem you attempt to address.
As a warm-up I suggest have a look at the following:
1. http://mathworld.wolfram.com/CirclePacking.html
2. http://mathworld.wolfram.com/SquarePacking.html
3. http://www2.stetson.edu/~efriedma/packing.html
You want to approximate a sphere with cuboids, the 3D version of
4. http://www2.stetson.edu/~efriedma/squincir/
I assume without 'piercing' the surface of the sphere, otherwise for r=1, 1 with square L=sqrt(pi) has exactly same area, same would apply to any higher amount of squares always reaching zero error.
From the Stetson series of interest, since there is only mathematical proof of the 1st 2 cases, 1 square and 2 squares, I have selected those showing symmetry around x axis:
amount of squares
s1=[1 2 4 5 7 12 14 21 26 31 33 35]
instead of growing circle, normalizing Stenson's to fixed circle radius, I get:
s2=[1 5^.5/2 2^.5 10^.5/2 13^.5/2 5^.5 2.5 34^.5/2 41^.5/2 5*2^.5/2 13^.5 53^.5/2]
I just search for a short time, but did match s2 any known series, but s2 looks quite the same as page 12 of http://www.winton.phy.cam.ac.uk/DanielFletcher
showing the volume of embryo cells over amount of cells, while inside the embryo.
Obviously, you want a more general case, using cuboids instead increases variables
In MATLAB you can build a cuboid with:
xc=0;yc=0 % define where you want it centered
Dx=1;Dy=2;Dz=3 % define dimensions
dx=[-Dx/2 Dx/2];dy=[-Dy/2 Dy/2];dz=[-Dz/2 Dz/3]
[x,y,z] = meshgrid(dx,dy,dz); % generate cuboid, not centered
x = [x(:);0];y = [y(:);0];z = [z(:);0];
DT = delaunayTriangulation(x,y,z);
h1=tetramesh(DT);
camorbit(20,0)
Another easier way to build a cuboid is:
vert = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];
fac = [1 2 6 5;2 3 7 6;3 4 8 7;4 1 5 8;1 2 3 4;5 6 7 8];
patch('Vertices', vert, 'Faces', fac, 'FaceVertexCData', hsv(6), 'FaceColor', 'flat')
view(3)
axis([-3 3 -3 3 -3 3]);grid on
The sphere can be approximated with slices, search contourslice
x = -2:0.2:2;
y = -2:0.25:2;
z = -2:0.16:2;
[X,Y,Z] = meshgrid(x,y,z);
[Xi,Yi,Zi] = sphere;
V = X.*exp(-X.^2-Y.^2-Z.^2);
contourslice(X,Y,Z,V,Xi,Yi,Zi)
view(3)
Volume of sphere is Vol_sphere=4/3*pi*R^3
Volume of your approximation, just with cubes all same L side: for each slice Vol(i)=N(i)*L^2
where N(i) is the amount of cubes in slice I.
Vol_total=sum(Vol)
error=abs(4/3*pi-Vol_total)
If you find anything useful to start, I would appreciate you let me know about your progress.
Regards
John
