AbstractLambertSolver <: SciMLBase.AbstractSciMLAlgorithm
├── GoodingSolver      # Gooding 1990 - robust, accurate
├── IzzoSolver         # Izzo 2015 - very fast
├── ValladoSolver      # Vallado 2013 - guaranteed convergence (bisection)
├── AroraSolver        # Arora 2013 - fast cosine transformation
├── BattinSolver       # Battin 1984 - handles 180° singularity
├── GaussSolver        # Gauss 1809 - historical, limited accuracy
└── AvanziniSolver     # Avanzini 2008 - eccentricity-based, single rev

AbstractAstroProblem <: SciMLBase.AbstractSciMLProblem
└── LambertProblem     # (μ, r1, r2, tof)

LambertSolution        # (v1, v2, numiter, retcode)
LambertIterator        # Stateful iterator for solve!

src/
  Lambert.jl            # Module, imports SciMLBase.solve/remake/init/solve!
  lambert_problem.jl    # LambertProblem, LambertSolution, LambertIterator
  utils.jl              # Stumpff functions, geometry, normalization, Lagrange coefficients
  porkchop.jl           # PorkchopGrid, porkchop_grid() with EnsembleProblem
  {solver_name}_solver.jl  # One file per solver algorithm
ext/
  LambertPlotsExt.jl    # Weak dep for porkchop plot visualization

# Object-oriented (recommended)
prob = LambertProblem(μ, r1, r2, tof)
sol = solve(prob, GoodingSolver())
# sol.v1, sol.v2, sol.numiter, sol.retcode

# Remake for parameter sweeps
prob2 = remake(prob; tof=new_tof)

v1, v2, numiter, converged = gooding1990(μ, r1, r2, tof; M=0, prograde=true)

@with_kw struct GoodingSolver <: AbstractLambertSolver
    M::Int = 0             # Number of complete revolutions
    prograde::Bool = true  # Prograde transfer
    low_path::Bool = true  # Low-energy path (multi-rev)
    maxiter::Int = 35
    atol::Float64 = 1e-5
    rtol::Float64 = 1e-7
end

grid = porkchop_grid(μ, r1_func, r2_func, dep_dates, arr_dates, solver;
    ensemble_method=EnsembleThreads())