Volume of a bead with square hole- Program approach?

[Albert Chan]

Code:

function hv2(d,a,b) -- sphere, rectanglar hole removed volume
    local d2, a2, b2 = d*d, a*a, b*b
    local k = sqrt(d2 - a2 - b2)
    d = a*b*k - d*d2*atan((a*b)/(d*k))
    return (2*d + a*(3*d2-a2)*atan(b/k) + b*(3*d2-b2)*atan(a/k)) / 6
end

-- fast estimate: m1 < hv2 < m2
function hv2_m2(d,a,b) return a*b * sqrt(d*d - (a*a+b*b)/3) end
function hv2_m1(d,a,b) return a*b * (2*d + sqrt(d*d-a*a-b*b))/3 end

Although estimation formula assumed square hole, inequality work with rectangular hole too.
For estimate, hv2_m2(d,a,b) is recommended: cheaper to compute, better accuracy.