Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save saimn/bc2190e6e5b7a03284ae1d2ac0eec37b to your computer and use it in GitHub Desktop.
Save saimn/bc2190e6e5b7a03284ae1d2ac0eec37b to your computer and use it in GitHub Desktop.
custom model to fit a 2D moffat model to data, equate that to F_OIII_xy, then pass it to a 1D custom Gaussian equation to get a model spectrum. The idea is to have x, y and l (wavelength) as inputs, then the amplitude is fixed as it is calculated. When I run it I normally get mismatched dimension errors: wavelength is a 271 long list, and then x…
def Moffat_3d_test(x, mean, stddev, Gauss_bkg, Gauss_grad,
Moffat_amplitude, x_0, y_0, gamma, alpha, Moffat_bkg):
# Moffat
rr_gg = ((x_fit - x_0)**2 + (y_fit - y_0)**2) / gamma**2
F_OIII_xy = Moffat_amplitude * (1 + rr_gg)**(-alpha) + Moffat_bkg
# Prep for Gauss 1D
Gauss_std = np.sqrt(stddev**2 + std_MUSE**2)
A_OIII_xy = F_OIII_xy / (np.sqrt(2*np.pi) * Gauss_std)
check_1.append(A_OIII_xy)
model_spectra = []
for Amp in A_OIII_xy:
model_spectrum = (Gauss_bkg + (Gauss_grad * x) + np.abs(Amp) * np.exp(- 0.5 * (x - mean)** 2 / Gauss_std**2.) +
(np.abs(Amp)/3) * np.exp(- 0.5 * (x - (mean - 47.9399))** 2 / Gauss_std**2.))
model_spectra.append(model_spectrum)
return model_spectra
#%%
new_model = Model(Moffat_3d_test, independant_vars=["x"])
test_data = np.array(PNe_spectra_list[0][202])
pars = new_model.make_params(mean=5007., stddev=0.9, Gauss_bkg=1., Gauss_grad=0.1,
Moffat_amplitude=20., x_0=8., y_0=8., gamma=4.47, alpha=2.39, Moffat_bkg=10)
pars["alpha"].vary = False
pars["gamma"].vary = False
pars["stddev"].vary = False
#pars["stddev"].min = 0.15
#pars["stddev"].max = 1.0
result = new_model.fit(test_data, x=wavelength, params=pars)
best_fit = result.best_fit
test_1 = best_fit[0]
plt.plot(wavelength, test_data, "k")
plt.plot(wavelength, test_1, "r-")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment