Numerical Integration Exercise - Common Region Volume Computation

I triple-integrated over the cone interior, using cylindrical coordinates $(r,\theta,z)$ . An infinitesimal volume from the cone is added to a running total if that volume also resides within the sphere. I ran trials with the cone divided into $10^6$ pieces in the lowest-resolution case, and $10^9$ pieces in the highest-resolution case. The common volume calculation converges to somewhere around $347$ or $348$ . Results are printed at the end of the code.

I also considered solving this again using Monte Carlo integration to double check, but I was too lazy.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

import math

N = 1000   # spatial resolution parameter (cone divided into N^3 pieces)

Rc = 5.0*math.sqrt(2.0)  # cone base radius
H = 20.0  # cone height
Rs = 6.0  # sphere radius

xs = 5.0  # sphere center coordinates
ys = 0.0
zs = 6.0

###################################

dz = H/N                # infinitesimals
dtheta = 2.0*math.pi/N

V = 0.0
z = 0.0

# Integrate over interior of cone
# in cylindrical coordinates (r,theta,z)

while z <= H:               

    rmax = Rc*(1.0 - z/H)  # max radius for a particular z

    dr = rmax/N

    r = 0.0

    while r <= rmax:

        theta = 0.0

        while theta <= 2.0*math.pi:

            x = r*math.cos(theta)
            y = r*math.sin(theta)

            # Square of distance from point inside cone to sphere center

            Dsq = (x-xs)**2.0 + (y-ys)**2.0 + (z-zs)**2.0

            # If cone interior point is within sphere, add infinitesimal volume

            if Dsq <= Rs**2.0:

                dV = r*dr*dtheta*dz
                V = V + dV

            theta = theta + dtheta

        r = r + dr

    z = z + dz

###################################

print N
print V

#>>> 
#100
#356.343972812
#>>> ================================ RESTART ================================
#>>> 
#200
#351.703775033
#>>> ================================ RESTART ================================
#>>> 
#400
#349.459395077
#>>> ================================ RESTART ================================
#>>> 
#500
#347.176253307
#>>> ================================ RESTART ================================
#>>> 
#1000
#348.099929886
#>>>