#!/usr/bin/python
# -*- coding: utf-8 -*-
#silly sphere minus hole problem for IFLS facebook group
# (shame on me for not being bothered to do this analytically)
from numpy import *
from matplotlib.pyplot import *

rs = linspace(3,100,100); #some trial radii, must be rs > 3 but limits to infinite

l = 6.0; #length of 'hole'
Vs = 4*pi*rs**3/3; # volume of sphere

rc = sqrt(rs**2 - (l/2)**2); #radius of hole
h = rs - l/2;  #height of end cap
Vc = l*pi*rc**2; # volume of cylinder 
Ve = pi*h**2/3 * (3*rs - h); # volume of one end cap

Vt = Vs -Vc -2*Ve;

clf();
plot(rs, Vs, "b");
plot(rs, Vc, "r");
plot(rs, 2*Ve, "m");
plot(rs, Vt, "k");
legend(["sphere volume","hole volume","missing end cap volume","sphere - cyld - ends"]);
draw();