6.3 龐加萊圓盤模型的 SymPy 實驗（1882 年）#

from sympy import *
from sympy.simplify.fu import fu

z, w = symbols('z w', complex=True)
xr, yr = symbols('xr yr', real=True)
theta = symbols('theta', real=True)

print("=== 3.1 莫比烏斯變換 w = i(1+z)/(1-z) ===")
W = I*(1 + z)/(1 - z)
print("w(z)  =", W)
print("dw/dz =", simplify(diff(W, z)))
u = 1/(z - 1)**2
print("|dw/dz|^2 - 4/|1-z|^4 =",
      simplify(((re(u))**2 + (im(u))**2 - 1/((1 - xr)**2 + yr**2)**2).subs(z, xr + I*yr)))

print()
print("=== 3.2 圓盤內 ↔ 上半平面 ===")
Wxy = W.subs(z, xr + I*yr)
print("Im(w) =", simplify(im(Wxy)))
print("圓盤邊界 |z|=1 時 Im(w) =", simplify(im(Wxy).subs(yr**2, 1 - xr**2)))
print("往返檢驗 w((w-i)/(w+i)) - w =", simplify(W.subs(z, (w - I)/(w + I)) - w))

print()
print("=== 3.3 兩個模型的度量被這張地圖保持 ===")
lhs = ((re(diff(W, z)))**2 + (im(diff(W, z)))**2)/im(W)**2
rhs = 4/(1 - xr**2 - yr**2)**2
print("|dw|^2/Im(w)^2 - 4/(1-|z|^2)^2 =", simplify(lhs.subs(z, xr + I*yr) - rhs))

print()
print("=== 3.4 測地線圓弧 → 半圓 / 垂直線 ===")
w1 = simplify(W.subs(z, cos(theta) + I*sin(theta)))
print("w(e^{iθ}) 的虛部 =", simplify(im(w1)), "  → 單位圓被送到實軸")
print("w(e^{iθ}) 的實部 + cot(θ/2) =", fu(re(w1) + cot(theta/2)), "  → w = -cot(θ/2)")

an = 2
xi = Rational(1, an)
yi = sqrt(1 - xi**2)
for zz in [xi + I*yi, xi - I*yi]:
    val = W.subs(z, zz).evalf()
    print("圓心 (2,0)、半徑 sqrt(3) 的測地線，與單位圓交點 z =", zz,
          " → w =", N(re(val)), "+", N(im(val)), "i")
print("→ 上半平面中以 0 為圓心、半徑 sqrt(3) =", sqrt(3).evalf(), " 的半圓")
print("直徑（實軸上的線段）→ 虛軸：z=0 -> w =", W.subs(z, 0),
      " ; z=1/2 -> w =", simplify(W.subs(z, Rational(1, 2))),
      " ; Re =", simplify(re(W.subs(z, Rational(1, 2)))))

=== 3.1 莫比烏斯變換 w = i(1+z)/(1-z) ===
w(z)  = I*(z + 1)/(1 - z)
dw/dz = 2*I/(z - 1)**2
|dw/dz|^2 - 4/|1-z|^4 = 0

=== 3.2 圓盤內 ↔ 上半平面 ===
Im(w) = (-yr**2 - (xr - 1)*(xr + 1))/(yr**2 + (xr - 1)**2)
圓盤邊界 |z|=1 時 Im(w) = 0
往返檢驗 w((w-i)/(w+i)) - w = 0

=== 3.3 兩個模型的度量被這張地圖保持 ===
|dw|^2/Im(w)^2 - 4/(1-|z|^2)^2 = 0

=== 3.4 測地線圓弧 → 半圓 / 垂直線 ===
w(e^{iθ}) 的虛部 = 0   → 單位圓被送到實軸
w(e^{iθ}) 的實部 + cot(θ/2) = 0   → w = -cot(θ/2)
圓心 (2,0)、半徑 sqrt(3) 的測地線，與單位圓交點 z = 1/2 + sqrt(3)*I/2  → w = -1.73205080756888 + 2.71050543121376e-20 i
圓心 (2,0)、半徑 sqrt(3) 的測地線，與單位圓交點 z = 1/2 - sqrt(3)*I/2  → w = 1.73205080756888 + 2.71050543121376e-20 i
→ 上半平面中以 0 為圓心、半徑 sqrt(3) = 1.73205080756888  的半圓
直徑（實軸上的線段）→ 虛軸：z=0 -> w = I  ; z=1/2 -> w = 3*I  ; Re = 0

from sympy import *

x, y, a, p = symbols('x y a p', positive=True)

print("=== 4.1 正交條件：|c1-c2|^2 = r1^2 + r2^2 ===")
rho = sqrt(a**2 - 1)                    # 由 a^2 = 1 + rho^2 得來
print("候選半徑 rho =", rho)
print("驗算 a^2 - (1 + rho^2) =", simplify(a**2 - (1 + rho**2)))

xi = Rational(1, 1)/a
yi = sqrt(1 - xi**2)
r1 = Matrix([xi, yi])                   # O -> 交點
r2 = Matrix([xi - a, yi])               # C2 -> 交點
print("兩半徑向量內積 =", simplify(r1.dot(r2)))
print("切線夾角 =", simplify(acos(r1.dot(r2)/(sqrt(r1.dot(r1))*sqrt(r2.dot(r2))))))

print()
print("=== 4.2 陷阱：圓心 (a,0)、半徑 sqrt(1+a^2) 的圓過哪些點？ ===")
rho_plus = sqrt(1 + a**2)
print("(x-a)^2+y^2 = 1+a^2  令 x=0 → y^2 =", simplify(rho_plus**2 - a**2))
print("→ 它通過單位圓上的兩個極點 (0, ±1)")
v1 = Matrix([0, 1]); v2 = Matrix([-a, 1])
print("在 (0,1) 處兩半徑向量內積 =", simplify(v1.dot(v2)), "（正交應為 0）")
ang = acos(v1.dot(v2)/(sqrt(v1.dot(v1))*sqrt(v2.dot(v2))))
print("它與單位圓的切線夾角 =", simplify(ang), " ; a=2 時 =",
      N(ang.subs(a, 2)), "rad =", N(deg(ang.subs(a, 2))), "度")

print()
print("=== 4.3 通過圓盤內定點 P=(p,0)、且圓心在 x 軸上的測地線 ===")
cx, rr = symbols('cx rr', positive=True)
sol = solve([Eq(cx**2, rr**2 + 1),                   # 與單位圓正交
             Eq((p - cx)**2, rr**2)],               # 通過 P
            [cx, rr], dict=True)
for s in sol:
    print("解:", {k: simplify(v) for k, v in s.items()})
print("圓心的另一個表達：半徑 rr 的測地線，其圓心距原點 = sqrt(1+rr^2) =",
      simplify(sqrt(1 + sol[0][rr]**2)))
print("與 x 軸的另一交點 =", simplify(sol[0][cx] + sol[0][rr]), " (= 1/p )")

=== 4.1 正交條件：|c1-c2|^2 = r1^2 + r2^2 ===
候選半徑 rho = sqrt(a**2 - 1)
驗算 a^2 - (1 + rho^2) = 0
兩半徑向量內積 = 0
切線夾角 = pi/2

=== 4.2 陷阱：圓心 (a,0)、半徑 sqrt(1+a^2) 的圓過哪些點？ ===
(x-a)^2+y^2 = 1+a^2  令 x=0 → y^2 = 1
→ 它通過單位圓上的兩個極點 (0, ±1)
在 (0,1) 處兩半徑向量內積 = 1 （正交應為 0）
它與單位圓的切線夾角 = acos(1/sqrt(a**2 + 1))  ; a=2 時 = 1.10714871779409 rad = 63.4349488229220 度

=== 4.3 通過圓盤內定點 P=(p,0)、且圓心在 x 軸上的測地線 ===
解: {cx: (p**2 + 1)/(2*p), rr: (1 - p**2)/(2*p)}
解: {cx: (p**2 + 1)/(2*p), rr: (p**2 - 1)/(2*p)}
圓心的另一個表達：半徑 rr 的測地線，其圓心距原點 = sqrt(1+rr^2) = sqrt(p**4 + 2*p**2 + 1)/(2*p)
與 x 軸的另一交點 = 1/p  (= 1/p )

