scipy
SciPy is an open-source Python scientific computing library built on top of NumPy arrays, providing a wide range of numerical algorithms and tools. It is widely used in cislunar orbit design, dynamics analysis, signal processing, and statistical data analysis.
Main Modules
| Module | Function | Cislunar Application |
|---|---|---|
scipy.integrate | Numerical integration, ODE solvers | Orbit integration, dynamics propagation, state transition matrix |
scipy.optimize | Optimization algorithms, nonlinear & least-squares solvers | Orbit design, parameter optimization, Lambert problem solving |
scipy.linalg | Linear algebra operations, matrix decompositions | Matrix computation, eigenvalue analysis, covariance propagation |
scipy.interpolate | 1D/ND interpolation and smoothing splines | Data interpolation, ephemeris smoothing, observation gap filling |
scipy.special | Special mathematical functions (spherical harmonics, elliptic integrals, Bessel, etc.) | Spherical harmonic expansions, gravity field modeling |
scipy.fft | Fast Fourier Transform (FFT) and related transforms | Orbit spectral analysis, period detection, signal filtering |
scipy.signal | Signal processing (filtering, spectral analysis, windowing) | Range/range-rate signal processing, filter design |
scipy.stats | Statistical distributions, hypothesis tests, descriptive statistics | Orbit determination error analysis, Monte Carlo simulation |
scipy.spatial | Spatial data structures (KDTree, convex hull, Voronoi) | Star pattern matching, spatial proximity search, coverage analysis |
scipy.sparse | Sparse matrices and sparse linear algebra | Large-scale orbit network analysis, finite element modeling |
scipy.differentiate | Finite-difference differentiation tools | Sensitivity analysis, gradient computation |
scipy.cluster | Clustering algorithms | Orbit classification, space object clustering |
Technical Features
- Languages: Python 3.12+ for the primary API; core performance-critical code written in C (C17), C++ (C++17), Cython (≥3.2.0), and Fortran
- Build system: Meson (≥1.5.0) + meson-python, supporting Debug/Release/ASan build modes
- BLAS/LAPACK: Supports OpenBLAS (default), MKL, Accelerate (macOS), and others for high-performance linear algebra
- Array interface: Full NumPy 2.0+ support with ongoing Array API standard compatibility
- Version: Current stable release is 1.17+
Installation
# Standard installation
pip install scipy
# Development version (with test dependencies)
pip install -e ".[dev,test]"
# Via pixi (recommended for SciPy development)
pixi install
Example Applications
Orbit Integration (CR3BP)
import numpy as np
from scipy.integrate import solve_ivp
def orbital_dynamics(t, state, mu):
"""Circular Restricted Three-Body Problem (CR3BP) dynamics"""
x, y, z, vx, vy, vz = state
r1 = np.sqrt((x + mu)**2 + y**2 + z**2)
r2 = np.sqrt((x - 1 + mu)**2 + y**2 + z**2)
ax = 2*vy + x - (1 - mu)*(x + mu)/r1**3 - mu*(x - 1 + mu)/r2**3
ay = -2*vx + y - (1 - mu)*y/r1**3 - mu*y/r2**3
az = -(1 - mu)*z/r1**3 - mu*z/r2**3
return [vx, vy, vz, ax, ay, az]
# Initial state (approximate Halo orbit initial condition near L1)
state0 = [0.8234, 0, 0.1009, 0, 0.1923, 0]
mu = 0.01215 # Earth-Moon mass ratio
# Integrate for one orbital period
t_span = [0, 2.8]
sol = solve_ivp(
orbital_dynamics, t_span, state0,
args=(mu,), method='RK45', rtol=1e-9, atol=1e-12
)
print(f"Integration complete, {len(sol.t)} time points, endpoint deviation: {np.linalg.norm(sol.y[:3, -1] - state0[:3]):.6e}")
Parameter Optimization
from scipy.optimize import minimize, differential_evolution
def objective(params):
"""Objective function: minimize orbit deviation"""
x, y = params
return (x - 1)**2 + (y - 2)**2
# Local optimization (BFGS)
result = minimize(objective, [0, 0], method='BFGS')
print(f"Local optimum: {result.x}, objective: {result.fun:.6e}")
# Global optimization (differential evolution)
result = differential_evolution(objective, bounds=[(-5, 5), (-5, 5)])
print(f"Global optimum: {result.x}, objective: {result.fun:.6e}")
FFT Spectral Analysis (Orbit Period Detection)
import numpy as np
from scipy.fft import fft, fftfreq
# Simulate noisy orbital radius time series
t = np.linspace(0, 100, 1000)
r = 1.0 + 0.1 * np.sin(2 * np.pi * 0.05 * t) + 0.05 * np.random.randn(len(t))
# FFT analysis
yf = fft(r)
xf = fftfreq(len(t), t[1] - t[0])
# Find dominant frequency
positive_freqs = xf[:len(xf)//2]
amplitude = 2.0/len(t) * np.abs(yf[:len(yf)//2])
peak_idx = np.argmax(amplitude[1:]) + 1 # Exclude DC component
print(f"Dominant frequency: {positive_freqs[peak_idx]:.4f} Hz")
print(f"Corresponding period: {1/positive_freqs[peak_idx]:.2f} s")
Sparse Matrices (Space Object Proximity Analysis)
import numpy as np
from scipy.sparse import csr_matrix
from scipy.spatial.distance import cdist
# Simulate positions of 1000 space objects (km)
np.random.seed(42)
positions = np.random.randn(1000, 3) * 400000 # LEO-scale
# Compute pairwise distances, keep only close-approach pairs (< 100 km)
threshold = 100.0
dists = cdist(positions, positions)
row, col = np.where((dists < threshold) & (dists > 0))
data = dists[row, col]
# Build sparse distance matrix (1000×1000 with few non-zero elements)
sparse_dist = csr_matrix((data, (row, col)), shape=(1000, 1000))
density = sparse_dist.nnz / (1000 * 1000)
print(f"Space objects: {1000}")
print(f"Close-approach pairs: {sparse_dist.nnz}")
print(f"Matrix density: {density:.4%} (sparse saves {1 - density:.2%} memory)")
Statistical Analysis (Orbit Determination Errors)
import numpy as np
from scipy import stats
# Simulated position error samples from orbit determination (km)
errors = np.random.normal(loc=0.5, scale=2.0, size=1000)
# Descriptive statistics
print(f"Mean: {np.mean(errors):.3f} km")
print(f"Std dev: {np.std(errors):.3f} km")
print(f"95% confidence interval: {stats.norm.interval(0.95, loc=np.mean(errors), scale=stats.sem(errors))}")
# Normality test (Shapiro-Wilk limited to n ≤ 5000; use kstest for larger samples)
sample = np.random.choice(errors, size=500)
stat, p_value = stats.shapiro(sample)
print(f"Shapiro-Wilk test p-value: {p_value:.4f}")