Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove SpeX fitting #101

Open
mdecleir opened this issue Jan 12, 2022 · 0 comments
Open

Remove SpeX fitting #101

mdecleir opened this issue Jan 12, 2022 · 0 comments
Labels
bug Something isn't working enhancement New feature or request

Comments

@mdecleir
Copy link
Collaborator

This function to fit SpeX extinction data should be removed from this repository.
A more extended fitting is done in the spex_nir_extinction repo.

I expect that in the future, most stars will not have SpeX data available, so this should not be the default way to calculate A(V).

def fit_spex_ext(
self, amp_bounds=(-1.5, 1.5), index_bounds=(0.0, 5.0), AV_bounds=(0.0, 6.0)
):
"""
Fit the observed NIR extinction curve with a powerlaw model, based on the SpeX spectra
Parameters
----------
amp_bounds : tuple [default=(-1.5,1.5)]
Model amplitude bounds to be used in the fitting
index_bounds : tuple [default=(0.0,5.0)]
Powerlaw index bounds to be used in the fitting
AV_bounds : tuple [default=(0.0,6.0)]
A(V) bounds to be used in the fitting
Returns
-------
Updates self.model["waves", "exts", "residuals", "params"] and self.columns["AV"] with the fitting results:
- waves: np.ndarray with the wavelengths used in the fitting
- exts: np.ndarray with the fitted powerlaw model to the extinction curve
- residuals: np.ndarray with the fractional residuals, i.e. (data-fit)/fit
- params: tuple with the parameters (amplitude, alpha) if data in A(lambda)/A(V) or (amplitude, alpha, A(V)) if data in E(lambda-V)
"""
# retrieve the SpeX data, and sort the curve from short to long wavelengths
(waves, exts, exts_unc) = self.get_fitdata(["SpeX_SXD", "SpeX_LXD"])
indx = np.argsort(waves)
waves = waves[indx].value
exts = exts[indx]
exts_unc = exts_unc[indx]
# fit a powerlaw to the spectrum
if self.type == "alav":
func = PowerLaw1D(
fixed={"x_0": True},
bounds={"amplitude": amp_bounds, "alpha": index_bounds},
)
else:
func = (
PowerLaw1D(
fixed={"x_0": True},
bounds={"amplitude": amp_bounds, "alpha": index_bounds},
)
| AxAvToExv(bounds={"Av": AV_bounds})
)
fit = LevMarLSQFitter()
fit_result = fit(func, waves, exts, weights=1 / exts_unc)
# save the fitting results
self.model["waves"] = waves
self.model["exts"] = fit_result(waves)
self.model["residuals"] = exts - self.model["exts"]
if self.type == "alav":
self.model["params"] = (fit_result.amplitude.value, fit_result.alpha.value)
else: # in this case, fitted amplitude has to be multiplied by A(V) to get the "combined" amplitude
self.model["params"] = (
fit_result.amplitude_0.value * fit_result.Av_1.value,
fit_result.alpha_0.value,
fit_result.Av_1.value,
)
self.columns["AV"] = fit_result.Av_1.value

@karllark karllark added bug Something isn't working enhancement New feature or request labels Dec 31, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working enhancement New feature or request
Projects
None yet
Development

No branches or pull requests

2 participants