from sympy import *
import random

xr, yr, xw, yw = symbols('xr yr xw yw', real=True)
r, t = symbols('r t', positive=True)

print("=== 5.1 關鍵恆等式：|1-w̄z|^2 = (1-|z|^2)(1-|w|^2) + |z-w|^2 ===")
lhs = (1 - (xw*xr + yw*yr))**2 + (yw*xr - xw*yr)**2
rhs = (1 - xr**2 - yr**2)*(1 - xw**2 - yw**2) + (xr-xw)**2 + (yr-yw)**2
print("lhs - rhs =", simplify(lhs - rhs))

print()
print("=== 5.2 三種距離公式等價：先用一個正數 u 把它們化開 ===")
u = symbols('u', positive=True)
x = u/sqrt(1 + u**2)
print("atanh(x) = asinh(u) 的關鍵：x/sqrt(1-x^2) - u =", simplify(x/sqrt(1 - x**2) - u))
print("cosh(2 asinh u) - (1 + 2 u^2) =",
      expand_trig(cosh(2*asinh(u)) - 1 - 2*sinh(asinh(u))**2))
print("（於是 atanh 形式 = asinh 形式 = arcosh 形式）")

print()
print("=== 5.3 精確有理數檢查：z=0, w=r = 1/2 ===")
rr = Rational(1, 2)
d1 = 2*atanh(rr)
d2 = 2*asinh(rr/sqrt(1 - rr**2))
d3 = acosh(1 + 2*rr**2/(1 - rr**2))
print("d1 =", d1, "= 2*atanh(1/2)")
print("d2 =", d2, "= 2*asinh(1/sqrt(3))")
print("d3 =", d3, "= acosh(5/3)")
print("d1.rewrite(log) - d2.rewrite(log) =", simplify(d1.rewrite(log) - d2.rewrite(log)))
print("d1.rewrite(log) - d3.rewrite(log) =", simplify(d1.rewrite(log) - d3.rewrite(log)))
print("三者共同的對數形式 =", simplify(d1.rewrite(log)))

print()
print("=== 5.4 沿半徑的距離：雙曲距離把圓盤撐開 ===")
print("沿半徑的度量：ds = 2|dz|/(1-|z|^2) ，沿 z = t 的原函数")
print("  d/dt log((1+t)/(1-t)) - 2/(1-t^2) =", simplify(diff(log((1+t)/(1-t)), t) - 2/(1-t**2)))
print("  log((1+r)/(1-r)) - 2*atanh(r) =",
      simplify(expand_log(log((1+r)/(1-r)), force=True) - (2*atanh(r)).rewrite(log)))
print()
print("   r        歐氏距離      龐加萊距離 2*atanh(r)    比值 d/r")
for rr in [0.1, 0.25, 0.5, 0.75, 0.9, 0.99]:
    d = 2*atanh(rr).evalf()
    print(f"  {rr:<8}{rr:<13}{d:<24.10f}{d/rr:.6f}")
print("  r -> 0 時 d/r ->", limit(2*atanh(r)/r, r, 0), " ; r -> 1 時 d/r ->",
      limit(2*atanh(r)/r, r, 1))

print()
print("=== 5.5 對稱性與三角不等式（數值實驗） ===")
def pd(z1, z2):
    return 2*atanh(abs(z1 - z2)/abs(1 - conjugate(z2)*z1))

pts = [0, 0.3+0.2*I, 0.7*I, -0.5+0.4*I, 0.9]
worst = max(abs(pd(A, B) - pd(B, A)) for A in pts for B in pts)
print("對稱性 d(z,w) = d(w,z)：最大誤差 =", worst)
random.seed(42)
bad = 0
for _ in range(3000):
    A = complex(random.uniform(-0.9, 0.9), random.uniform(-0.9, 0.9))
    B = complex(random.uniform(-0.9, 0.9), random.uniform(-0.9, 0.9))
    C = complex(random.uniform(-0.9, 0.9), random.uniform(-0.9, 0.9))
    if max(abs(A), abs(B), abs(C)) >= 1:
        continue
    if pd(A, C) > pd(A, B) + pd(B, C) + 1e-12:
        bad += 1
