Models with units.ipynb

For the Gaussian1D case, the following definition does not make sense (i.e., setting stddev to a unit that is incompatible with mean) but it is allowed, only to cause other exceptions to be raised on evaluation:

>>> g1 = Gaussian1D(amplitude=1*u.Jy, mean=5000*u.AA, stddev=0.1*u.Jy)
>>> g1(5000)
UnitsError: Can only apply 'subtract' function to ...
>>> g1(5000 * u.AA)
TypeError: Can only apply 'exp' function to dimensionless quantities

While I understand that compound model is future work, we also have to be careful such that massive refactoring won't be required in the future in order to make units work with them. In real life, most of us deal with compound models, so if we want to advertise that units are possible in models, we also have to make sure they will indeed be possible for compound models (albeit not implemented in this PR).

Currently:

>>> g1 = Gaussian1D(amplitude=1*u.Jy, mean=5000*u.AA, stddev=10*u.AA)
>>> g2 = Gaussian1D()
>>> gg = g1 + g2
>>> gg(5000)
UnitsError: Can only apply 'subtract' function to ...
>>> gg(5000 * u.AA)
TypeError: Can only apply 'exp' function to dimensionless quantities

or

>>> g1 = Gaussian1D(amplitude=1*u.Jy, mean=5000*u.AA, stddev=10*u.AA)
>>> g2 = Gaussian1D(amplitude=300*u.mJy, mean=30*u.MHz, stddev=10*u.Hz)
>>> gg = g1 + g2
>>> gg(5000 * u.AA, equivalencies={'x': u.spectral(), 'mean':u.spectral(), 'stddev':u.spectral()})
UnitConversionError: Can only apply 'subtract' function to ...