I don't understand why the operator is so not caring about the type of other that is allowed for proper math. self.add() here is all the way in astropy.nddata and it has no knowledge of what a spectral axis is, so any spectral axis check has to happen before we call self.add(). But it is not clear to me why you accept seemingly any random kind of data. How am I supposed to know where the spectral axis is in any given kind of things accepted here? The arithmetric tests only testing Spectrum1D inputs.

Also, this is out of scope, right?

My proposal for __add__ and __sub__:

• Throw error if other is not Spectrum1D.
• As Ricky suggested, strict check for spectral axis values (match must be 1-to-1 within tolerance). Is sys.float_info.epsilon good here?
• For unit compatibility, maybe we can delegate to nddata but I'll have to double check.
• Write better tests also to cover edge cases.

Is __mul__ and __div__ in scope here? The original post above only mentioned - operator. Those will have very different logic from __add__ and __sub__.

• Other must be unitless or equivalent (e.g., maybe sr is okay if it doesn't mess with the rest of the specutils calculations). Unitless scalar should be acceptable. Otherwise, it must be Spectrum1D and spectral axis must match as above.
• For division only, if other is not unitless, then it must have the same unit as self. The result here would give you a unitless Spectrum1D. Do you support unitless Spectrum1D at all? I am not talking about count. It has to be u.dimensionless_unscaled.
• Write better tests also to cover edge cases.

Do I have to worry about magnitude flux units in specutils?

Thanks for writing up a proposal @pllim. Here are my thoughts:

My proposal for __add__ and __sub__:

* Throw error if `other` is not `Spectrum1D`.

I think we should also allow other to be a Quantity of matching unit/shape to the flux array. I'm thinking about the case where someone wants to do continuum-subtraction - I don't think the user should necessarily have to put their continuum into a Spectrum1D after average a flux array or something. If the other operand is also a Spectrum1D, then we should check that the axes match.

* As Ricky suggested, strict check for spectral axis values (match must be 1-to-1 within tolerance). Is `sys.float_info.epsilon` good here?

I haven't used this before, from a quick look it's probably sufficient.

Is __mul__ and __div__ in scope here? The original post above only mentioned - operator. Those will have very different logic from __add__ and __sub__.

These are in scope, but the unit compatibility checks you're proposing below aren't. If you want to add those now, go for it, but the scope of this issue is checking that the spectral axes match.

* Other must be unitless or equivalent (e.g., maybe `sr` is okay if it doesn't mess with the rest of the `specutils` calculations). Unitless scalar should be acceptable. Otherwise, it must be `Spectrum1D` and spectral axis must match as above.

* For division only, if other is not unitless, then it must have the same unit as `self`. The result here would give you a unitless `Spectrum1D`. Do you support unitless `Spectrum1D` at all? I am not talking about `count`. It has to be `u.dimensionless_unscaled`.

We do use u.dimensionless_unscaled in a few places.

Do I have to worry about magnitude flux units in specutils?

No.

I think we should also allow other to be a Quantity of matching unit/shape

Can you please clarify? It makes sense to allow scalar Quantity, but if the Quantity has an array, it feels wrong to just assume user knows the arrays maps to the same spectral axis.

the unit compatibility checks you're proposing below aren't

Well, I guess they can also be addressed in a different PR...

I think we should also allow other to be a Quantity of matching unit/shape

Can you please clarify? It makes sense to allow scalar Quantity, but if the Quantity has an array, it feels wrong to just assume user knows the arrays maps to the same spectral axis.

Here's my theoretical workflow: I have a 2D or 3D Spectrum1D, I take an average of some background region, I create a numpy array of the same shape as flux and assign that average to every array element and turn it into a Quantity using the flux unit. Now I want to subtract it from the original spectrum - I know it's something that I want to apply to every element of the spectrum and might be annoyed that I have to take the extra step of getting the spectral axis from the original to turn this into a Spectrum1D just to subtract it.

It seems unlikely (thought not impossible) that I would just coincidentally have something of the exact same shape as my original spectrum and be bamboozled by addition/subtraction working when it shouldn't.

the unit compatibility checks you're proposing below aren't

Well, I guess they can also be addressed in a different PR...

Right, I'm not saying it shouldn't be done, I'm just saying if you want to minimize the scope of your PR for this issue, the units stuff can be separate work.

To enforce good hygiene, i.e. to avoid adding apples and oranges, I think we should be restrictive on the quantities that can be added to/subtracted from Spectrum 1D's:

1. They should have matching flux units
2. They should have matching spectral axes and units.

For multiplication and division:

1. flux units may be present, or not.
2. They should have matching spectral axes and units.

While it may be an added burden to attach a spectral axis for the cases discussed by @rosteen above, it only takes one extra step.

@eteq @nmearl @keflavich Any thoughts you want to throw in the pile here?

Yeah, I agree with @pllim's suggestion that Spectrum1Ds should be required for non-scalars. This definitely makes it harder for users, so we need clear, explicit documentation about how to turn an array, or a model, into an appropriate Spectrum1D.

Scalar addition & subtraction should be allowed with matching units. Scalar multiplication and division should be allowed with or without units.

Scalar multiplication and division should be allowed with or without units

What does it even mean physically when you multiply Jy by any random unit, like another Jy? Or Jy with Kelvin? Multiplication/division only makes sense to me when the other operand is unitless. Am I missing something?

What does it even mean physically when you multiply Jy by any random unit, like another Jy? Or Jy with Kelvin? Multiplication/division only makes sense to me when the other operand is unitless. Am I missing something?

You're right, that's a pretty rare case. There are some where I do that, though. For example, if you want to integrate over a spectrum, you might do so by multiplying by the spectral width of the pixel and then summing. If the pixels have varying width, you need to be able to multiply by a vector of pixel widths. So multiplying by, e.g., km/s or nm or Hz is sometimes reasonable.

There may be other cases - I don't think we should strictly exclude multiplication with unit'd quantities.

Ok, sounds like I'm outvoted on requiring Spectrum1D for non-scalars 👍

I think I have implemented all the requests here. Please review #998 . The tests in #998 should give you a good idea of the expected new behaviors. Thanks!