print("3000 組隨機點違反三角不等式的次數 =", bad)

=== 5.1 關鍵恆等式：|1-w̄z|^2 = (1-|z|^2)(1-|w|^2) + |z-w|^2 ===
lhs - rhs = 0

=== 5.2 三種距離公式等價：先用一個正數 u 把它們化開 ===
atanh(x) = asinh(u) 的關鍵：x/sqrt(1-x^2) - u = 0
cosh(2 asinh u) - (1 + 2 u^2) = 0
（於是 atanh 形式 = asinh 形式 = arcosh 形式）

=== 5.3 精確有理數檢查：z=0, w=r = 1/2 ===
d1 = 2*atanh(1/2) = 2*atanh(1/2)
d2 = 2*asinh(sqrt(3)/3) = 2*asinh(1/sqrt(3))
d3 = acosh(5/3) = acosh(5/3)
d1.rewrite(log) - d2.rewrite(log) = 0
d1.rewrite(log) - d3.rewrite(log) = 0
三者共同的對數形式 = log(3)

=== 5.4 沿半徑的距離：雙曲距離把圓盤撐開 ===
沿半徑的度量：ds = 2|dz|/(1-|z|^2) ，沿 z = t 的原函数
  d/dt log((1+t)/(1-t)) - 2/(1-t^2) = 0
  log((1+r)/(1-r)) - 2*atanh(r) = 0

   r        歐氏距離      龐加萊距離 2*atanh(r)    比值 d/r
  0.1     0.1          0.2006706955            2.006707
  0.25    0.25         0.5108256238            2.043302
  0.5     0.5          1.0986122887            2.197225
  0.75    0.75         1.9459101491            2.594547
  0.9     0.9          2.9444389792            3.271599
  0.99    0.99         5.2933048247            5.346773
  r -> 0 時 d/r -> 2  ; r -> 1 時 d/r -> oo

=== 5.5 對稱性與三角不等式（數值實驗） ===
對稱性 d(z,w) = d(w,z)：最大誤差 = 0
3000 組隨機點違反三角不等式的次數 = 0

import numpy as np
from sympy import *

print("=== 6.1 圓盤模型的面板度量是共形的 → 雙曲角 = 歐氏角 ===")
x, y = symbols('x y', real=True)
g = Matrix([[4/(1 - x**2 - y**2)**2, 0],
            [0, 4/(1 - x**2 - y**2)**2]])
print("g =", g)
print("g12 =", g[0, 1], "  g11 - g22 =", simplify(g[0, 0] - g[1, 1]))

# ------------------------------------------------------------------ 工具
def geodesic_circle(P, Q):
    """過 P、Q 且與單位圓正交的圓（圓心、半徑）"""
    M = np.array([P, Q], dtype=float)
    rhs = np.array([(1 + P@P)/2, (1 + Q@Q)/2])
    c = np.linalg.solve(M, rhs)
    return c, np.sqrt(c @ c - 1)

def tangent_at(P, c, rho, toward):
    r = (P - c)/rho
    t = np.array([-r[1], r[0]])
    return -t if t @ (toward - P) < 0 else t

def angle(u, v):
    u = u/np.linalg.norm(u); v = v/np.linalg.norm(v)
    return np.arccos(np.clip(u @ v, -1, 1))

def triangle_data(V):
    centers = [geodesic_circle(V[i], V[(i+1) % 3]) for i in range(3)]
    angs = []
    for i in range(3):
        P = V[i]
        cp, rp = centers[(i-1) % 3]
        cn, rn = centers[i]
        angs.append(angle(tangent_at(P, cp, rp, V[(i-1) % 3]),
                          tangent_at(P, cn, rn, V[(i+1) % 3])))
    return angs, centers

