Implement 2D sky subtraction using bivariate splines

[tepickering]

This PR implements the Kelson (2003) technique for subtracting the sky background from 2D spectra without the need for rectification or interpolation.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 81.61%. Comparing base (aea9d50) to head (4441b8d).
⚠️ Report is 432 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #193      +/-   ##
==========================================
- Coverage   81.91%   81.61%   -0.30%     
==========================================
  Files          10       10              
  Lines        1001     1001              
==========================================
- Hits          820      817       -3     
- Misses        181      184       +3     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[chris-simpson]

I really don't think this is an implementation of Kelson's method. You're fitting the BivariateSpline as a function of IMAGE pixel coordinates

spline_fit = LSQBivariateSpline(x.ravel(), y.ravel(), ...)

whereas Kelson says the spline should be fit as a function of x_r and y_r (page 694), coordinates which are orthogonal and aligned with the spatial and dispersion axes, using the known transformation from measuring the distortion from arc lamps or sky lines (page 691).