Search code examples
pythongaussianastronomyspectrum

Fitting data to a gaussian profile


I have been trying to fit a gaussian into my spectrum. (intensity v/s velocity spectrum)

spectrum

New fitted spectrum

I used the following code to fit the data to a gaussian profile. However, as seen in the result only one data point was included in the fit. Is there anything I can do so that i can include more points into the gaussian.

from numpy import exp, linspace, random
import matplotlib.pyplot as plt
def gaussian(x, amp, cen, wid):
    return amp * exp(-(x-cen)**2 / wid)

from scipy.optimize import curve_fit
x = velocity
y = data
print(x)
print(y)

init_vals = [0.00950554, 60000, 35]  # for [amp, cen, wid]
best_vals, covar = curve_fit(gaussian, x, y, p0=init_vals)
print(best_vals)
print(repr(covar))

ym = gaussian(x, best_vals[0], best_vals[1], best_vals[2])
fig = plt.figure()
ax = fig.add_subplot(211)
ax.plot(x,y)
ax = fig.add_subplot(212)
ax.plot(x, y, c='k', label='Function')
ax.plot(x, ym, c='r', label='Best fit')
plt.legend()
plt.show()

Data points:

x: [-5.99999993e+04 -4.99999993e+04 -3.99999993e+04 -2.99999993e+04
 -1.99999993e+04 -9.99999934e+03  6.65010004e-04  1.00000007e+04
  2.00000007e+04  3.00000007e+04  4.00000007e+04  5.00000007e+04
  6.00000007e+04  7.00000007e+04  8.00000007e+04  9.00000007e+04
  1.00000001e+05  1.10000001e+05  1.20000001e+05  1.30000001e+05
  1.40000001e+05]

y: [ 0.00056511 -0.00098584 -0.00325616 -0.00101042  0.00168894 -0.00097406
 -0.00134408  0.00128847 -0.00111633 -0.00151621  0.00299326  0.00916455
  0.00960554  0.00317363  0.00311124 -0.00080881  0.00215932  0.00596419
 -0.00192256 -0.00190138 -0.00013216]

These are the data points of the spectrum. Can anybody help me get a better gaussian fit into the data. I have been trying everything I can.

Thanks a lot.


Solution

  • The first step is always to plot the data, which you already did. The next is to guess initial values. Those for amp and cen look reasonably, if compared to the plot. But what about wid? It is 2 times the width of the distribution SQUARED. From the plot, the width itself must be of odrer of several thousand. If squared, it may reach 10^7, times 2 gives 2*10^7 as the initial value. Very far from your 35!

    One possible solution:

    amp = 0.0106283
    cen = 55784
    wid = 1.92911e+08

    Plot:

    Plot of itting function