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

Jump to bottom

Issue 1005 & 1073: Bypassing comp stars that fall of detector & adding Rp/R* to AAVSO priors #1079

Merged
merged 4 commits into from
Oct 28, 2022
Merged

Issue 1005 & 1073: Bypassing comp stars that fall of detector & adding Rp/R* to AAVSO priors #1079

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9c11b84
Saving null values in array for comp stars that drift off detector
tamimfatahi Oct 25, 2022
6baef2f
Excluding comparison stars that fall of the detector during photometr…
tamimfatahi Oct 28, 2022
9c79dde
Merge branch 'develop' into issue_1005
tamimfatahi Oct 28, 2022
26da7cb
Fixed rprsUnc in json file
tamimfatahi Oct 28, 2022
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
44 changes: 27 additions & 17 deletions exotic/exotic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1928,7 +1928,10 @@ def main():

pix_coords = wcs_hdr.world_to_pixel_values(comp_radec[j][0], comp_radec[j][1])
cx, cy = pix_coords[0].take(0), pix_coords[1].take(0)
psf_data[ckey][i] = fit_centroid(imageData, [cx, cy])
if cx > 0 and cy > 0:
psf_data[ckey][i] = fit_centroid(imageData, [cx, cy])
else:
psf_data[ckey][i] = np.empty(7) * np.nan