def area_grid(centers, N=2000):
    """在圓盤內對三角形做 Riemann 網格積分（面積元素 4 dx dy/(1-r^2)^2）"""
    xs = np.linspace(-1, 1, N+1); ys = np.linspace(-1, 1, N+1)
    X, Y = np.meshgrid(xs, ys, indexing='ij')
    dx = xs[1]-xs[0]; dy = ys[1]-ys[0]
    r2 = X**2 + Y**2
    mask = r2 < 1 - 1e-12
    for (c, r) in centers:
        mask &= ((X-c[0])**2 + (Y-c[1])**2) >= r**2
    dens = np.where(r2 < 1-1e-12, 4.0/np.maximum(1-r2, 1e-12)**2, 0.0)
    return dens[mask].sum()*dx*dy

def area_grid_right(A, B, C, c, rho, N=2000):
    """直角型三角形：兩條直邊 + 一條圓弧；區域 = 盒子 ∩ 圓弧外側"""
    xs = np.linspace(A[0], B[0], N+1); ys = np.linspace(C[1], A[1], N+1)
    X, Y = np.meshgrid(xs, ys, indexing='ij')
    dx = xs[1]-xs[0]; dy = ys[1]-ys[0]
    mask = ((X - c[0])**2 + (Y - c[1])**2) >= rho**2
    dens = 4.0/(1 - X**2 - Y**2)**2
    return dens[mask].sum()*dx*dy

# ------------------------------------------------------------------ 6.2
print()
print("=== 6.2 一個具體的雙曲三角形 A=(0,0), B=(0.5,0), C=(0,-0.5) ===")
A = np.array([0.0, 0.0]); B = np.array([0.5, 0.0]); C = np.array([0.0, -0.5])
c0, rho0 = geodesic_circle(B, C)
print("弧 BC 的圓心 =", c0, " 半徑 =", rho0)
angA = angle(B - A, C - A)                     # A 處兩邊是直徑，夾角 90 度
angB = angle(A - B, tangent_at(B, c0, rho0, C))
angC = angle(A - C, tangent_at(C, c0, rho0, B))
for nm, a_ in [("A", angA), ("B", angB), ("C", angC)]:
    print(f"  角 {nm} = {np.degrees(a_):.6f} 度  (= {a_:.10f} rad)")
defect = np.pi - (angA + angB + angC)
print("角虧 π-(A+B+C) =", defect, " =", np.degrees(defect), "度")

print()
print("=== 6.3 用網格積分直接算面積，與 Hurwitz 公式對照 ===")
for N in [500, 1000, 2000, 4000]:
    num = area_grid_right(A, B, C, c0, rho0, N)
    print(f"  N={N:<6} 網格面積 = {num:.10f}   Hurwitz = {defect:.10f}   差 = {abs(num-defect):.2e}")

# ------------------------------------------------------------------ 6.4
print()
print("=== 6.4 等邊三角形族：頂點在 Euclidean 半徑 r 的圓上 ===")
print("   r      角 A(度)     角虧=Hurwitz 面積    網格積分面積        差")
for r in [0.3, 0.5, 0.7, 0.8, 0.9, 0.95, 0.99]:
    th = np.array([0.0, 2*np.pi/3, 4*np.pi/3])
    V = np.column_stack([r*np.cos(th), r*np.sin(th)])
    angs, centers = triangle_data(V)
    defect = np.pi - sum(angs)
    num = area_grid(centers, 1200)
    print(f"  {r:<6}{np.degrees(angs[0]):<12.6f}{defect:<20.10f}{num:<20.10f}{abs(num-defect):.2e}")

# ------------------------------------------------------------------ 6.5
print()
print("=== 6.5 理想三角形（頂點在邊界上）：面積 = π ===")
th = np.array([0.0, 2*np.pi/3, 4*np.pi/3])
V = np.column_stack([np.cos(th), np.sin(th)])
angs, centers = triangle_data(V)
print("三條邊的圓心與半徑 =", [(tuple(np.round(c, 6)), round(r, 6)) for (c, r) in centers])
print("三個內角 =", angs, "（理想頂點：兩邊相切 → 角為 0）")
print("Hurwitz 面積 = π - 0 =", np.pi)

# 用極座標做數值積分（θ = s² 消去尖點的可積奇性）
def r_max(t):
    c = 2*np.cos(t - np.pi/3)
    return c - np.sqrt(np.maximum(c**2 - 1, 0.0))

