PAHFIT Covariance Kinematics

Other considerations:

Covariance matrix

scipy.optimize.leastsq often (almost always?) returns None for cov_x, which leads to no param_cov returned from astropy.modeling. As stated in the leastsq docs:

A value of None indicates a singular matrix, which means the curvature in parameters x is numerically flat.

It’s not entirely clear why such numerical flatness occurs. LevMarLSQFitter simply sets up and drives leastsq, forcibly truncating parameters at their boundaries during the fit, which may interfere with the numerically computed Jacobean/Hessian (this is also how MPFIT handled bounds). One idea is to use Hessian estimated by numdifftool after the model has converged at the best fit, and invert it for the covariance matrix. Would need to check normalization, as in this stackexchange post.

Kinematic line structure

At R~3000 (100km/s resolution), non-negligible kinematic broadening or sub-structure may be present in complex regions. This is clearly not the domain of the instrument pack, but squarely in the science pack regime. But it woudl also add too much complexity to PAHFIT. A simple handling would be to allow (at model construction time) some sigma_velocity, for quadrature addition to the instrument-pack-returned widths, to slightly dilate all lines/dust features. Or this could be set per line or group. Note that this doesn’t belong in the science pack, since it varies source to source.

Much (much) fancier methods exist. Many other tools specialize, and discriminate between various species in their kinematics. More sensible to zoom in on individual (sets of) lines and use one of those bespoke kinematic decomposition tools.