Skip to content

Commit

Permalink
Merge pull request #1079 from rzellem/issue_1005
Browse files Browse the repository at this point in the history 
Issue 1005 & 1073: Bypassing comp stars that fall of detector & adding Rp/R* to AAVSO priors
  • Loading branch information
@rzellem
rzellem authored Oct 28, 2022
2 parents a6df3d2 + 26da7cb commit 5562e5c
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 20 deletions.
44 changes: 27 additions & 17 deletions exotic/exotic.py
View file
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
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

0 comments on commit 5562e5c

Please sign in to comment.