def integrand(t):
    r = min(r_max(t), 1 - 1e-13)
    return 2.0/(1 - r**2) - 2.0

for N in [1000, 4000, 16000]:
    s = np.linspace(0, np.sqrt(np.pi/3), N+1)
    val = np.trapezoid([2*ss*integrand(ss**2) for ss in s], s)
    print(f"  N={N:<7} 極座標數值積分面積 = {6*val:.10f}   π = {np.pi:.10f}   差 = {6*val-np.pi:.2e}")

=== 6.1 圓盤模型的面板度量是共形的 → 雙曲角 = 歐氏角 ===
g = Matrix([[4/(-x**2 - y**2 + 1)**2, 0], [0, 4/(-x**2 - y**2 + 1)**2]])
g12 = 0   g11 - g22 = 0

=== 6.2 一個具體的雙曲三角形 A=(0,0), B=(0.5,0), C=(0,-0.5) ===
弧 BC 的圓心 = [ 1.25 -1.25]  半徑 = 1.4577379737113252
  角 A = 90.000000 度  (= 1.5707963268 rad)
  角 B = 30.963757 度  (= 0.5404195003 rad)
  角 C = 30.963757 度  (= 0.5404195003 rad)
角虧 π-(A+B+C) = 0.48995732625372845  = 28.072486935852968 度

=== 6.3 用網格積分直接算面積，與 Hurwitz 公式對照 ===
  N=500    網格面積 = 0.4924496802   Hurwitz = 0.4899573263   差 = 2.49e-03
  N=1000   網格面積 = 0.4911888904   Hurwitz = 0.4899573263   差 = 1.23e-03
  N=2000   網格面積 = 0.4905653005   Hurwitz = 0.4899573263   差 = 6.08e-04
  N=4000   網格面積 = 0.4902601377   Hurwitz = 0.4899573263   差 = 3.03e-04

=== 6.4 等邊三角形族：頂點在 Euclidean 半徑 r 的圓上 ===
   r      角 A(度)     角虧=Hurwitz 面積    網格積分面積        差
  0.3   51.468879   0.4466884268        0.4467728554        8.44e-05
  0.5   38.213211   1.1407536201        1.1409671109        2.13e-04
  0.7   22.357125   1.9709763351        1.9712668237        2.90e-04
  0.8   14.445796   2.3852125290        2.3853185920        1.06e-04
  0.9   6.936432    2.7784019173        2.7793066645        9.05e-04
  0.95  3.389575    2.9641149156        2.9732292002        9.11e-03
  0.99  0.664895    3.1067788318        3.2167591677        1.10e-01

=== 6.5 理想三角形（頂點在邊界上）：面積 = π ===
三條邊的圓心與半徑 = [((np.float64(1.0), np.float64(1.732051)), np.float64(1.732051)), ((np.float64(-2.0), np.float64(0.0)), np.float64(1.732051)), ((np.float64(1.0), np.float64(-1.732051)), np.float64(1.732051))]
三個內角 = [np.float64(0.0), np.float64(0.0), np.float64(0.0)] （理想頂點：兩邊相切 → 角為 0）
Hurwitz 面積 = π - 0 = 3.141592653589793
  N=1000    極座標數值積分面積 = 3.1382949545   π = 3.1415926536   差 = -3.30e-03
  N=4000    極座標數值積分面積 = 3.1407680021   π = 3.1415926536   差 = -8.25e-04
  N=16000   極座標數值積分面積 = 3.1413864765   π = 3.1415926536   差 = -2.06e-04

import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Arc as MplArc, Circle as MplCircle
from sympy import *
from sympy.plotting import plot_parametric

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'PingFang SC', 'Heiti TC', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

t = symbols('t', real=True)

def arc_angles(cx, cy, rho):
    """圓心 (cx,cy)、半徑 rho 的圓落在單位圓內的那一段弧的角度範圍"""
    phi0 = np.arctan2(cy, cx)
    psi = np.arccos(-rho/np.hypot(cx, cy))         # 交點相對圓心的角度偏角
    return np.degrees(phi0 + psi), np.degrees(phi0 + 2*np.pi - psi)

