Cube-Sphere Intersection Volume

This week's math problem was particularly difficult. I wanted to know the intersection volume between a cube and a sphere. Imagine a sphere contained in a cube. The sphere will first touch the cube when its radius is equal to half the width of the cube. Until that point, its volume is just the formula for the volume of a sphere. Once the sphere's radius is 2 sqrt(3) times the width of the cube, then it completely engulfs the cube and the intersection of the 2 objects is just the volume of a cube. In between, I've found it looks like this:

The gray vertical lines represent when to use which equation. There are in fact 2 other regimes besides the two I discussed. I'm going to discuss these regimes as related to the ratio of the rho and d where the first is the radius of the sphere and the second is half the cube width. Let's see them visually.

The first regime is the spherical, where rho / d < 1. It looks like a sphere:

The second regime has circular "lenses" cut-off from the sphere and occurs while rho / d < sqrt(2):

The third regime is characterized by filling the 8 corners of the cube:

Finally, the last regime is just a cube. This occurs at rho / d > sqrt(3):


The equation for the sphere regime is:

The equation for the "lens" regime may be found by calculating the volume of the removed lenses and then subtracting that from the volume of the sphere. First, look at the cross section of the lens:

We can use those rotation transformation we learned in calculus. That is described by the equation:

The last regime is much more difficult to deal with. The sphere is first broken into 16ths and then each piece of that is divided into 5 different regions. The 16th piece may be viewed as an animation here to visualize it. I'll break down each region with an animation showing how it changes as a function of the ratio of the radius to the cube width and the resulting volume of the area.
 The first region is defined by:

Here is a gif of the region:

The second region is defined by:

Here is a gif of the region:

The third region is defined by:

Here is a gif of the region:

The fourth region is defined by:

Here is a gif of the region:

The fifth region is defined by:

Here is a gif of the region:

Finally, all the regimes/regions may be combined into a python script which calculates the cube-sphere intersection volume. I should mention that I'm not that good at calculus and I think the two unevaluated integrals have analytic forms which I just can't find yet. That's one of the problems of doing integration by hand; you have to be better than me. I'm considering investing in a Mathematica license after doing these integrals. Here's the python code:

	#! /usr/bin/python
	from scipy.integrate import quad
	import numpy as np
	from math import pi
	#square regions
	def region_1(rho, d):
	return 0.5 * (d**2 * np.sqrt(rho**2 - 2 * d ** 2))
	def region_1_2_theta(rho, d):
	return np.arccos(d / np.sqrt(rho ** 2 - d ** 2))
	#curved square regions
	def region_2_integrand(theta, rho, d):
	return np.sqrt(np.cos(theta) ** 2 - (d / rho)**2) / (np.cos(theta) ** 3)
	def region_2(rho, d):
	i4 = d**3 / 6. * (rho**2 / d **2 - 1) * (pi / 4 - region_1_2_theta(rho, d))
	i3 = d ** 2 * rho / 3. * quad(region_2_integrand, region_1_2_theta(rho,d), pi / 4, args=(rho, d))[0]
	return i4 + i3
	#spherical region
	def region_3_integrand(theta, rho, d):
	return np.sqrt(np.cos(theta) ** 2 - (d / rho)**2) / np.cos(theta)
	def region_3(rho, d):
	return rho ** 3 / 3. * (d / rho * (pi / 4 - region_1_2_theta(rho, d)) - quad(region_3_integrand, region_1_2_theta(rho, d), pi / 4, args=(rho, d))[0])
	def calc_volume(rho, d):
	alpha = rho / d
	if(alpha <= 1):
	return 4./3 * pi * (rho) ** 3
	if(alpha <= np.sqrt(2)):
	return 4. / 3 * pi * (rho) ** 3 - 6. * pi * (2 * (rho)**3 - 3 * d * (rho)**2 + d**3) / 3.
	if(alpha < np.sqrt(3)):
	return 16. * (region_1(rho,d) + region_2(rho, d) + region_3(rho,d))
	return 8. * d ** 3
	for rho in np.arange(1, np.sqrt(3), 0.001):
	print "%g %g" % (rho, calc_volume(rho, 1))

view raw gistfile1.py hosted with ❤ by GitHub