Volume of a bead with square hole- Program approach?

[Albert Chan]

Hi, DM48

Nice to know there is another way to get volume:

sphere - hole = (x-strip - hole) + (2 spherical caps)

Since we have direct formula for spherical cap, 2nd integral can be skipped
2 caps = 2 * pi/3*h^2*(3r-h) = pi/12*(2h)^2*(3d-2h) = pi/12*(d-a)^2*(2d+a)

10 DEF FNH(X) @ N=N+1 @ FNH=SQRT(H-X*X) @ END DEF
20 DEF FNZ(Y) @ H=D*D-Y*Y @ FNZ=INTEGRAL(0,FNH(B),P,FNH(IVAR)-B) @ END DEF
30 INPUT "D,A,B = ";D,A,B @ N=0 @ P=1E-5
40 I1=INTEGRAL(0,A,P,FNZ(IVAR)) @ I2=PI/12*(D-A)^2*(2*D+A)
50 DISP I1;"+";I2;"=";I1+I2,N @ GOTO 30

>run
D,A,B = 12,2,2
 176.547653773 + 680.678408277 = 857.22606205       960
D,A,B = 12,2,10
 17.2338339713 + 680.678408277 = 697.912242248      992
D,A,B = 12,10,2
 662.307552043 + 35.6047167408 = 697.912268784      1888

Numerically, thinner, less curvy strip (A ≤ B) converge better.

For comparison, below patch directly get (-hole) + sphere

>20 DEF FNX(Y) @ H=D*D-Y*Y @ FNX=INTEGRAL(0,B,P,FNH(IVAR)) @ END DEF
>40 I1=-INTEGRAL(0,A,P,FNX(IVAR)) @ I2=PI/6*D*D*D

>run
D,A,B = 12,2,2
-47.552634005  + 904.778684234 = 857.226050229      225
D,A,B = 12,2,10
-206.866458199 + 904.778684234 = 697.912226035      465
D,A,B = 12,10,2
-206.866453454 + 904.778684234 = 697.91223078       465