1827 - Gauss 內在幾何（絕妙定理：曲率不需要離開曲面）#

import numpy as np

# 1) 圓柱面與平面等距（K = 0，可展開）
def cylinder(u, v, R=1.0):
    """圓柱面：X(u,v) = (R cos v, R sin v, u)"""
    return np.array([R*np.cos(v), R*np.sin(v), u])

def gauss_curvature(X_func, u0, v0, h=1e-4):
    """數值計算高斯曲率 K = (LN - M²)/(EG - F²)"""
    Xu = (X_func(u0+h, v0) - X_func(u0-h, v0)) / (2*h)
    Xv = (X_func(u0, v0+h) - X_func(u0, v0-h)) / (2*h)
    Xuu = (X_func(u0+h, v0) - 2*X_func(u0, v0) + X_func(u0-h, v0)) / h**2
    Xuv = (X_func(u0+h, v0+h) - X_func(u0+h, v0-h) - X_func(u0-h, v0+h) + X_func(u0-h, v0-h)) / (4*h**2)
    Xvv = (X_func(u0, v0+h) - 2*X_func(u0, v0) + X_func(u0, v0-h)) / h**2
    n = np.cross(Xu, Xv); n = n / np.linalg.norm(n)
    E, F, G = Xu@Xu, Xu@Xv, Xv@Xv
    L, M, N = Xuu@n, Xuv@n, Xvv@n
    return (L*N - M**2) / (E*G - F**2)

print("絕妙定理驗證：")
K_cyl = gauss_curvature(cylinder, 0.5, 0.7)
print(f"  圓柱面 K = {K_cyl:.6f}（應 = 0，可展開成平面）\n")

# 2) 球面 K = 1/R²（不可展開）
def sphere(u, v, R=1.0):
    """球面：X(u,v) = (R sin u cos v, R sin u sin v, R cos u)"""
    return np.array([R*np.sin(u)*np.cos(v), R*np.sin(u)*np.sin(v), R*np.cos(u)])

for R in [1.0, 2.0, 5.0]:
    K_sp = gauss_curvature(lambda u, v, R=R: sphere(u, v, R), 0.7, 0.7)
    print(f"  球面 R={R}: K = {K_sp:.6f}（理論 = {1/R**2:.6f}）")
print("  → 球面 K > 0：不能展開成平面（地圖必失真）\n")

# 3) 主曲率與 Euler 公式
print("Euler 公式驗證（主曲率）：")
# 圓柱面：k1 = 0（沿軸），k2 = 1/R（環向）
k1, k2, R = 0.0, 1.0, 1.0
for theta in [0, np.pi/6, np.pi/4, np.pi/2]:
    kn = k1*np.cos(theta)**2 + k2*np.sin(theta)**2
    print(f"  θ={np.degrees(theta):5.1f}°: k_n = {kn:.6f}")
print(f"  K = k1·k2 = {k1*k2:.6f}（圓柱面恆為 0）\n")

# 4) Gauss–Bonnet：球面總曲率恆為 4π
print("Gauss–Bonnet 定理：")
# 數值積分球面曲率
us = np.linspace(0.001, np.pi-0.001, 100)
vs = np.linspace(0, 2*np.pi, 200)
total = 0.0
for i in range(len(us)-1):
    for j in range(len(vs)-1):
        du, dv = us[i+1]-us[i], vs[j+1]-vs[j]
        u0, v0 = (us[i]+us[i+1])/2, (vs[j]+vs[j+1])/2
        Xu = (sphere(u0+1e-5, v0) - sphere(u0-1e-5, v0)) / 2e-5
        Xv = (sphere(u0, v0+1e-5) - sphere(u0, v0-1e-5)) / 2e-5
        area_el = np.linalg.norm(np.cross(Xu, Xv))
        K = 1.0  # 球面 K = 1
        total += K * area_el * du * dv  # dA = |Xu×Xv| du dv
print(f"  球面總曲率 ∫K dA = {total:.4f}")
print(f"  理論值 4π = {4*np.pi:.4f}")
print("  → 無論怎麼變形球面，總曲率恆為 4π（χ=2 的必然）")

絕妙定理驗證：
  圓柱面 K = 0.000000（應 = 0，可展開成平面）

  球面 R=1.0: K = 1.000000（理論 = 1.000000）
  球面 R=2.0: K = 0.250000（理論 = 0.250000）
  球面 R=5.0: K = 0.040000（理論 = 0.040000）
  → 球面 K > 0：不能展開成平面（地圖必失真）

Euler 公式驗證（主曲率）：
  θ=  0.0°: k_n = 0.000000
  θ= 30.0°: k_n = 0.250000
  θ= 45.0°: k_n = 0.500000
  θ= 90.0°: k_n = 1.000000
  K = k1·k2 = 0.000000（圓柱面恆為 0）

Gauss–Bonnet 定理：
  球面總曲率 ∫K dA = 12.5669
  理論值 4π = 12.5664
  → 無論怎麼變形球面，總曲率恆為 4π（χ=2 的必然）