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.

Sign up for GitHub

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

Save results #75

Merged
merged 4 commits into from
Mar 5, 2021
Merged

Save results #75

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
faf2eca
Save AV, RV and parameter uncertainties
mdecleir Feb 24, 2021
d703fae
:small update
mdecleir Feb 24, 2021
7dd9214
Added uncertainties to read function
mdecleir Feb 25, 2021
f978d76
Saving AV uncertainties, and give preference to columns given as argu…
mdecleir Mar 4, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 52 additions & 12 deletions measure_extinction/extdata.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -475,7 +475,7 @@ def calc_RV(self):
# obtain or calculate A(V)
if "AV" not in self.columns.keys():
self.calc_AV()
av = self.columns["AV"]
av = _get_column_val(self.columns["AV"])

# obtain E(B-V)
dwaves = np.absolute(self.waves["BAND"] - 0.438 * u.micron)
Expand Down Expand Up @@ -657,6 +657,7 @@ def save(
"rv": ("RV", "R(V)"),
"rvunc": ("RV_UNC", "R(V) uncertainty"),
}
# give preference to the column info that is given as argument to the save function
if column_info is not None:
for ckey in column_info.keys():
if ckey in ext_col_info.keys():
Expand All @@ -666,6 +667,29 @@ def save(
hval.append(column_info[ckey])
else:
print(ckey + " not supported for saving extcurves")
else: # save the column info if available in the extdata object
if "AV" in self.columns.keys():
hname.append("AV")
hcomment.append("V-band extinction A(V)")
if isinstance(self.columns["AV"], tuple):
hval.append(self.columns["AV"][0])
if len(self.columns["AV"]) == 2:
hname.append("AV_UNC")
hcomment.append("A(V) uncertainty")
hval.append(self.columns["AV"][1])
elif len(self.columns["AV"]) == 3:
hname.append("AV_L")
hcomment.append("A(V) lower uncertainty")
hval.append(self.columns["AV"][1])
hname.append("AV_U")
hcomment.append("A(V) upper uncertainty")
hval.append(self.columns["AV"][2])
else:
hval.append(self.columns["AV"])
if "RV" in self.columns.keys():
hname.append("RV")
hcomment.append("total-to-selective extintion R(V)")
hval.append(self.columns["RV"])

# legacy save param keywords
if fm90_best_params is not None:
Expand Down Expand Up @@ -765,6 +789,7 @@ def save(
)
columns = fits.ColDefs([col1, col2, col3])
tbhdu = fits.BinTableHDU.from_columns(columns)
# add the paramaters and their uncertainties
for param in self.model["params"]:
tbhdu.header.set(
param.name[:8],
Expand All @@ -775,7 +800,20 @@ def save(
+ " | fixed="
+ str(param.fixed),
)
tbhdu.header.set(
param.name[:6] + "_L",
param.unc_minus,
param.name + " lower uncertainty",
)
tbhdu.header.set(
param.name[:6] + "_U",
param.unc_plus,
param.name + " upper uncertainty",
)
tbhdu.header.set("MOD_TYPE", self.model["type"], "Type of fitted model")
tbhdu.header.set(
"chi2", self.model["chi2"], "Chi squared for the fitted model"
)
tbhdu.header.set("EXTNAME", "MODEXT", "Fitted model extinction")
hdulist.append(tbhdu)

Expand Down Expand Up @@ -841,9 +879,11 @@ def read(self, ext_filename):

# get the fitted model if available
if "MODEXT" in extnames:
self.model["waves"] = hdulist["MODEXT"].data["MOD_WAVE"]
self.model["exts"] = hdulist["MODEXT"].data["MOD_EXT"]
self.model["residuals"] = hdulist["MODEXT"].data["RESIDUAL"]
data = hdulist["MODEXT"].data
hdr = hdulist["MODEXT"].header
self.model["waves"] = data["MOD_WAVE"]
self.model["exts"] = data["MOD_EXT"]
self.model["residuals"] = data["RESIDUAL"]
self.model["params"] = []
paramkeys = [
"AMPLITUD",
Expand All @@ -855,19 +895,19 @@ def read(self, ext_filename):
"ASYM",
"AV",
]
hdr = hdulist["MODEXT"].header
self.model["type"] = hdr["MOD_TYPE"]
for paramkey in paramkeys:
if paramkey in list(hdr.keys()):
comment = hdr.comments[paramkey].split(" |")
self.model["params"].append(
Parameter(
name=comment[0],
default=hdr[paramkey],
bounds=comment[1].split("=")[1],
fixed=comment[2].split("=")[1],
)
param = Parameter(
name=comment[0],
default=hdr[paramkey],
bounds=comment[1].split("=")[1],
fixed=comment[2].split("=")[1],
)
param.unc_minus = hdr[paramkey[:6] + "_L"]
param.unc_plus = hdr[paramkey[:6] + "_U"]
self.model["params"].append(param)

# get the columns p50 +unc -unc fit parameters if they exist
if pheader.get("AV_p50"):
Expand Down