#!/usr/bin/env python3
"""
微积分求解地球表面积演示
Calculating Earth's Surface Area using Calculus
"""

import math

# 地球参数
EARTH_RADIUS = 6371000  # 米 (平均半径)

def spherical_surface_element(theta, phi, r):
    """
    球面微元面积: dS = r² sin(θ) dθ dφ
    θ: 极角 (0 到 π)
    φ: 方位角 (0 到 2π)
    r: 半径
    """
    return r**2 * math.sin(theta)

def numerical_integration_2d(func, theta_range, phi_range, n_theta=1000, n_phi=1000):
    """二重积分数值计算"""
    theta_min, theta_max = theta_range
    phi_min, phi_max = phi_range
    
    d_theta = (theta_max - theta_min) / n_theta
    d_phi = (phi_max - phi_min) / n_phi
    
    total = 0
    for i in range(n_theta):
        theta = theta_min + (i + 0.5) * d_theta
        for j in range(n_phi):
            phi = phi_min + (j + 0.5) * d_phi
            total += func(theta, phi, EARTH_RADIUS) * d_theta * d_phi
    
    return total

def analytical_solution(r):
    """
    解析解推导:
    S = ∫∫ r² sin(θ) dθ dφ
    = r² ∫₀^π sin(θ) dθ ∫₀^2π dφ
    = r² × [-cos(θ)]₀^π × [φ]₀^2π
    = r² × [-cos(π) + cos(0)] × [2π - 0]
    = r² × [1 + 1] × 2π
    = 4πr²
    """
    return 4 * math.pi * r**2

def main():
    print("🌍 微积分求解地球表面积")
    print("=" * 40)
    
    print(f"🔧 参数设置:")
    print(f"  地球半径 R = {EARTH_RADIUS:,} m = {EARTH_RADIUS/1000:.0f} km")
    
    # 积分区域
    theta_range = (0, math.pi)      # 极角: 0 到 π
    phi_range = (0, 2*math.pi)      # 方位角: 0 到 2π
    
    print(f"\n📐 积分区域:")
    print(f"  极角 θ: {theta_range[0]:.1f} 到 {theta_range[1]:.2f} (0 到 π)")
    print(f"  方位角 φ: {phi_range[0]:.1f} 到 {phi_range[1]:.2f} (0 到 2π)")
    
    # 数值积分计算
    print(f"\n🧮 积分计算:")
    numerical_result = numerical_integration_2d(
        spherical_surface_element, theta_range, phi_range, 500, 500
    )
    
    # 解析解
    analytical_result = analytical_solution(EARTH_RADIUS)
    
    print(f"  数值积分: {numerical_result:.2e} m²")
    print(f"  解析解:   {analytical_result:.2e} m²")
    print(f"  相对误差: {abs(numerical_result - analytical_result) / analytical_result * 100:.6f}%")
    
    # 单位转换
    numerical_km2 = numerical_result / 1e6
    analytical_km2 = analytical_result / 1e6
    
    print(f"\n🗺️  表面积结果:")
    print(f"  数值积分: {numerical_km2:,.0f} km²")
    print(f"  解析解:   {analytical_km2:,.0f} km²")
    
    # 微积分推导过程
    print(f"\n📚 微积分推导过程:")
    print("1. 球坐标系微元面积: dS = r² sin(θ) dθ dφ")
    print("2. 表面积积分: S = ∫∫ dS = ∫₀^π ∫₀^2π r² sin(θ) dθ dφ")
    print("3. 分离变量: S = r² ∫₀^π sin(θ) dθ × ∫₀^2π dφ")
    print("4. 计算内积分: ∫₀^π sin(θ) dθ = [-cos(θ)]₀^π = 2")
    print("5. 计算外积分: ∫₀^2π dφ = 2π")
    print("6. 最终结果: S = r² × 2 × 2π = 4πr²")
    
    # 验证步骤
    print(f"\n🔍 分步验证:")
    inner_integral = 2  # ∫₀^π sin(θ) dθ = 2
    outer_integral = 2 * math.pi  # ∫₀^2π dφ = 2π
    step_result = EARTH_RADIUS**2 * inner_integral * outer_integral
    
    print(f"  内积分 ∫sin(θ)dθ = {inner_integral}")
    print(f"  外积分 ∫dφ = {outer_integral:.4f}")
    print(f"  r² × 内积分 × 外积分 = {step_result:.2e} m²")
    
    # 与已知数据对比
    known_earth_surface = 510.1e6  # km² (已知地球表面积)
    calculated_km2 = analytical_result / 1e6
    
    print(f"\n🌐 与实际数据对比:")
    print(f"  计算结果: {calculated_km2:,.1f} km²")
    print(f"  实际数据: {known_earth_surface:,.1f} km²")
    print(f"  差异: {abs(calculated_km2 - known_earth_surface) / known_earth_surface * 100:.1f}%")
    print("  (差异主要由于地球非完美球体)")

if __name__ == "__main__":
    main()