run_demo.py

#!/usr/bin/env python3
"""
Main Demo Script for Multipole Expansion

This is the primary entry point for demonstrating all features of the
multipole expansion package.

Usage:
    uv run run_demo.py
"""

from multipole_expansion import MultipoleExpansion


def print_header(title):
    """Print a formatted header."""
    print("\n" + "=" * 80)
    print(f"  {title}")
    print("=" * 80 + "\n")


def main():
    """Main demonstration function."""

    print_header("MULTIPOLE EXPANSION WITH SYMBOLICA")

    print("""
This package implements generic multipole moment calculations using pure Symbolica.

The multipole expansion decomposes the potential
    φ(x) = 1/|x - x_a|
    
into a series
    φ(x) = φ^(0) + φ^(1) + φ^(2) + ...

where each term captures different geometric features of the source distribution.

We implement TWO equivalent formulations:
    1. Taylor expansion: φ^(n) = ((-1)^n/n!) x_a^{i₁}...x_a^{iₙ} ∂^n/∂x^{i₁}...∂x^{iₙ}(1/r)
    2. Q tensor form: φ^(n) = (1/n!) Q^{i₁...iₙ} n^{i₁}...n^{iₙ} / r^{n+1}

And verify they are equivalent.

See the phi_from_Q function in `multipole_expansion/multipole_moments.py` for the 2nd formulation,
and the phi_n function in `multipole_expansion/taylor_expansion.py` for the 1st formulation.

Detailed explanation can be found at `report/main.pdf`.

The output formulae below are generated by code, not hard-wired text. See run_demo.py for details.
""")

    # Initialize
    print("Initializing multipole expansion framework...")
    mp = MultipoleExpansion(max_order=2)
    print("✓ Initialization complete\n")

    # Run complete computation and verification
    mp.compute_all(max_order=3, run_verification=True)

    # Additional detailed examples
    print_header("DETAILED EXAMPLE: Quadrupole Term")

    print("The quadrupole (n=2) is particularly interesting.")
    print("Let's examine it in detail:\n")

    from symbolica import S

    print("1. Cartesian Q tensor:")
    i, j = S("i"), S("j")
    Q_2 = mp.compute_Q_tensor(2, [i, j])
    print(f"   Q^{{ij}} = {Q_2}")
    print("   This is symmetric: Q^{ij} = Q^{ji}")
    print("   And traceless: Q^{ii} = 0")

    print("\n2. Specific components (using i,j = x,y,z):")
    for idx1 in ["x", "y", "z"]:
        for idx2 in ["x", "y", "z"]:
            if idx1 <= idx2:  # Only upper triangle (symmetric)
                Q_comp = mp.compute_Q_tensor(2, [S(idx1), S(idx2)])
                label = f"Q^{{{idx1}{idx2}}}"
                print(f"   {label:8s} = {Q_comp}")

    print("\n3. Quadrupole potential:")
    phi_2 = mp.taylor_term(2)
    print(f"   φ^(2) = {phi_2}")
    print("   This is (3(x_a·n)² - x_a²)/(2r³)")

    print("\n4. Spherical representation:")
    print("   Converting to spherical harmonics basis:")
    for m in range(-2, 3):
        try:
            q_2m = mp.spherical_moment(2, m)
            print(f"   q_{{2,{m:+d}}} = {q_2m}")
        except Exception as e:
            print(f"   q_{{2,{m:+d}}} : {e}")

    print_header("SUMMARY")

    print("""
✓ Successfully implemented multipole expansion in pure Symbolica
✓ Verified all theoretical properties
✓ Demonstrated both Cartesian and spherical representations

This implementation can be used for:
  - Electrostatics problems
  - Gravitational potentials  
  - Quantum mechanics (angular momentum)
  - Molecular physics
  - Any system with 1/r potentials!

For more details, see:
  - README.md for theory
  - multipole_expansion/ for implementation
  - examples/ for more demonstrations
    """)

    print("=" * 80)
    print("Demo complete! Thank you for using the Multipole Expansion package.")
    print("=" * 80)


if __name__ == "__main__":
    try:
        main()
    except KeyboardInterrupt:
        print("\n\nDemo interrupted by user.")
    except Exception as e:
        print(f"\n\nError during demo: {e}")
        import traceback

        traceback.print_exc()