def draw_arc(ax, cx, cy, rho, **kw):
    t1, t2 = arc_angles(cx, cy, rho)
    ax.add_patch(MplArc((cx, cy), 2*rho, 2*rho, theta1=t1, theta2=t2, **kw))

def geodesic_circle(P, Q):
    M = np.array([P, Q], dtype=float)
    rhs = np.array([(1 + P@P)/2, (1 + Q@Q)/2])
    c = np.linalg.solve(M, rhs)
    return c, np.sqrt(c @ c - 1)

# ------------------------------------------------- 圖 1：圓盤與測地線
p = plot_parametric((cos(t), sin(t)), (t, 0, 2*pi))     # sympy 畫邊界
fig = p.fig                                          # 取用 sympy 產生的 figure/axes
ax = p.ax
ax.add_patch(MplCircle((0, 0), 1, fill=False, ls='--', lw=1.2, color='0.45'))
for (cx, cy) in [(2, 0), (0, 2), (-2, 0), (1.5, 0.5), (0.5, -1.5)]:
    draw_arc(ax, cx, cy, np.sqrt(cx**2 + cy**2 - 1), color='C0', lw=1.7)
ax.plot([-1, 1], [0, 0], color='C3', lw=1.8)
ax.plot([0, 0], [-1, 1], color='C3', lw=1.8)
ax.set_xlim(-1.15, 1.15); ax.set_ylim(-1.15, 1.15)
ax.set_aspect('equal'); ax.grid(alpha=.25)
ax.set_title("龐加萊圓盤：藍＝與邊界正交的圓弧（測地線），\n紅＝直徑（退化情形），灰虛線＝邊界＝無窮遠")
fig.savefig("fig1_disk.png", dpi=140, bbox_inches='tight')
print("fig1_disk.png 完成")

# ------------------------------------------------- 圖 2：雙曲三角形
A = np.array([0.0, 0.0]); B = np.array([0.5, 0.0]); C = np.array([0.0, -0.5])
c0, rho0 = geodesic_circle(B, C)
fig, ax = plt.subplots(figsize=(6.2, 6.2))
ax.add_patch(MplCircle((0, 0), 1, fill=False, ls='--', lw=1.2, color='0.45'))
ax.plot([A[0], B[0]], [A[1], B[1]], 'C3', lw=2.2)
ax.plot([A[0], C[0]], [A[1], C[1]], 'C3', lw=2.2)
draw_arc(ax, c0[0], c0[1], rho0, color='C0', lw=2.2)
ax.plot(*zip(A, B, C, A), 'ko', ms=5)
for P, lab in [(A, "A"), (B, "B"), (C, "C")]:
    ax.annotate(lab, P, textcoords="offset points", xytext=(6, 6), fontsize=12)
ax.set_xlim(-1.1, 1.1); ax.set_ylim(-1.1, 1.1); ax.set_aspect('equal'); ax.grid(alpha=.25)
ax.set_title("雙曲三角形：紅＝直徑，藍＝正交圓弧；角要用『切線』量")
fig.savefig("fig2_triangle.png", dpi=140, bbox_inches='tight')
print("fig2_triangle.png 完成")

# ------------------------------------------------- 圖 3：莫比烏斯變換
fig, ax = plt.subplots(figsize=(6.4, 4.4))
for (cx, cy) in [(2, 0), (0, 2), (-2, 0)]:
    rho = np.sqrt(cx**2 + cy**2 - 1)
    t1, t2 = arc_angles(cx, cy, rho)
    tt = np.radians(np.linspace(t1, t2, 500))
    zs = (cx + rho*np.cos(tt)) + 1j*(cy + rho*np.sin(tt))     # 圓弧上的點 z
    ws = 1j*(1 + zs)/(1 - zs)                                # 莫比烏斯變換 w = i(1+z)/(1-z)
    ax.plot(np.real(ws), np.imag(ws), 'C0', lw=2)
ax.plot([0, 0], [0, 4], 'C3', lw=2)                  # 直徑 → 虛軸
ax.axhline(0, color='0.45', ls='--')
ax.set_xlim(-6, 6); ax.set_ylim(-0.4, 4.2)
ax.set_title("莫比烏斯變換 $w=i\\frac{1+z}{1-z}$：圓弧 → 半圓，直徑 → 垂直線")
fig.savefig("fig3_mobius.png", dpi=140, bbox_inches='tight')
print("fig3_mobius.png 完成")

