#!/usr/bin/env python3
"""
地球表面积微积分计算 - 多种方法演示
Earth Surface Area Calculation - Multiple Calculus Methods
"""

import math

EARTH_RADIUS = 6371000  # 米

def method1_spherical_coordinates():
    """方法1: 球坐标系二重积分"""
    print("📐 方法1: 球坐标系二重积分")
    print("   dS = r² sin(θ) dθ dφ")
    print("   S = ∫₀^π ∫₀^2π r² sin(θ) dθ dφ")
    
    # 解析计算
    r = EARTH_RADIUS
    inner_integral = 2  # ∫₀^π sin(θ) dθ = [-cos(θ)]₀^π = 2
    outer_integral = 2 * math.pi  # ∫₀^2π dφ = 2π
    result = r**2 * inner_integral * outer_integral
    
    print(f"   结果: {result:.2e} m² = {result/1e6:,.0f} km²")
    return result

def method2_revolution_surface():
    """方法2: 旋转曲面积分"""
    print("\n🔄 方法2: 旋转曲面积分")
    print("   半圆 x² + z² = r² 绕z轴旋转")
    print("   x = r sin(θ), z = r cos(θ)")
    print("   dS = 2πx √(1 + (dx/dz)²) dz")
    
    # 参数方程: x = r sin(θ), z = r cos(θ)
    # dx/dθ = r cos(θ), dz/dθ = -r sin(θ)
    # ds = r dθ (弧长微元)
    # 旋转表面积: S = ∫₀^π 2π × r sin(θ) × r dθ
    
    r = EARTH_RADIUS
    # S = 2πr² ∫₀^π sin(θ) dθ = 2πr² × 2 = 4πr²
    result = 4 * math.pi * r**2
    
    print(f"   S = 2πr² ∫₀^π sin(θ) dθ = 4πr²")
    print(f"   结果: {result:.2e} m² = {result/1e6:,.0f} km²")
    return result

def method3_parametric_surface():
    """方法3: 参数曲面积分"""
    print("\n📊 方法3: 参数曲面积分")
    print("   r⃗(θ,φ) = (r sin(θ)cos(φ), r sin(θ)sin(φ), r cos(θ))")
    print("   dS = |∂r⃗/∂θ × ∂r⃗/∂φ| dθ dφ")
    
    # 偏导数计算
    print("   ∂r⃗/∂θ = (r cos(θ)cos(φ), r cos(θ)sin(φ), -r sin(θ))")
    print("   ∂r⃗/∂φ = (-r sin(θ)sin(φ), r sin(θ)cos(φ), 0)")
    
    # 叉积模长 = r² sin(θ)
    print("   |∂r⃗/∂θ × ∂r⃗/∂φ| = r² sin(θ)")
    
    r = EARTH_RADIUS
    result = 4 * math.pi * r**2
    
    print(f"   S = ∫₀^π ∫₀^2π r² sin(θ) dθ dφ = 4πr²")
    print(f"   结果: {result:.2e} m² = {result/1e6:,.0f} km²")
    return result

def method4_cartesian_coordinates():
    """方法4: 直角坐标系积分"""
    print("\n📏 方法4: 直角坐标系积分")
    print("   球面方程: x² + y² + z² = r²")
    print("   上半球: z = √(r² - x² - y²)")
    print("   dS = √(1 + (∂z/∂x)² + (∂z/∂y)²) dx dy")
    
    # 偏导数: ∂z/∂x = -x/z, ∂z/∂y = -y/z
    # 1 + (∂z/∂x)² + (∂z/∂y)² = 1 + x²/z² + y²/z² = (x² + y² + z²)/z² = r²/z²
    # dS = (r/z) dx dy = r/√(r² - x² - y²) dx dy
    
    print("   √(1 + (∂z/∂x)² + (∂z/∂y)²) = r/√(r² - x² - y²)")
    print("   上半球面积 = ∫∫ r/√(r² - x² - y²) dx dy")
    print("   转换为极坐标: x = ρcos(φ), y = ρsin(φ)")
    print("   = ∫₀^2π ∫₀^r r/√(r² - ρ²) × ρ dρ dφ")
    
    # 内积分: ∫₀^r rρ/√(r² - ρ²) dρ
    # 设 u = r² - ρ², du = -2ρ dρ
    # = -r/2 ∫ᵣ²^0 u^(-1/2) du = -r/2 × 2√u |ᵣ²^0 = r²
    
    r = EARTH_RADIUS
    hemisphere_area = 2 * math.pi * r**2
    total_area = 2 * hemisphere_area  # 上下两个半球
    
    print(f"   上半球面积 = 2πr²")
    print(f"   总表面积 = 2 × 2πr² = 4πr²")
    print(f"   结果: {total_area:.2e} m² = {total_area/1e6:,.0f} km²")
    return total_area

def numerical_verification():
    """数值验证"""
    print("\n🔢 数值积分验证:")
    
    def integrand(theta, phi):
        return EARTH_RADIUS**2 * math.sin(theta)
    
    # 蒙特卡洛积分
    import random
    n_samples = 1000000
    total = 0
    
    for _ in range(n_samples):
        theta = random.uniform(0, math.pi)
        phi = random.uniform(0, 2*math.pi)
        total += integrand(theta, phi)
    
    # 积分区域体积 = π × 2π = 2π²
    monte_carlo_result = total / n_samples * math.pi * 2 * math.pi
    analytical_result = 4 * math.pi * EARTH_RADIUS**2
    
    print(f"   蒙特卡洛积分 ({n_samples:,} 样本): {monte_carlo_result:.2e} m²")
    print(f"   解析解: {analytical_result:.2e} m²")
    print(f"   相对误差: {abs(monte_carlo_result - analytical_result) / analytical_result * 100:.2f}%")

def main():
    print("🌍 地球表面积微积分计算 - 多种方法演示")
    print("=" * 50)
    print(f"地球半径: {EARTH_RADIUS:,} m = {EARTH_RADIUS/1000} km")
    
    # 四种不同的微积分方法
    result1 = method1_spherical_coordinates()
    result2 = method2_revolution_surface()
    result3 = method3_parametric_surface()
    result4 = method4_cartesian_coordinates()
    
    # 数值验证
    numerical_verification()
    
    # 总结
    print(f"\n📋 方法总结:")
    print(f"   所有方法都得到相同结果: 4πr² = {4*math.pi*EARTH_RADIUS**2:.2e} m²")
    print(f"   地球表面积: {4*math.pi*EARTH_RADIUS**2/1e6:,.0f} km²")
    
    print(f"\n🎯 微积分核心概念:")
    print("   1. 微元选择: 根据坐标系选择合适的面积微元")
    print("   2. 积分区域: 确定参数的取值范围")
    print("   3. 坐标变换: 不同坐标系间的转换")
    print("   4. 分离变量: 将二重积分分解为两个单积分")
    print("   5. 几何直观: 理解积分的几何意义")

if __name__ == "__main__":
    main()