Rewrite crab spectrum as class

[adonath]

This PR introduces a CrabSpectrum class, with the different reference models represented as SpectrumModel. This has the advantage, that all the methods of SpectrumModel, such as .integral()
and energy_flux() can be reused. It has now full Quantity support as well.

[cdeil]

Thanks!

I've left two inline comments.

Concerning API: I probably would have made this a factory function, and just attached whatever other info you want to the SpectralModel objects after construction. But if you prefer the CrabSpectrum class where one accesses the model all the time, OK.

In any case: add one or two examples how to call it to the docstring?

[cdeil]

approximisation -> approximation ?

[adonath]

Fixed

[cdeil]

No need to indent lists, putting them at the same level as the surrounding text gives nice HTML output.

[adonath]

Fixed

[cdeil]

Does such a small eps yield numerically stable results?
I think I would have put something like 1e-3.

But I don't know, if you think or know the spectral index with that small eps gives stable output, OK to keep as-is.

[adonath]

I made it an option to the method. Default is now 1E-5

[cdeil]

After thinking about this some more, I'm thinking maybe keeping it a class it a good idea?
If we want to add extra methods like "convert given flux to Crab units" we could add it there and the method docstring would show up in the docs.

[adonath]

@cdeil The use case with the unit is what I had in mind as well. Right now the following works, which isn't too bad:

from astropy import units as u
from gammapy.spectrum import CrabSpectrum
from gammapy.spectrum.models import PowerLaw

# define power law model
amplitude = 1E-12 * u.Unit('1 / (cm2 s TeV)')
index = 2.3
reference = 1 * u.TeV
pwl = PowerLaw(index, amplitude, reference)

# define crab reference
crab = CrabSpectrum('hess_pl').model

# compute flux at 10 TeV in crab units
energy = 10 * u.TeV
flux_cu = (pwl(energy) / crab(energy)).to('%')
print(flux_cu)

# compute integral flux in crab units
emin, emax = [1, 10] * u.TeV
flux_int_cu = (pwl.integral(emin, emax) / crab.integral(emin, emax)).to('%')
print(flux_int_cu)

I think it's hard to implement a better API, because the flux unit is parametric and depends on what you do with the spectrum you want to measure the flux in crab units for...so something along the lines of:

pwl(energy).to('crab')
pwl.integral(emin, emax).to('crab')

Is very hard to implement and probably not worth the effort.

Possibly one could have something like:

pwl_cu = CrabSpectrum('meyer').use_as_unit_for(pwl)
pwl_cu(energy)

But it's not a super nice API either...

[cdeil]

I think what you have now, explaining how to do it via the docs, is a good solution.
Merge when travis-ci passes, right?