//! Example demonstrating the use of Lebedev quadrature for spherical integration
//!
//! This example shows various use cases for Lebedev quadrature:
//! - Basic integration over a unit sphere
//! - Integration of spherical harmonics
//! - Computation of moments of distributions on the sphere
//! - Application to physics problems (electrostatic potential)

use scirs2_integrate::lebedev::{lebedev_integrate, lebedev_rule, LebedevOrder};
use std::f64::consts::PI;

#[allow(dead_code)]
fn main() {
    println!("Lebedev Quadrature Examples");
    println!("==========================\n");

    // Example 1: Basic Lebedev rule information
    println!("Example 1: Basic Lebedev Rule Information\n");
    println!("Available Lebedev quadrature orders:");

    let orders = [
        LebedevOrder::Order6,
        LebedevOrder::Order14,
        LebedevOrder::Order26,
        LebedevOrder::Order38,
        LebedevOrder::Order50,
        LebedevOrder::Order74,
        LebedevOrder::Order86,
        LebedevOrder::Order110,
    ];

    for &order in &orders {
        print!("  {:?}: {} points", order, order.num_points());

        // Try to generate the rule to see if it's implemented
        match lebedev_rule::<f64>(order) {
            Ok(rule) => println!(" (implemented, degree {})", rule.degree),
            Err(_) => println!(" (not yet implemented)"),
        }
    }

    // Example 2: Integrating functions over the unit sphere
    println!("\nExample 2: Integrating Functions Over the Unit Sphere\n");

    // Test with constant function f(x,y,z) = 1
    // The integral over the unit sphere should be 4π (the surface area)
    let constant_result: f64 = lebedev_integrate(|___| 1.0, LebedevOrder::Order14).unwrap();
    println!("Integrating f(x,y,z) = 1:");
    println!("  Result: {constant_result:.10}");
    println!("  Expected: {:.10} (4π)", 4.0 * PI);
    println!("  Error: {:.10e}", (constant_result - 4.0 * PI).abs());

    // Test with a function dependent on coordinates
    // Here we use f(x,y,z) = x^2, which should integrate to 4π/3
    let x2_result: f64 = lebedev_integrate(|x, y, _z| x * x, LebedevOrder::Order14).unwrap();
    println!("\nIntegrating f(x,y,z) = x²:");
    println!("  Result: {x2_result:.10}");
    println!("  Expected: {:.10} (4π/3)", 4.0 * PI / 3.0);
    println!("  Error: {:.10e}", (x2_result - 4.0 * PI / 3.0).abs());

    // Test with a more complex function
    // f(x,y,z) = x^2 * y^2 * z^2, which should integrate to 4π/15
    let xyz_result: f64 =
        lebedev_integrate(|x, y, z| x * x * y * y * z * z, LebedevOrder::Order26).unwrap();
    println!("\nIntegrating f(x,y,z) = x² * y² * z²:");
    println!("  Result: {xyz_result:.10}");
    println!("  Expected: {:.10} (4π/15)", 4.0 * PI / 15.0);
    println!("  Error: {:.10e}", (xyz_result - 4.0 * PI / 15.0).abs());

    // Example 3: Spherical harmonics integration
    println!("\nExample 3: Spherical Harmonics Integration\n");

    // Spherical harmonics should integrate to zero over the sphere
    // due to orthogonality conditions

    // Y₁₀ ∝ z
    let y10_result: f64 = lebedev_integrate(|__, z| z, LebedevOrder::Order14).unwrap();
    println!("Integrating Y₁₀ ∝ z:");
    println!("  Result: {y10_result:.10e}");
    println!("  Expected: 0 (by orthogonality)");

    // Y₂₀ ∝ (3z² - 1)
    let y20_result: f64 =
        lebedev_integrate(|__, z| 3.0 * z * z - 1.0, LebedevOrder::Order14).unwrap();
    println!("\nIntegrating Y₂₀ ∝ (3z² - 1):");
    println!("  Result: {y20_result:.10e}");
    println!("  Expected: 0 (by orthogonality)");

    // Y₂₂ ∝ (x² - y²)
    let y22_result: f64 =
        lebedev_integrate(|x, y, _z| x * x - y * y, LebedevOrder::Order14).unwrap();
    println!("\nIntegrating Y₂₂ ∝ (x² - y²):");
    println!("  Result: {y22_result:.10e}");
    println!("  Expected: 0 (by orthogonality)");

    // Example 4: Comparison of different Lebedev orders
    println!("\nExample 4: Comparing Different Lebedev Orders\n");

    // Define a function with known integral (x^4 + y^4 + z^4 = 3/5 * (x² + y² + z²)² on the sphere)
    // The integral of (x² + y² + z²)² = 1² = 1 for points on the unit sphere
    // So the expected result is 3/5 * 4π
    let test_func = |x: f64, y: f64, z: f64| x.powi(4) + y.powi(4) + z.powi(4);
    let expected = 3.0 / 5.0 * 4.0 * PI;

    println!("Integrating f(x,y,z) = x⁴ + y⁴ + z⁴:");
    println!("Expected result: {expected:.10} (3/5 * 4π)");

    for &order in &[
        LebedevOrder::Order6,
        LebedevOrder::Order14,
        LebedevOrder::Order26,
        LebedevOrder::Order38,
        LebedevOrder::Order50,
    ] {
        // Skip orders that aren't implemented
        if let Ok(result) = lebedev_integrate::<f64>(test_func, order) {
            let error = (result - expected).abs();
            println!(
                "  {:?} ({} points): {:.10} (error: {:.10e})",
                order,
                order.num_points(),
                result,
                error
            );
        }
    }

    // Example 5: Physical application - Electrostatic potential
    println!("\nExample 5: Physical Application - Electrostatic Potential\n");

    println!("Calculating the electrostatic potential from a spherical charge distribution");

    // Simulate a point outside a uniformly charged sphere
    // The potential outside the sphere is equivalent to a point charge at the center
    // We'll place our observation point at (0, 0, 2) - twice the sphere radius away

    let observation_point = [0.0, 0.0, 2.0];

    // Define our charge distribution - uniform over the sphere
    let charge_density = 1.0 / (4.0 * PI); // Normalized to total charge of 1

    // Calculate potential using Lebedev quadrature
    // V(x') = ∫[ρ(x) / |x - x'|] dΩ
    let potential: f64 = lebedev_integrate(
        |x, y, z| {
            let dx: f64 = x - observation_point[0];
            let dy: f64 = y - observation_point[1];
            let dz: f64 = z - observation_point[2];
            let distance: f64 = (dx * dx + dy * dy + dz * dz).sqrt();

            charge_density / distance
        },
        LebedevOrder::Order26,
    )
    .unwrap();

    // For a unit sphere with unit total charge, the expected potential at distance 2
    // is simply 1/2 (from Coulomb's law: V = q/r)
    println!(
        "  Observation point: ({}, {}, {})",
        observation_point[0], observation_point[1], observation_point[2]
    );
    println!("  Calculated potential: {potential:.10}");
    println!("  Expected potential: {:.10} (1/2)", 0.5);
    println!(
        "  Relative error: {:.10e}",
        ((potential - 0.5).abs() / 0.5) as f64
    );
}
