{}
import numpy as np
# Parameters
R = 20.37 # Sphere radius
h = 16.5 # Cap height
a = 20 # Base radius
phi_max = 1.38 # Maximum polar angle for the cap
# Discretization
num_theta = 1000 # Number of steps in theta
num_phi = 500 # Number of steps in phi
dtheta = 2 * np.pi / num_theta
dphi = phi_max / num_phi
# Initialize total surface area
total_area = 0
# Numerical integration
for i in range(num_theta):
theta = i * dtheta
for j in range(num_phi):
phi = j * dphi
# Compute x, y, z
x = 20 + R * np.sin(phi) * np.cos(theta)
y = 20 + R * np.sin(phi) * np.sin(theta)
z = 3.87 + R * np.cos(phi)
# Apply second cut (y in [0, 6.5] and x not in [0, 9.5] or [30.5, 40])
if 0 <= y < 6.5:
if not (0 <= x <= 9.5 or 30.5 <= x <= 40):
continue
# Apply first cut (y <= 27)
if y >= 27:
continue
# Add patch area to total
total_area += R**2 * np.sin(phi) * dtheta * dphi
print("Total Surface Area:", total_area)