Proposed spectrum arithmetic tools

(I'll need to fix some of the commit messages)

This PR includes a proposed spectrum arithmetic toolset.  It will allow you to do something like:

```
# assume Spectrum1D() magically creates a valid spectrum...
sp1 = Spectrum1D()
sp2 = Spectrum1D()
# this works iff all(sp1.disp == sp2.disp)
sp3 = sp1+sp2
# this works if sp1.disp ~ sp2.disp
sp1._arithmetic_threshold(1,'angstrom')
sp3 = sp1-sp2
```

There probably needs to be a `force` option as well, to override the dispaxis agreement entirely.

This PR includes at least one assumption that is not necessarily satisfied: that `disp` has a `units` attribute (it should, but perhaps not by that name)

In the comments, let's discuss first whether this is a good idea at all (I think so), then implementation if we give it a go-ahead.

Also, tell me if `nddata` will do all this implicitly... I am still unfamiliar with its capabilities.