def gaussian_model(xaxis, amplitude, offset, width):
    amplitude = u.Quantity(amplitude, u.K)
    offset = u.Quantity(offset, u.km/u.s)
    width = u.Quantity(width, u.km/u.s)
   
    return amplitude*np.exp(-(xaxis-offset)**2/(2.*width**2))

from specutils.io import fits
spec = fits.read_fits_spectrum1d('gbt_1d.fits')

from astropy import units as u

%%bash
which conda

import specutils
import numpy
import astropy
specutils.__version__, astropy.__version__, numpy.__version__, astropy.__path__

spec.velocity

model = gaussian_model(spec.velocity, amplitude=5*u.K, offset=20*u.km/u.s, width=5*u.km/u.s)

%matplotlib inline
import pylab as pl

pl.plot(spec.velocity, spec.flux, 'k-')
pl.plot(spec.velocity, model)

def cost_function(params):
    return ((spec.flux*u.K-gaussian_model(spec.velocity, *params))**2).sum().value
    

from scipy.optimize import curve_fit, minimize

result = minimize(cost_function, (-5, 20, 20))

result

best_fit_parameters = result.x

best_fit_model = gaussian_model(spec.velocity, *best_fit_parameters)
best_fit_model

pl.plot(spec.velocity, spec.flux, 'k-')
pl.plot(spec.velocity, best_fit_model)
pl.xlim(-30,50)

"""
Explanation: Fitting data to models

Build a model
Create a "fitness function", i.e. something that returns a scalar "distance" between the model and the data
Apply an "optimizer" to get the best-fit parameters
End of explanation
"""


from scipy.optimize import curve_fit

# curve_fit does not play well with units
#result_curve_fit = curve_fit(gaussian_model, spec.velocity, spec.flux*u.K, p0=(-5, 20, 20))

"""
Explanation: Fitting Tools

scipy.optimize.curve_fit: simpler fitter, lets you skip the cost-fitting section
astropy.models
pyspeckit.specfit

scipy.optimize.curve_fit
End of explanation
"""


import pyspeckit
sp = pyspeckit.Spectrum('gbt_1d.fits')

sp.plotter(xmin=-40, xmax=50)
sp.specfit(guesses=(-5, 20, 20))

"""
Explanation: astropy.modeling
pyspeckit.specfit
End of explanation
"""