# ------------------------------------------------- 圖 4：三種幾何的距離
fig, ax = plt.subplots(figsize=(6.2, 4.2))
rs = np.linspace(0.001, 0.999, 400)
ax.plot(rs, 2*np.arctanh(rs), 'C0', lw=2, label="雙曲 K=-1：$2\\operatorname{artanh}(r)$")
ax.plot(rs, rs, 'k--', lw=1.5, label="歐氏 K=0：$r$")
ax.set_xlabel("從圓心出發的 Euclidean 半徑 $r$")
ax.set_ylabel("距離")
ax.set_title("同一條半徑上，雙曲距離把圓盤「撐開」")
ax.legend(); ax.grid(alpha=.3)
fig.savefig("fig4_distance.png", dpi=140, bbox_inches='tight')
print("fig4_distance.png 完成")

fig1_disk.png 完成
fig2_triangle.png 完成
fig3_mobius.png 完成
fig4_distance.png 完成

from sympy import *

a = symbols('a', positive=True)          # 等邊三角形的邊長（K = ±1 時單位半徑）

print("=== 8.1 同樣的邊長 a = 1，在三種幾何裡給出不同的角與面積 ===")
A_eu = pi/3
A_hyp = acos(cosh(a)/(cosh(a) + 1))       # 雙曲：cos A = cosh a/(cosh a + 1)
A_sph = acos(cos(a)/(1 + cos(a)))         # 球面：cos A = cos a/(1 + cos a)

for name, A in [("歐氏 K=0", A_eu), ("雙曲 K=-1", A_hyp), ("球面 K=+1", A_sph)]:
    ang = A.subs(a, 1)
    if name.startswith("雙曲"):
        area = pi - 3*ang
    elif name.startswith("球面"):
        area = 3*ang - pi
    else:
        area = sqrt(3)/4
    print(f"  {name:<10} 內角 A = {N(deg(ang)):.6f} 度    面積 = {N(area):.10f}")

print()
print("=== 8.2 反過來：角固定為 60-60-60 時，邊長是多少？ ===")
ch = symbols('ch', positive=True)
eq_h = Eq(cos(pi/3), -cos(pi/3)*cos(pi/3) + sin(pi/3)*sin(pi/3)*ch)
print("雙曲：cos A = -cos B cos C + sin B sin C cosh a  =>", solve(eq_h, ch))
print("      ⇒ cosh a =", solve(eq_h, ch)[0], " ⇒ a = acosh(1) =", acosh(1), "（退化成一點！）")
cs = symbols('cs', real=True)
eq_s = Eq(cos(pi/3), -cos(pi/3)*cos(pi/3) + sin(pi/3)*sin(pi/3)*cs)
print("球面：cos A = -cos B cos C + sin B sin C cos a   =>", solve(eq_s, cs))
print("      ⇒ cos a =", solve(eq_s, cs)[0], " ⇒ a = 0（也退化！）")
print("歐氏：a 任意正數皆可 ⇒ 存在相似但不全等的等邊三角形 ✓")

=== 8.1 同樣的邊長 a = 1，在三種幾何裡給出不同的角與面積 ===
  歐氏 K=0     內角 A = 60.000000 度    面積 = 0.4330127019
  雙曲 K=-1    內角 A = 52.643240 度    面積 = 0.3851990371
  球面 K=+1    內角 A = 69.465165 度    面積 = 0.4955948957

=== 8.2 反過來：角固定為 60-60-60 時，邊長是多少？ ===
雙曲：cos A = -cos B cos C + sin B sin C cosh a  => [1]
      ⇒ cosh a = 1  ⇒ a = acosh(1) = 0 （退化成一點！）
球面：cos A = -cos B cos C + sin B sin C cos a   => [1]
      ⇒ cos a = 1  ⇒ a = 0（也退化！）
歐氏：a 任意正數皆可 ⇒ 存在相似但不全等的等邊三角形 ✓