Rebranding vsini as resolution uncertainty in low resolution spectra

At (sufficiently) high spectral resolution, and with (sufficiently) rapidly rotating stars, the line broadening can be dominated by the (projected) rotational broadening, often parameterized by vsini. At low spectral resolution, or with low vsini sources, the instrumental resolution dominates the spectral broadening.

In the latter case, imperfections in spectral resolution estimates can add to model-data mismatches. Some of this residual structure can be taken up with the global Gaussian Process covariance structure. To first order, a finite vsini approximates additional broadening from under-estimated instrumental resolution. However, the right thing to do would be to rebrand vsini as a small, extra instrumental spectral resolution.

From a practical standpoint, this rebranding would require one small change in the update_theta() step-- modify the rotational broadening kernel to an instrumental resolution taper, mimicking the kernel used in the grid_tools.py load_flux() routine.

I have implemented a demonstration of this rebranding, and found one subtle bug. At some point in the code there is a self.ss[0] = 0.01 hack with the note # junk so we don't get a divide by zero error. My (admittedly limited, but seemingly reasonable) interpretation of this term is that it assigns a finite DC power to avoid numerical artifacts. In practice the relative value may affect your spectra, depending on the absolute mean level of your spectra. Indeed, in an experiment with-and-without this hardcode toggled, I see that the mean level of my model spectrum changes.

I'm not sure this affects many users-- I'm working with a customized and hacky version of Starfish (i.e. my fork's mix_model_omega2 branch... no shame in my bad branch hygiene though). So I wouldn't sound any alarms about this. For what it's worth, I do not get any divide by zero errors, so something seems to be working fine without the hack.

Finally, if this procedure were to be followed, it's conceivable that you could have either positive or negative perturbations to the kernel. Since we are using multiplications in Fourier space, in principle we could divide by a taper to get back high frequency structure, meaning that we would support either over- or under- estimates of the spectral resolution! I haven't really worked this out, but I think it would work. Cool, huh!

Start in update_theta() and compare its rotational broadening kernel with the instrumental-resolution taper used by grid_tools.py's load_flux() routine. Then investigate the self.ss[0] = 0.01 workaround and validate whether changing it alters model mean levels; done requires a defined kernel rebranding and a resolution of the divide-by-zero behavior.