Cleanup and refactor eps_fit_lorentzian.py

[oskooi]

What changed

nlopt →SciPy

• Replaced nlopt.opt(...) + the MMA/L-BFGS local-optimizer setup with scipy.optimize.least_squares(..., method="trf", bounds=(0, np.inf)) — the idiomatic tool for bound-constrained nonlinear least squares, and the natural fit for this residual-minimization problem.
• Dropped the nlopt-specific in-place-grad idiom. lorentzerr is replaced by a clean _residuals (real/imag-stacked residual vector) plus an analytic _jacobian passed to least_squares (verified against finite differences to ~1e-9).
• lorentzfit now returns the optimal params and the residual sum of squares Σ|Δ|² (= 2 * result.cost), preserving the original error metric.

Cleanup / refactor

• Split the monolithic __main__ block into focused functions: fit_medium (random-restart fitting → builds the Medium) and plot_fit (comparison plot), plus main()/parse_args().
• Added argparse so the hard-coded fitting parameters are now CLI flags: csv, --eps-inf, --num-lorentzians, --wl-min, --wl-max, --num-repeat, --seed, --output (defaults preserve the original behavior; --wl-min/--wl-max default to 400/700).
• Added a --seed so restarts are reproducible (the original used unseeded np.random).
• Kept the public lorentzfunc/lorentzfit names (nothing in the repo imports them, but they're the recognizable API), added type hints and docstrings, and preserved the Drude special-case (ω = 0).

Validation

• Analytic Jacobian correct (≈1e-9 vs. finite difference).
• Running the refactored script on the Sarkar TiO₂ data recovers the same fit (to be added to materials.py in a separate PR): γ→0 (lossless), and A/ω² gives σ = 2.6134 and 0.5625 — matching results from the unmodified script exactly, max |Δn| = 6.0×10⁻⁵
• Compiles, no nlopt references remain, -h works.

[oskooi]

Output of python3 eps_fit_lorentzian.py --help:

Using MPI version 3.1, 1 processes
usage: eps_fit_lorentzian.py [-h] [--eps-inf EPS_INF] [--num-lorentzians NUM_LORENTZIANS] [--wl-min WL_MIN]
                             [--wl-max WL_MAX] [--num-repeat NUM_REPEAT] [--seed SEED] [--output OUTPUT]
                             [csv]

Fit a complex refractive index profile to a sum of Lorentzian (and,
optionally, Drude) susceptibility terms and build a `Medium` from the result.

The complex permittivity is modeled as

    ε(f) = ε_inf + Σ_n  A_n / (ω_n² − f² − i f γ_n)

and the parameters (A_n, ω_n, γ_n) are found by nonlinear least-squares fitting
to tabulated data using SciPy (scipy.optimize.least_squares). The fit is
repeated from several random starting points to avoid poor local minima.

Usage:
    python eps_fit_lorentzian.py mymaterial.csv --num-lorentzians 2 \
        --eps-inf 1.0 --wl-min 400 --wl-max 700

The input CSV has three comma-separated columns and no header:
    wavelength (nm), real(n), imag(n)

positional arguments:
  csv                   input CSV: wavelength (nm), real(n), imag(n) (default: mymaterial.csv)

options:
  -h, --help            show this help message and exit
  --eps-inf EPS_INF     instantaneous (infinite-frequency) permittivity. Must be > 1.0 for stability; choose so
                        that min(real(eps - eps_inf)) is ~1.0 (default: 1.1)
  --num-lorentzians NUM_LORENTZIANS
                        number of Lorentzian terms in the fit (default: 2)
  --wl-min WL_MIN       minimum wavelength to fit, in nm (default: 400)
  --wl-max WL_MAX       maximum wavelength to fit, in nm (default: 700)
  --num-repeat NUM_REPEAT
                        number of random restarts of the local optimization (default: 30)
  --seed SEED           seed for the random restarts, for reproducitiblity (default: 0)
  --output OUTPUT       filename for the comparison plot (default: eps_fit_sample.png)

Elapsed run time = 0.0124 s