Cube-sphere intersection

gistfile1.py

#! /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))