import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plot_3d_figure():
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')

    R = 3 # Radius of the sphere and cylinder
    color = '#4C99D1'
    alpha = 0.85

    # The object can be interpreted as a combination of a sphere and a half-cylinder.
    # Specifically:
    # 1. A full hemisphere on top (Z >= R)
    # 2. A quarter-sphere in the bottom-front (Z < R, Y < 0)
    # 3. A half-cylinder in the bottom-back (Z < R, Y >= 0)
    
    # 1. Top Hemisphere
    theta1 = np.linspace(0, 2*np.pi, 60)
    phi1 = np.linspace(0, np.pi/2, 30)
    T1, P1 = np.meshgrid(theta1, phi1)
    X1 = R * np.sin(P1) * np.cos(T1)
    Y1 = R * np.sin(P1) * np.sin(T1)
    Z1 = R + R * np.cos(P1)
    ax.plot_surface(X1, Y1, Z1, color=color, alpha=alpha, shade=True, edgecolor='none')

    # 2. Bottom-Front Quarter-Sphere
    theta2 = np.linspace(np.pi, 2*np.pi, 30)
    phi2 = np.linspace(np.pi/2, np.pi, 30)
    T2, P2 = np.meshgrid(theta2, phi2)
    X2 = R * np.sin(P2) * np.cos(T2)
    Y2 = R * np.sin(P2) * np.sin(T2)
    Z2 = R + R * np.cos(P2)
    ax.plot_surface(X2, Y2, Z2, color=color, alpha=alpha, shade=True, edgecolor='none')

    # 3. Bottom-Back Half-Cylinder (Curved Surface)
    theta3 = np.linspace(0, np.pi, 30)
    z3 = np.linspace(0, R, 30)
    T3, Z3_mesh = np.meshgrid(theta3, z3)
    X3 = R * np.cos(T3)
    Y3 = R * np.sin(T3)
    ax.plot_surface(X3, Y3, Z3_mesh, color=color, alpha=alpha, shade=True, edgecolor='none')

    # 4. Bottom flat face of the half-cylinder (Z = 0)
    r4 = np.linspace(0, R, 20)
    theta4 = np.linspace(0, np.pi, 30)
    R4, T4 = np.meshgrid(r4, theta4)
    X4 = R4 * np.cos(T4)
    Y4 = R4 * np.sin(T4)
    Z4 = np.zeros_like(X4)
    ax.plot_surface(X4, Y4, Z4, color=color, alpha=alpha, shade=True, edgecolor='none')

    # 5. Front flat faces (Exposed corners of the half-cylinder at Y = 0)
    # These are the flat vertical areas visible in the front view that form the "U" shape
    x5 = np.linspace(-R, R, 60)
    v5 = np.linspace(0, 1, 20)
    X5, V5 = np.meshgrid(x5, v5)
    Y5 = np.zeros_like(X5)
    # The height of the flat face goes from Z=0 up to the boundary of the quarter-sphere
    Z_max = R - np.sqrt(np.clip(R**2 - X5**2, 0, None))
    Z5 = V5 * Z_max
    ax.plot_surface(X5, Y5, Z5, color=color, alpha=alpha, shade=True, edgecolor='none')

    # Set axes limits
    ax.set_xlim([-R, R])
    ax.set_ylim([-R, R])
    ax.set_zlim([0, 2*R])

    # Set labels
    ax.set_xlabel('X (Width)')
    ax.set_ylabel('Y (Depth)')
    ax.set_zlabel('Z (Height)')

    # Set viewing angle to clearly see the spherical front and the cylindrical base corners
    ax.view_init(elev=25, azim=-55)

    # Force equal proportions for the axes
    ax.set_box_aspect([1, 1, 1])
    
    plt.title("3D Reconstruction: Sphere with Half-Cylinder Base", fontsize=14, pad=20)
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    plot_3d_figure()