if i != 0:
if not (tar_comp_dist[ckey][0] - 1 <= abs(int(psf_data[ckey][i][0]) - int(psf_data['target'][i][0])) <= tar_comp_dist[ckey][0] + 1 and
Expand All @@ -1951,7 +1954,10 @@ def main():
ckey = f"comp{j + 1}"

cx, cy = tform(coord)[0]
psf_data[ckey][i] = fit_centroid(imageData, [cx, cy])
if cx > 0 and cy > 0:
psf_data[ckey][i] = fit_centroid(imageData, [cx, cy])
else:
psf_data[ckey][i] = np.empty(7) * np.nan

if i == 0:
tar_comp_dist[ckey][0] = abs(int(psf_data[ckey][0][0]) - int(psf_data['target'][0][0]))
Expand All @@ -1973,10 +1979,13 @@ def main():
# loop through comp stars
for j in range(len(compStarList)):
ckey = f"comp{j + 1}"
aper_data[ckey][i][a][an], aper_data[f"{ckey}_bg"][i][a][an] = aperPhot(imageData,
psf_data[ckey][i, 0],
psf_data[ckey][i, 1],
aper, annulus)
if not np.isnan(psf_data[ckey][i][0]):
aper_data[ckey][i][a][an], \
aper_data[f"{ckey}_bg"][i][a][an] = aperPhot(imageData, psf_data[ckey][i, 0],
psf_data[ckey][i, 1], aper, annulus)
else:
aper_data[ckey][i][a][an] = np.nan
aper_data[f"{ckey}_bg"][i][a][an] = np.nan

# close file + delete from memory
hdul.close()
Expand Down Expand Up @@ -2099,9 +2108,10 @@ def main():
# try to fit data with comp star
for j in range(len(compStarList)):
ckey = f"comp{j + 1}"
cFlux = aper_data[ckey][:, a, an]
aper_mask = np.isfinite(aper_data[ckey][:, a, an])
cFlux = aper_data[ckey][aper_mask][:, a, an]

myfit, tFlux1, cFlux1 = fit_lightcurve(times, tFlux, cFlux, airmass, ld, pDict)
myfit, tFlux1, cFlux1 = fit_lightcurve(times[aper_mask], tFlux[aper_mask], cFlux, airmass[aper_mask], ld, pDict)

if j in vsp_num:
temp_ref_flux[j] = {
Expand Down Expand Up @@ -2134,17 +2144,17 @@ def main():
nonBJDTimes = np.copy(myfit.time)
# nonBJDPhases = np.copy(myfit.phase)
goodAirmasses = np.copy(myfit.airmass)
goodTargets = tFlux
goodTargets = tFlux[aper_mask]
goodReferences = cFlux
goodTUnc = tFlux ** 0.5
goodTUnc = tFlux[aper_mask] ** 0.5
goodRUnc = cFlux ** 0.5
# goodResids = myfit.residuals
bestlmfit = myfit

finXTargCent = psf_data["target"][:, 0]
finYTargCent = psf_data["target"][:, 1]
finXRefCent = psf_data[ckey][:, 0]
finYRefCent = psf_data[ckey][:, 1]
finXTargCent = psf_data["target"][aper_mask][:, 0]
finYTargCent = psf_data["target"][aper_mask][:, 1]
finXRefCent = psf_data[ckey][aper_mask][:, 0]
finYRefCent = psf_data[ckey][aper_mask][:, 1]

if ref_flux_opt or ref_flux_opt2:
if j in vsp_num:
Expand Down Expand Up @@ -2389,8 +2399,8 @@ def main():
log_info("\n*********************************************************")
log_info("FINAL PLANETARY PARAMETERS\n")
log_info(f" Mid-Transit Time [BJD_TDB]: {round_to_2(myfit.parameters['tmid'], myfit.errors['tmid'])} +/- {round_to_2(myfit.errors['tmid'])}")
log_info(f" Radius Ratio (Planet/Star) [Rp/Rs]: {round_to_2(myfit.parameters['rprs'], myfit.errors['rprs'])} +/- {round_to_2(myfit.errors['rprs'])}")
log_info(f" Transit depth [(Rp/Rs)^2]: {round_to_2(100. * (myfit.parameters['rprs'] ** 2.))} +/- {round_to_2(100. * 2. * myfit.parameters['rprs'] * myfit.errors['rprs'])} [%]")
log_info(f" Radius Ratio (Planet/Star) [Rp/R*]: {round_to_2(myfit.parameters['rprs'], myfit.errors['rprs'])} +/- {round_to_2(myfit.errors['rprs'])}")
log_info(f" Transit depth [(Rp/R*)^2]: {round_to_2(100. * (myfit.parameters['rprs'] ** 2.))} +/- {round_to_2(100. * 2. * myfit.parameters['rprs'] * myfit.errors['rprs'])} [%]")
log_info(f" Semi Major Axis/ Star Radius [a/Rs]: {round_to_2(myfit.parameters['ars'], myfit.errors['ars'])} +/- {round_to_2(myfit.errors['ars'])}")
log_info(f" Airmass coefficient 1: {round_to_2(myfit.parameters['a1'], myfit.errors['a1'])} +/- {round_to_2(myfit.errors['a1'])}")
log_info(f" Airmass coefficient 2: {round_to_2(myfit.parameters['a2'], myfit.errors['a2'])} +/- {round_to_2(myfit.errors['a2'])}")
Expand All @@ -2417,7 +2427,7 @@ def main():
f"Triangle_{pDict['pName']}_{exotic_infoDict['date']}.png")

if vsp_params:
VSPoutput_files = VSPOutputFiles(myfit, pDict, exotic_infoDict, durs, vsp_params)
VSPoutput_files = VSPOutputFiles(myfit, pDict, exotic_infoDict, vsp_params)
output_files = OutputFiles(myfit, pDict, exotic_infoDict, durs)
error_txt = "\n\tPlease report this issue on the Exoplanet Watch Slack Channel in #data-reductions."

Expand Down
10 changes: 7 additions & 3 deletions exotic/output_files.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ def final_planetary_params(self, phot_opt, vsp_params, comp_star=None, comp_coor
params_num = {
"Mid-Transit Time (Tmid)": f"{round_to_2(self.fit.parameters['tmid'], self.fit.errors['tmid'])} +/- "
f"{round_to_2(self.fit.errors['tmid'])} BJD_TDB",
"Ratio of Planet to Stellar Radius (Rp/Rs)": f"{round_to_2(self.fit.parameters['rprs'], self.fit.errors['rprs'])} +/- "
"Ratio of Planet to Stellar Radius (Rp/R*)": f"{round_to_2(self.fit.parameters['rprs'], self.fit.errors['rprs'])} +/- "
f"{round_to_2(self.fit.errors['rprs'])}",
"Transit depth (Rp/Rs)^2": f"{round_to_2(100. * (self.fit.parameters['rprs'] ** 2.))} +/- "
f"{round_to_2(100. * 2. * self.fit.parameters['rprs'] * self.fit.errors['rprs'])} [%]",
Expand Down Expand Up @@ -95,6 +95,7 @@ def aavso(self, comp_star, airmasses, ld0, ld1, ld2, ld3):
f"#FILTER={self.i_dict['filter']}\n"
f"#FILTER-XC={dumps(filter_dict)}\n"
f"#PRIORS=Period={round_to_2(self.p_dict['pPer'], self.p_dict['pPerUnc'])} +/- {round_to_2(self.p_dict['pPerUnc'])}"
f",Rp/R*={round_to_2(self.p_dict['rprs'], self.p_dict['rprsUnc'])} +/- {round_to_2(self.p_dict['rprsUnc'])}"
f",a/R*={round_to_2(self.p_dict['aRs'], self.p_dict['aRsUnc'])} +/- {round_to_2(self.p_dict['aRsUnc'])}"
f",inc={round_to_2(self.p_dict['inc'], self.p_dict['incUnc'])} +/- {round_to_2(self.p_dict['incUnc'])}"
f",ecc={round_to_2(self.p_dict['ecc'])}"
Expand Down Expand Up @@ -129,11 +130,10 @@ def aavso(self, comp_star, airmasses, ld0, ld1, ld2, ld3):


class VSPOutputFiles:
def __init__(self, fit, p_dict, i_dict, durs, vsp_params):
def __init__(self, fit, p_dict, i_dict, vsp_params):
self.fit = fit
self.p_dict = p_dict
self.i_dict = i_dict
self.durs = durs
self.dir = Path(self.i_dict['save'])
self.vsp_params = vsp_params

Expand Down Expand Up @@ -169,6 +169,10 @@ def aavso_dicts(planet_dict, fit, info_dict, durs, ld0, ld1, ld2, ld3):
'uncertainty': str(round_to_2(planet_dict['pPerUnc'])) if planet_dict['pPerUnc'] else planet_dict['pPerUnc'],
'units': "days"
},
'Rp/R*': {
'value': str(round_to_2(planet_dict['rprs'], planet_dict['rprsUnc'])),
'uncertainty': str(round_to_2(planet_dict['rprsUnc'])) if planet_dict['rprsUnc'] else planet_dict['rprsUnc'],
},
'a/R*': {
'value': str(round_to_2(planet_dict['aRs'], planet_dict['aRsUnc'])),
'uncertainty': str(round_to_2(planet_dict['aRsUnc'])) if planet_dict['aRsUnc'] else planet_dict['aRsUnc'],
Expand Down