Thank you Jack, for sharing your code.
While GILDAS/CLASS do a god job in working with single spectra, it is too cumbersome to use on large number of spectra. At least I have not yet found a good method to process a large number of spectra. I should probably dig deeper into scripting within CLASS.
Due to this, I have written a Python script which is somewhat similar to yours. In my case I am using "lmfit" to do the baseline fitting. There is one challenge with fitting the baseline which can also be seen in your example: One needs to exclude the range of the actual signal from the fit, otherwise one will get a "bump" in the area of the actual signal. This can be seen in your 3m_Lime-H1_spec_10994_1.png plot. I avoid this by setting a window around the signal where the weighting of data points is set to zero. This is pretty much the same as one does when using CLASS.
The problem with that, however, is that the window needs to be placed differently depending on the doppler shift of the line. I still need to do some work to see how this could be done automatically in a reliable fashion.
My code is not in a state yet where I can share it. There are too many things which are specific to our data and too many places where the code needs to be adapted depending on circumstances.
So there is still some work ahead.....
Wolfgang