Skip to content

Commit

Permalink
Merge pull request #475 from kevin218/dev/kbs
Browse files Browse the repository at this point in the history 
Write S5 LC files to a single Xarray file in Stage 6
  • Loading branch information
@kevin218
kevin218 authored Oct 21, 2022
2 parents cda5df5 + b22fbf0 commit 3389d9b
Show file tree
Hide file tree
Showing 2 changed files with 100 additions and 23 deletions.
7 changes: 5 additions & 2 deletions src/eureka/S5_lightcurve_fitting/s5_fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ def fitlc(eventlabel, ecf_path=None, s4_meta=None, input_meta=None):
Enabled Astraeus
'''
print("\nStarting Stage 5: Light Curve Fitting\n")

if input_meta is None:
# Load Eureka! control file and store values in Event object
ecffile = 'S5_' + eventlabel + '.ecf'
Expand Down Expand Up @@ -158,6 +158,9 @@ def fitlc(eventlabel, ecf_path=None, s4_meta=None, input_meta=None):
else:
time_units = lc.data.attrs['time_units']
meta.time = lc.time.values
# Record units for Stage 6
meta.time_units = time_units
meta.wave_units = lc.data.attrs['wave_units']

# If any of the parameters' ptypes are set to 'white_free', enforce
# a Gaussian prior based on a white-light light curve fit. If any
Expand Down Expand Up @@ -230,7 +233,7 @@ def fitlc(eventlabel, ecf_path=None, s4_meta=None, input_meta=None):
try:
ld_coeffs = np.loadtxt(ld_fix_file)
except FileNotFoundError:
raise Exception("The limb-darkening file " + ld_fix_file +
raise Exception("The limb-darkening file " + ld_fix_file +
" could not be found.")
else:
ld_coeffs = None
Expand Down
116 changes: 95 additions & 21 deletions src/eureka/S6_planet_spectra/s6_spectra.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -287,7 +287,7 @@ def plot_spectra(eventlabel, ecf_path=None, s5_meta=None, input_meta=None):
meta.y_label_unit = ' '+meta.y_label_unit

if (rcParams['text.usetex'] and
(meta.y_label_unit.count(r'\%') !=
(meta.y_label_unit.count(r'\%') !=
meta.y_label_unit.count('%'))):
# Need to escape % with \ for LaTeX
meta.y_label_unit = meta.y_label_unit.replace('%', r'\%')
Expand Down Expand Up @@ -321,6 +321,7 @@ def plot_spectra(eventlabel, ecf_path=None, s5_meta=None, input_meta=None):
scale_height, meta.planet_R0)

save_table(meta, log)
convert_s5_LC(meta, log)

# Save results
log.writelog('Saving results')
Expand Down Expand Up @@ -465,7 +466,7 @@ def parse_unshared_saves(meta, log, fit_methods):
return meta
spectrum_median.extend(median)
spectrum_err.extend(err.T)

meta.spectrum_median = np.array(spectrum_median)
meta.spectrum_err = np.array(spectrum_err).T

Expand Down Expand Up @@ -523,6 +524,79 @@ def compute_timescale(meta):
return meta


def convert_s5_LC(meta, log):
'''
Loads spectroscopic light curves save files from S5 and write as
single Xarray save file.
'''
event_ap_bg = (meta.eventlabel+"_ap"+str(meta.spec_hw_val)+'_bg' +
str(meta.bg_hw_val))

if meta.sharedp:
niter = 1
else:
niter = meta.nspecchan
wavelengths = np.zeros(niter)
bin_widths = np.zeros(niter)
for ch in range(niter):
# Get the channel key
if meta.sharedp:
channel_key = 'shared'
else:
nzfill = int(np.floor(np.log10(meta.nspecchan))+1)
channel_key = 'ch'+str(ch).zfill(nzfill)

# Load text file
fname = f'S5_{event_ap_bg}_Table_Save_{channel_key}.txt'
full_fname = meta.inputdir+fname
lc_table = astropytable.readtable(full_fname)

# Assign known values to array
wavelengths[ch] = lc_table['wavelength'][0]
bin_widths[ch] = lc_table['bin_width'][0]
lc_table.remove_column('wavelength')
lc_table.remove_column('bin_width')
if ch == 0:
# Record time array
time = lc_table['time']
lc_table.remove_column('time')
# Get remaining column names and number
colnames = lc_table.colnames
n_col = len(colnames)
n_int = len(time)
# Create numpy array to hold data
lc_array = np.zeros((n_col, niter, n_int))
else:
lc_table.remove_column('time')
# Assign remaining values to array
for i, col in enumerate(lc_table.itercols()):
lc_array[i, ch] = col.value

# Create Xarray DataArrays and dictionary
flux_units = 'Normalized'
if hasattr(meta, 'wave_units'):
wave_units = meta.wave_units
else:
wave_units = 'microns'
if hasattr(meta, 'time_units'):
time_units = meta.time_units
else:
time_units = 'BMJD'
lc_da = []
dict = {}
for i in range(n_col):
lc_da.append(xrio.makeLCDA(lc_array[i], wavelengths, time, flux_units,
wave_units, time_units, name=colnames[i]))
dict[colnames[i]] = lc_da[-1]

# Create Xarray Dataset
ds = xrio.makeDataset(dict)
# Write to file
meta.lc_filename_s6 = (meta.outputdir+'S6_'+event_ap_bg + "_LC")
xrio.writeXR(meta.lc_filename_s6, ds)
return meta


def load_s5_saves(meta, log, fit_methods):
if 'dynesty' in fit_methods:
fitter = 'dynesty'
Expand Down Expand Up @@ -550,7 +624,7 @@ def load_s5_saves(meta, log, fit_methods):
channel_key = 'ch'+str(ch).zfill(nzfill)

fname = f'S5_{fitter}_samples_{channel_key}'

# Load HDF5 files
full_fname = meta.inputdir+fname+'.h5'
ds = xrio.readXR(full_fname, verbose=False)
Expand All @@ -560,7 +634,7 @@ def load_s5_saves(meta, log, fit_methods):
sample = hf['samples'][:]
# Need to figure out which columns are which
fname = f'S5_{fitter}_fitparams_{channel_key}.csv'
fitted_values = pd.read_csv(meta.inputdir+fname,
fitted_values = pd.read_csv(meta.inputdir+fname,
escapechar='#',
skipinitialspace=True)
full_keys = np.array(fitted_values["Parameter"])
Expand All @@ -581,7 +655,7 @@ def load_s5_saves(meta, log, fit_methods):
samples = np.array(meta.spectrum_median)
if all(x is None for x in samples):
samples = np.zeros((meta.nspecchan, 0))

return np.array(samples)


Expand All @@ -594,7 +668,7 @@ def compute_offset(meta, log, fit_methods, nsamp=1e4):
suffix = '2'
else:
suffix = '1'

# Load sine amplitude
meta.y_param = 'AmpSin'+suffix
ampsin = load_s5_saves(meta, log, fit_methods)
Expand All @@ -604,7 +678,7 @@ def compute_offset(meta, log, fit_methods, nsamp=1e4):
'fitted parameters')
log.writelog(f' Skipping {y_param}')
return meta

# Load cosine amplitude
meta.y_param = 'AmpCos'+suffix
ampcos = load_s5_saves(meta, log, fit_methods)
Expand All @@ -614,7 +688,7 @@ def compute_offset(meta, log, fit_methods, nsamp=1e4):
'fitted parameters')
log.writelog(f' Skipping {y_param}')
return meta

# Reset meta.y_param
meta.y_param = y_param

Expand All @@ -630,14 +704,14 @@ def compute_offset(meta, log, fit_methods, nsamp=1e4):
offset[2] = offset[0]-offset[2]
meta.spectrum_median.append(offset[0])
meta.spectrum_err.append(offset[1:])

# Convert the lists to an array
meta.spectrum_median = np.array(meta.spectrum_median)
if meta.fitter == 'lsq':
meta.spectrum_err = np.ones((2, meta.nspecchan))*np.nan
else:
meta.spectrum_err = np.array(meta.spectrum_err).T

return meta


Expand Down Expand Up @@ -681,28 +755,28 @@ def compute_amp(meta, log, fit_methods, nsamp=1e4):
'fitted parameters')
log.writelog(f' Skipping {y_param}')
return meta

# Reset meta.y_param
meta.y_param = y_param

meta.spectrum_median = []
meta.spectrum_err = []

for i in range(meta.nspecchan):
amps = fp[i]*np.sqrt(ampcos[i]**2+ampsin[i]**2)*2
amp = np.percentile(np.array(amps), [16, 50, 84])[[1, 2, 0]]
amp[1] -= amp[0]
amp[2] = amp[0]-amp[2]
meta.spectrum_median.append(amp[0])
meta.spectrum_err.append(amp[1:])

# Convert the lists to an array
meta.spectrum_median = np.array(meta.spectrum_median)
if meta.fitter == 'lsq':
meta.spectrum_err = np.ones((2, meta.nspecchan))*np.nan
else:
meta.spectrum_err = np.array(meta.spectrum_err).T

return meta


Expand All @@ -719,7 +793,7 @@ def compute_fn(meta, log, fit_methods, nsamp=1e4):
'fitted parameters')
log.writelog(f' Skipping {y_param}')
return meta

# Load cosine amplitude
meta.y_param = 'AmpCos1'
ampcos = load_s5_saves(meta, log, fit_methods)
Expand All @@ -729,7 +803,7 @@ def compute_fn(meta, log, fit_methods, nsamp=1e4):
'fitted parameters')
log.writelog(f' Skipping {y_param}')
return meta

# Reset meta.y_param
meta.y_param = y_param

Expand All @@ -743,14 +817,14 @@ def compute_fn(meta, log, fit_methods, nsamp=1e4):
flux[2] = flux[0]-flux[2]
meta.spectrum_median.append(flux[0])
meta.spectrum_err.append(flux[1:])

# Convert the lists to an array
meta.spectrum_median = np.array(meta.spectrum_median)
if meta.fitter == 'lsq':
meta.spectrum_err = np.ones((2, meta.nspecchan))*np.nan
else:
meta.spectrum_err = np.array(meta.spectrum_err).T

return meta


Expand Down Expand Up @@ -921,7 +995,7 @@ def save_table(meta, log):
meta.spectrum_err)

transit_latex_table(meta, log)

return


Expand Down Expand Up @@ -1064,7 +1138,7 @@ def transit_latex_table(meta, log):
# Trim off the leading space
y_unit = meta.y_label_unit[1:]
# Need to make sure to escape % with \ for LaTeX
if (meta.y_label_unit.count(r'\%') !=
if (meta.y_label_unit.count(r'\%') !=
meta.y_label_unit.count('%')):
y_unit = y_unit.replace('%', r'\%')
out += "\\colhead{("+xunit+")} & \\colhead{"+y_unit+"} &"
Expand Down Expand Up @@ -1113,7 +1187,7 @@ def transit_latex_table(meta, log):
# End the table
out += "\\enddata\n"
out += "\\end{deluxetable}"

# Save the table as a txt file
meta.tab_filename_s6_latex = meta.tab_filename_s6[:-4]+'_LaTeX.txt'
with open(meta.tab_filename_s6_latex, 'w') as file:
Expand Down

0 comments on commit 3389d9b

Please sign in to comment.