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[bot]]

Codecov Report

All modified and coverable lines are covered by tests :white_check_mark:

Comparison is base (aea9d50) 81.56% compared to head (4441b8d) 81.61%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #193      +/-   ##
==========================================
+ Coverage   81.56%   81.61%   +0.05%     
==========================================
  Files          10       10              
  Lines         998     1001       +3     
==========================================
+ Hits          814      817       +3     
  Misses        184      184              

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[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).