diff --git a/exotic/api/colab.py b/exotic/api/colab.py index 7041b499..33b27e81 100644 --- a/exotic/api/colab.py +++ b/exotic/api/colab.py @@ -1,39 +1,78 @@ - -# If the user presses enter to run the sample data, download sample data if needed and -# put it into a sample-data directory at the top level of the user's Gdrive. Count -# the .fits files (images) and .json files (inits files) in the directory entered -# by the user (or in the sample-data directory if the user pressed enter). If -# there are at least 20 .fits files, assume this is a directory of images and display -# the first one in the series. If there is exactly one inits file in the directory, -# show the specified target and comp coords so that the user can check these against -# the displayed image. Otherwise, prompt for target / comp coords and make an inits -# file based on those (save this new inits file in the folder with the output files -# so that the student can consult it later). Finally, run EXOTIC with the newly-made +# ########################################################################### # +# Copyright (c) 2019-2020, California Institute of Technology. +# All rights reserved. Based on Government Sponsored Research under +# contracts NNN12AA01C, NAS7-1407 and/or NAS7-03001. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# 3. Neither the name of the California Institute of +# Technology (Caltech), its operating division the Jet Propulsion +# Laboratory (JPL), the National Aeronautics and Space +# Administration (NASA), nor the names of its contributors may be +# used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE CALIFORNIA +# INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +# ########################################################################### # +# EXOplanet Transit Interpretation Code (EXOTIC) +# # NOTE: See companion file version.py for version info. +# ########################################################################### # +# IMPORTANT RUNTIME INFORMATION ABOUT DATA ACCESS AND STORAGE +# +# If the user presses enter to run the sample data, download sample data if +# needed and put it into a sample-data directory at the top level of the +# user's Gdrive. Count the .fits files (images) and .json files (inits files) +# in the directory entered by the user (or in the sample-data directory if the +# user pressed enter). If there are at least 20 .fits files, assume this is a +# directory of images and display the first one in the series. If there is +# exactly one inits file in the directory, show the specified target and comp +# coords so that the user can check these against the displayed image. +# Otherwise, prompt for target / comp coords and make an inits file based on +# those (save this new inits file in the folder with the output files so that +# the student can consult it later). Finally, run EXOTIC with the newly-made # or pre-existing inits file, plus any other inits files in the directory. - +# ######################################################### -from IPython.display import display, HTML +from astropy.io import fits from astropy.time import Time from barycorrpy import utc_tdb -import numpy as np -from io import BytesIO -from astropy.io import fits -from scipy.ndimage import label -from bokeh.plotting import figure, output_file, show +# import bokeh.io +# from bokeh.io import output_notebook from bokeh.palettes import Viridis256 -from bokeh.models import ColorBar, LinearColorMapper, LogColorMapper, LogTicker -from bokeh.models import BoxZoomTool,WheelZoomTool,ResetTool,HoverTool,PanTool,FreehandDrawTool -#import bokeh.io -#from bokeh.io import output_notebook -from pprint import pprint -#from IPython.display import Image -#from ipywidgets import widgets, HBox -from skimage.transform import rescale, resize, downscale_local_mean -#import copy +from bokeh.plotting import figure, output_file, show +from bokeh.models import BoxZoomTool, ColorBar, FreehandDrawTool, HoverTool, LinearColorMapper, LogColorMapper, \ + LogTicker, PanTool, ResetTool, WheelZoomTool +# import copy +from io import BytesIO +from IPython.display import display, HTML +# from IPython.display import Image +# from ipywidgets import widgets, HBox +import json +import numpy as np import os +from pprint import pprint import re -import json -#import subprocess +from scipy.ndimage import label +from skimage.transform import rescale, resize, downscale_local_mean +# import subprocess import time diff --git a/exotic/api/elca.py b/exotic/api/elca.py index eef283f4..35768b2d 100644 --- a/exotic/api/elca.py +++ b/exotic/api/elca.py @@ -35,21 +35,19 @@ # EXOplanet Transit Interpretation Code (EXOTIC) # # NOTE: See companion file version.py for version info. # ########################################################################### # -# ########################################################################### # # Exoplanet light curve analysis # # Fit an exoplanet transit model to time series data. # ########################################################################### # - from astropy.time import Time import copy from itertools import cycle import matplotlib.pyplot as plt import numpy as np +from pylightcurve.models.exoplanet_lc import transit as pytransit from scipy import spatial from scipy.optimize import least_squares from scipy.signal import savgol_filter - try: from ultranest import ReactiveNestedSampler except ImportError: @@ -57,13 +55,12 @@ import dynesty.plotting from dynesty.utils import resample_equal from scipy.stats import gaussian_kde + try: from plotting import corner -except: +except ImportError: from .plotting import corner -from pylightcurve.models.exoplanet_lc import transit as pytransit - def weightedflux(flux, gw, nearest): return np.sum(flux[nearest] * gw, axis=-1) diff --git a/exotic/api/ephemeris.py b/exotic/api/ephemeris.py index 81349b1b..f9bff13b 100644 --- a/exotic/api/ephemeris.py +++ b/exotic/api/ephemeris.py @@ -41,16 +41,31 @@ # a periodogram analysis. # # ########################################################################### # - -import numpy as np -from itertools import cycle -import statsmodels.api as sm -import matplotlib.pyplot as plt -from exotic.api.plotting import corner -from ultranest import ReactiveNestedSampler from astropy.timeseries import LombScargle from astropy import units as u import copy +from itertools import cycle +import matplotlib.pyplot as plt +import numpy as np +import rebound +import statsmodels.api as sm +try: + from ultranest import ReactiveNestedSampler +except ImportError: + import dynesty + import dynesty.plotting + from dynesty.utils import resample_equal + from scipy.stats import gaussian_kde + +try: + from nested_linear_fitter import linear_fitter, non_linear_fitter +except ImportError: + from .nested_linear_fitter import linear_fitter, non_linear_fitter +try: + from plotting import corner +except ImportError: + from .plotting import corner + class ephemeris_fitter(object): @@ -910,363 +925,117 @@ def plot_triangle(self): return fig - -def empty_data(N): # orbit parameters in Nbody simulation - return { - 'x':np.zeros(N), - 'y':np.zeros(N), - 'z':np.zeros(N), - 'P':np.zeros(N), - 'a':np.zeros(N), - 'e':np.zeros(N), - 'inc':np.zeros(N), - 'Omega':np.zeros(N), - 'omega':np.zeros(N), - 'M':np.zeros(N), - } - -def integrate(sim, objects, Ndays, Noutputs): - # ps is now an array of pointers and will change as the simulation runs - ps = sim.particles - - times = np.linspace(0., Ndays, Noutputs) # units of days - - pdata = [empty_data(Noutputs) for i in range(len(objects)-1)] - star = {'x':np.zeros(Noutputs),'y':np.zeros(Noutputs),'z':np.zeros(Noutputs) } - - # run simulation time steps - for i,time in enumerate(times): - sim.integrate(time) - - # record data - for k in star.keys(): - star[k][i] = getattr(ps[0], k) - - for j in range(1,len(objects)): - for k in pdata[j-1].keys(): - pdata[j-1][k][i] = getattr(ps[j],k) - - sim_data = ({ - 'pdata':pdata, - 'star':star, - 'times':times, - 'objects':objects, - 'dt': Noutputs/(Ndays*24*60*60) # conversion factor to get seconds for RV semi-amp - }) - - return sim_data - -def generate(objects, Ndays=None, Noutputs=None): - # create rebound simulation - # for object parameters see: - # https://rebound.readthedocs.io/en/latest/_modules/rebound/particle.html - sim = rebound.Simulation() - - # sim.integrator = "tes" - # sim.ri_tes.dq_max = 1e-2 - # sim.ri_tes.recti_per_orbit = 1.61803398875 - # sim.ri_tes.epsilon = 1e-6 - - # sim.integrator = "whfast" - # sim.ri_whfast.corrector = 5 - # sim.ri_whfast.safe_mode = 0 - sim.units = ('day', 'AU', 'Msun') - for i in range(len(objects)): - sim.add( **objects[i] ) - sim.move_to_com() - - if Ndays and Noutputs: - return integrate(sim, objects, Ndays, Noutputs) - else: - return sim - -def find_zero(t1,dx1, t2,dx2): - # find zero with linear interpolation - m = (dx2-dx1)/(t2-t1) - T0 = -dx1/m + t1 - return T0 - -def transit_times(xp,xs,times): - # check for sign change in position difference - dx = xp-xs - tt = [] - for i in range(1,len(dx)): - if dx[i-1] >= 0 and dx[i] <= 0: - tt.append( find_zero(times[i-1],dx[i-1], times[i],dx[i]) ) - return np.array(tt) - -def TTV(epochs, tt): - N = len(epochs) - A = np.vstack([np.ones(N), epochs]).T - b, m = np.linalg.lstsq(A, tt, rcond=None)[0] - ttv = (tt-m*np.array(epochs)-b) - return [ttv,m,b] - -class nbody_fitter(): - - def __init__(self, data, prior=None, bounds=None, verbose=True): - self.data = data - self.bounds = bounds - self.prior = prior - self.verbose = verbose - self.fit_nested() - - def fit_nested(self): - - # set up some arrays for mapping sampler output - freekeys = [] - boundarray = [] - for i,planet in enumerate(self.bounds): - for bound in planet: - freekeys.append(f"{i}_{bound}") - boundarray.append(planet[bound]) - - # find min and max time for simulation - min_time = np.min(self.data[1]['Tc']) - max_time = np.max(self.data[1]['Tc']) - # todo extend for multiplanet systems - - sim_time = (max_time-min_time)*1.05 - # TODO extend for multiplanet systems - Tc_norm = self.data[1]['Tc'] - min_time # normalize the data to the first observation - self.orbit = np.rint(Tc_norm / self.prior[1]['P']).astype(int) # number of orbits since first observation (rounded to nearest integer) - - # numpify - boundarray = np.array(boundarray) - bounddiff = np.diff(boundarray,1).reshape(-1) - - # create queue and save simulation to - def loglike(pars): - chi2 = 0 - - # set parameters - for i,par in enumerate(pars): - idx,key = freekeys[i].split('_') - idx = int(idx) - if key == 'tmid': - continue - # this dict goes to REBOUND and needs specific keys - self.prior[idx][key] = par - - # take difference between data and simulation - # run N-body simulation - sim_data = generate(self.prior, sim_time, int(sim_time*24)) # uses linspace behind the scenes - - # json structure with analytics of interest from the simulation - # ttv_data = analyze(sim_data) # slow - - sim_shift = 0 - # loop over planets, check for data - for i,planet in enumerate(self.prior): - if self.data[i]: - # compute transit times from N-body simulation - Tc_sim = transit_times( sim_data['pdata'][i-1]['x'], sim_data['star']['x'], sim_data['times'] ) - - # derive an offset in time from the first planet - if i-1==0: - sim_shift = Tc_sim.min() - - # shift the first mid-transit in simulation to observation - Tc_sim -= sim_shift # first orbit has tmid at 0 - - # scale Tc_sim to data - try: - residual = self.data[i]['Tc'] - Tc_sim[self.orbit] - except: - #import pdb; pdb.set_trace - #ttv,m,b = TTV(np.arange(Tc_sim.shape[0]), Tc_sim) - #ttv1,m1,b1 = TTV(self.orbit, self.data[i]['Tc']) - # could be unstable orbit or not enough data - # switch to average error and clip by max epoch? - print(self.prior) - chi2 += -1e6 - continue - - Tc_sim += residual.mean() - - # take difference between data and simulation - try: - chi2 += -0.5*np.sum(((self.data[i]['Tc'] - Tc_sim[self.orbit])/self.data[i]['Tc_err'])**2) - except: - chi2 += -1e6 - print(self.prior) - # usually unstable orbit - - return chi2 - - - def prior_transform(upars): - return (boundarray[:,0] + bounddiff*upars) - - if self.verbose: - self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=1e5) - else: - self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=1e5, show_status=self.verbose, viz_callback=self.verbose) - - self.errors = {} - self.quantiles = {} - self.parameters = copy.deepcopy(self.prior) - - #for i, key in enumerate(freekeys): - # self.parameters[key] = self.results['maximum_likelihood']['point'][i] - # self.errors[key] = self.results['posterior']['stdev'][i] - # self.quantiles[key] = [ - # self.results['posterior']['errlo'][i], - # self.results['posterior']['errup'][i]] - - -#def main(): if __name__ == "__main__": Tc = np.array([ # measured mid-transit times #Exoplanet Watch transit times array -2457025.7867,2457347.7655,2457694.8263,2458112.8495,2458492.6454,2459532.7834,2459604.81036,2459604.8137,2459614.6365,2459616.8209, -2459651.7466,2459903.8588,2459914.7783,2459915.8684,2459925.6921,2459939.8793,2459949.7047,2459959.5249,2459962.8037,2459973.7129, -2459974.798,2459986.8057,2459994.4489,2460009.7312,2458867.01587,2459501.1295, - -#ExoClock transit times array -2454508.977,2454515.525,2454801.477,2454840.769,2455123.447,2455140.91,2455147.459,2455147.459,2455148.551,2455148.552,2455158.372, -2455158.373,2455159.464,2455159.467,2455160.556,2455163.831,2455172.562,2455192.205,2455209.669,2455210.762,2455215.129,2455230.407, -2455238.046,2455254.419,2455265.331,2455494.53,2455494.531,2455509.81,2455510.902,2455530.547,2455542.553,2455548.011,2455565.472, -2455566.563,2455575.296,2455575.297,2455578.569,2455589.483,2455590.576,2455598.216,2455600.397,2455600.398,2455600.398,2455601.49, -2455601.49,2455601.49,2455603.673,2455612.404,2455623.318,2455623.319,2455624.41,2455624.41,2455663.702,2455842.693,2455842.694, -2455876.528,2455887.442,2455887.442,2455887.442,2455888.533,2455888.534,2455888.534,2455889.625,2455890.716,2455901.63,2455903.814, -2455903.814,2455920.185,2455923.459,2455926.733,2455939.829,2455945.287,2455946.378,2455946.379,2455947.47,2455948.561,2455951.835, -2455952.927,2455959.475,2455960.567,2455970.39,2455971.481,2455981.304,2455982.395,2455983.487,2455984.578,2455985.67,2455996.584, -2456000.95,2456005.315,2456006.406,2456014.046,2456175.577,2456245.427,2456249.794,2456273.805,2456282.536,2456284.719,2456297.816, -2456302.182,2456305.455,2456319.644,2456328.376,2456329.467,2456604.505,2456605.596,2456606.688,2456607.779,2456629.607,2456630.699, -2456654.71,2456659.076,2456662.35,2456663.441,2456664.533,2456674.356,2456677.63,2456688.544,2456694.002,2456703.824,2456711.464, -2456719.104,2456721.287,2456722.378,2456986.502,2457010.513,2457012.696,2457033.434,2457045.438,2457046.53,2457059.627,2457060.718, -2457067.267,2457068.358,2457103.284,2457345.579,2457368.5,2457378.321,2457390.327,2457391.418,2457426.343,2457426.344,2457427.435, -2457451.446,2457474.364,2457486.372,2457671.913,2457691.559,2457703.564,2457706.838,2457726.484,2457727.575,2457760.318,2457765.775, -2457766.866,2457772.324,2457773.415,2457776.689,2457786.512,2457788.695,2457800.7,2457808.34,2457809.432,2457809.434,2457810.523, -2457832.348,2457832.351,2458026.624,2458050.635,2458060.459,2458072.462,2458073.555,2458074.647,2458077.921,2458123.76,2458124.852, -2458132.491,2458134.675,2458136.858,2458155.41,2458155.412,2458156.503,2458159.778,2458166.326,2458178.331,2458214.347,2458411.895, -2458467.557,2458471.923,2458478.472,2458479.564,2458490.477,2458494.843,2458501.39,2458501.392,2458506.848,2458536.315,2458537.408, -2458871.382,2458880.113,2458883.384,2458893.209,2458903.032,2459088.575,2459148.602,2459156.242,2459157.334,2459190.076,2459194.441, -2459194.443,2459196.623,2459197.717,2459204.264, - -#ETD mid transit times array -2454508.976854,2454836.403393,2454840.768603,2454840.771193,2454908.437992,2454931.358182,2455164.923956,2455164.924316,2455172.561816,2455172.562786, -2455198.756735,2455230.406730,2455254.418870,2455256.596033,2455265.334053,2455506.533285,2455509.808915,2455519.632074,2455532.729964,2455541.461084, -2455576.387242,2455577.475172,2455577.476782,2455577.477112,2455578.568422,2455593.850192,2455600.398282,2455601.489211,2455603.671161,2455612.400801, -2455612.401321,2455615.675911,2455647.328760,2455857.973783,2455857.973783,2455867.797832,2455873.253612,2455877.618482,2455877.619162,2455889.623302, -2455899.448441,2455899.451251,2455900.539781,2455904.905271,2455905.995241,2455912.544561,2455912.544561,2455913.636131,2455914.726741,2455922.370121, -2455928.913321,2455936.552560,2455936.561520,2455937.645430,2455938.738730,2455946.383030,2455947.468480,2455949.650900,2455949.650900,2455950.742510, -2455957.287600,2455957.294400,2455957.296780,2455957.296780,2455961.659130,2455962.749930,2455970.392369,2455982.394919,2455993.309869,2455993.313119, -2455994.398719,2455994.399359,2456003.129689,2456003.130169,2456005.315328,2456222.509695,2456246.516835,2456249.793845,2456269.425195,2456271.623265, -2456296.724384,2456304.364534,2456304.365804,2456304.367454,2456304.367784,2456329.469294,2456364.392544,2456584.860173,2456595.772813,2456596.866383, -2456604.503073,2456608.869564,2456616.511354,2456626.333324,2456627.429034,2456628.515004,2456628.517074,2456630.699754,2456639.430344,2456639.432894, -2456640.522554,2456641.610274,2456650.341584,2456663.443614,2456688.544734,2456722.372535,2456722.382615,2456734.383505,2456734.388675,2456746.385515, -2456746.391885,2456927.546138,2456950.485119,2456960.306869,2456963.578439,2456963.578959,2456963.582479,2456986.501140,2456998.508180,2456998.508180, -2457032.333311,2457033.433561,2457069.448912,2457080.367572,2457080.367572,2457092.368353,2457344.486723,2457344.488693,2457345.578813,2457345.579833, -2457356.493733,2457357.585013,2457357.586493,2457368.498964,2457368.499654,2457369.589474,2457379.414424,2457379.414424,2457414.334726,2457426.343356, -2457451.447237,2457474.364078,2457474.364348,2457486.371248,2457668.637786,2457726.484138,2457749.405049,2457750.494949,2457751.585649,2457751.586939, -2457757.043281,2457758.134851,2457760.317471,2457760.318171,2457761.408951,2457770.140361,2457772.327372,2457773.415162,2457774.506522,2457809.431273, -2457809.435043,2457832.350874,2458087.745843,2458096.474134,2458109.571734,2458131.400315,2458396.615162,2458396.617562,2458411.896052,2458455.554593, -2458457.734253,2458479.568743,2458489.383454,2458489.387804,2458490.477904,2458490.482824,2458501.391674,2458501.392014,2458501.396544,2458513.402254, -2458537.407634,2458538.507734,2458837.550156,2458848.461686,2458859.375956,2458860.465876,2458860.469616,2458871.379136,2458871.382556,2458872.470796, -2458883.387165,2458883.389555,2458884.478175,2458884.478425,2458884.479455,2458906.306925,2458930.332215,2459147.510251,2459171.522921,2459172.608741, -2459172.617851,2459183.531910,2459193.352480,2459194.442340,2459194.443610,2459196.623950,2459197.719250,2459220.632699,2459242.461499,2459265.384258, -2459265.384318,2459505.495609,2459553.514627,2459553.515137,2459564.427746,2459565.520376,2459566.615006,2459576.438096]) + 2457025.7867,2457347.7655,2457694.8263,2458112.8495,2458492.6454,2459532.7834,2459604.81036,2459604.8137,2459614.6365,2459616.8209, + 2459651.7466,2459903.8588,2459914.7783,2459915.8684,2459925.6921,2459939.8793,2459949.7047,2459959.5249,2459962.8037,2459973.7129, + 2459974.798,2459986.8057,2459994.4489,2460009.7312,2458867.01587,2459501.1295, + + #ExoClock transit times array + 2454508.977,2454515.525,2454801.477,2454840.769,2455123.447,2455140.91,2455147.459,2455147.459,2455148.551,2455148.552,2455158.372, + 2455158.373,2455159.464,2455159.467,2455160.556,2455163.831,2455172.562,2455192.205,2455209.669,2455210.762,2455215.129,2455230.407, + 2455238.046,2455254.419,2455265.331,2455494.53,2455494.531,2455509.81,2455510.902,2455530.547,2455542.553,2455548.011,2455565.472, + 2455566.563,2455575.296,2455575.297,2455578.569,2455589.483,2455590.576,2455598.216,2455600.397,2455600.398,2455600.398,2455601.49, + 2455601.49,2455601.49,2455603.673,2455612.404,2455623.318,2455623.319,2455624.41,2455624.41,2455663.702,2455842.693,2455842.694, + 2455876.528,2455887.442,2455887.442,2455887.442,2455888.533,2455888.534,2455888.534,2455889.625,2455890.716,2455901.63,2455903.814, + 2455903.814,2455920.185,2455923.459,2455926.733,2455939.829,2455945.287,2455946.378,2455946.379,2455947.47,2455948.561,2455951.835, + 2455952.927,2455959.475,2455960.567,2455970.39,2455971.481,2455981.304,2455982.395,2455983.487,2455984.578,2455985.67,2455996.584, + 2456000.95,2456005.315,2456006.406,2456014.046,2456175.577,2456245.427,2456249.794,2456273.805,2456282.536,2456284.719,2456297.816, + 2456302.182,2456305.455,2456319.644,2456328.376,2456329.467,2456604.505,2456605.596,2456606.688,2456607.779,2456629.607,2456630.699, + 2456654.71,2456659.076,2456662.35,2456663.441,2456664.533,2456674.356,2456677.63,2456688.544,2456694.002,2456703.824,2456711.464, + 2456719.104,2456721.287,2456722.378,2456986.502,2457010.513,2457012.696,2457033.434,2457045.438,2457046.53,2457059.627,2457060.718, + 2457067.267,2457068.358,2457103.284,2457345.579,2457368.5,2457378.321,2457390.327,2457391.418,2457426.343,2457426.344,2457427.435, + 2457451.446,2457474.364,2457486.372,2457671.913,2457691.559,2457703.564,2457706.838,2457726.484,2457727.575,2457760.318,2457765.775, + 2457766.866,2457772.324,2457773.415,2457776.689,2457786.512,2457788.695,2457800.7,2457808.34,2457809.432,2457809.434,2457810.523, + 2457832.348,2457832.351,2458026.624,2458050.635,2458060.459,2458072.462,2458073.555,2458074.647,2458077.921,2458123.76,2458124.852, + 2458132.491,2458134.675,2458136.858,2458155.41,2458155.412,2458156.503,2458159.778,2458166.326,2458178.331,2458214.347,2458411.895, + 2458467.557,2458471.923,2458478.472,2458479.564,2458490.477,2458494.843,2458501.39,2458501.392,2458506.848,2458536.315,2458537.408, + 2458871.382,2458880.113,2458883.384,2458893.209,2458903.032,2459088.575,2459148.602,2459156.242,2459157.334,2459190.076,2459194.441, + 2459194.443,2459196.623,2459197.717,2459204.264, + + #ETD mid transit times array + 2454508.976854,2454836.403393,2454840.768603,2454840.771193,2454908.437992,2454931.358182,2455164.923956,2455164.924316,2455172.561816,2455172.562786, + 2455198.756735,2455230.406730,2455254.418870,2455256.596033,2455265.334053,2455506.533285,2455509.808915,2455519.632074,2455532.729964,2455541.461084, + 2455576.387242,2455577.475172,2455577.476782,2455577.477112,2455578.568422,2455593.850192,2455600.398282,2455601.489211,2455603.671161,2455612.400801, + 2455612.401321,2455615.675911,2455647.328760,2455857.973783,2455857.973783,2455867.797832,2455873.253612,2455877.618482,2455877.619162,2455889.623302, + 2455899.448441,2455899.451251,2455900.539781,2455904.905271,2455905.995241,2455912.544561,2455912.544561,2455913.636131,2455914.726741,2455922.370121, + 2455928.913321,2455936.552560,2455936.561520,2455937.645430,2455938.738730,2455946.383030,2455947.468480,2455949.650900,2455949.650900,2455950.742510, + 2455957.287600,2455957.294400,2455957.296780,2455957.296780,2455961.659130,2455962.749930,2455970.392369,2455982.394919,2455993.309869,2455993.313119, + 2455994.398719,2455994.399359,2456003.129689,2456003.130169,2456005.315328,2456222.509695,2456246.516835,2456249.793845,2456269.425195,2456271.623265, + 2456296.724384,2456304.364534,2456304.365804,2456304.367454,2456304.367784,2456329.469294,2456364.392544,2456584.860173,2456595.772813,2456596.866383, + 2456604.503073,2456608.869564,2456616.511354,2456626.333324,2456627.429034,2456628.515004,2456628.517074,2456630.699754,2456639.430344,2456639.432894, + 2456640.522554,2456641.610274,2456650.341584,2456663.443614,2456688.544734,2456722.372535,2456722.382615,2456734.383505,2456734.388675,2456746.385515, + 2456746.391885,2456927.546138,2456950.485119,2456960.306869,2456963.578439,2456963.578959,2456963.582479,2456986.501140,2456998.508180,2456998.508180, + 2457032.333311,2457033.433561,2457069.448912,2457080.367572,2457080.367572,2457092.368353,2457344.486723,2457344.488693,2457345.578813,2457345.579833, + 2457356.493733,2457357.585013,2457357.586493,2457368.498964,2457368.499654,2457369.589474,2457379.414424,2457379.414424,2457414.334726,2457426.343356, + 2457451.447237,2457474.364078,2457474.364348,2457486.371248,2457668.637786,2457726.484138,2457749.405049,2457750.494949,2457751.585649,2457751.586939, + 2457757.043281,2457758.134851,2457760.317471,2457760.318171,2457761.408951,2457770.140361,2457772.327372,2457773.415162,2457774.506522,2457809.431273, + 2457809.435043,2457832.350874,2458087.745843,2458096.474134,2458109.571734,2458131.400315,2458396.615162,2458396.617562,2458411.896052,2458455.554593, + 2458457.734253,2458479.568743,2458489.383454,2458489.387804,2458490.477904,2458490.482824,2458501.391674,2458501.392014,2458501.396544,2458513.402254, + 2458537.407634,2458538.507734,2458837.550156,2458848.461686,2458859.375956,2458860.465876,2458860.469616,2458871.379136,2458871.382556,2458872.470796, + 2458883.387165,2458883.389555,2458884.478175,2458884.478425,2458884.479455,2458906.306925,2458930.332215,2459147.510251,2459171.522921,2459172.608741, + 2459172.617851,2459183.531910,2459193.352480,2459194.442340,2459194.443610,2459196.623950,2459197.719250,2459220.632699,2459242.461499,2459265.384258, + 2459265.384318,2459505.495609,2459553.514627,2459553.515137,2459564.427746,2459565.520376,2459566.615006,2459576.438096]) Tc_error = np.array([ 0.0039,0.0058,0.0054,0.0022,0.0065,0.0024,0.00091,0.00097,0.0022,0.0018, -0.002,0.0021,0.0065,0.0036,0.0071,0.0031,0.0044,0.0021,0.0028,0.0023, -0.0039,0.0071,0.0015,0.0055,0.00092,0.00092, - -#ExoClock error array -0.0002,0.00013,0.00041,0.00047,0.0007,0.000417,0.0009,0.00042,0.001,0.0013,0.00057, -0.0013,0.00095,0.00061,0.0007,0.000324,0.00044,0.0005,0.000463,0.000405,0.00082,0.00011, -0.00066,0.00014,0.0016,0.00048,0.00063,0.00037,0.000313,0.00077,0.00029,0.00024,0.0011, -0.00022,0.00063,0.00079,0.0008,0.00094,0.001,0.00036,0.00083,0.00026,0.00078,0.000289, -0.00017,0.00039,0.001,0.00088,0.0009,0.00054,0.00086,0.0014,0.00063,0.00093,0.0013, -0.0003,0.00038,0.00062,0.00075,0.00059,0.0003,0.00038,0.0009,0.00023,0.00083,0.00093, -0.000324,0.00046,0.0002,0.00092,0.0019,0.00079,0.00036,0.00034,0.00022,0.00028,0.00011, -0.0001,0.00018,0.0004,0.00029,0.00029,0.00094,0.00047,0.00029,0.000324,0.000417,0.00037, -0.0004,0.00038,0.0004,0.00023,0.00033,0.00033,0.000394,0.000301,0.0003,0.000301,0.000301, -0.00046,0.00026,0.000382,0.00027,0.00029,0.0002,0.0003,0.00034,0.000706,0.00019,0.00043, -0.000336,0.00034,0.00019,0.00019,0.00032,0.00028,0.000324,0.00041,0.00029,0.00029,0.00026, -0.00034,0.00034,0.00046,0.00043,0.00039,0.000486,0.0005,0.00049,0.00049,0.000347,0.000359, -0.00022,0.00021,0.0003,0.00042,0.0004,0.0013,0.00034,0.00033,0.00055,0.0006,0.00023,0.00021, -0.0007,0.0013,0.00035,0.00025,0.00034,0.00037,0.00028,0.00023,0.0006,0.00028,0.00039, -0.00024,0.00022,0.00029,0.00026,0.00048,0.00032,0.0004,0.00018,0.0009,0.00021,0.0006, -0.0006,0.00056,0.00023,0.0003,0.0003,0.00022,0.00034,0.00028,0.00027,0.00035,0.00031, -0.00032,0.00033,0.0005,0.00031,0.00032,0.00091,0.00034,0.00038,0.0017,0.0004,0.0005, -0.00026,0.0006,0.0006,0.0008,0.0003,0.0009,0.0003,0.00044,0.0008,0.0007,0.0009, -0.0003,0.0007,0.0005,0.0003,0.0013,0.0007,0.0003,0.0004,0.0003,0.0012,0.0006, -0.0005,0.0013,0.0004, - -#ETD error array -0.000200,0.000600,0.000470,0.001000,0.001000,0.000980,0.001490,0.001290,0.000440,0.000140, -0.001410,0.000110,0.000140,0.000590,0.001290,0.001120,0.000870,0.000700,0.000350,0.000500, -0.001350,0.000380,0.000690,0.000850,0.000820,0.000470,0.000420,0.001090,0.000940,0.000830, -0.000780,0.000540,0.001730,0.000710,0.000710,0.000750,0.001260,0.000920,0.001290,0.000730, -0.000630,0.000570,0.000380,0.001030,0.001350,0.000570,0.000570,0.000780,0.000470,0.000900, -0.000730,0.000910,0.001230,0.001190,0.000720,0.000770,0.001020,0.000590,0.000590,0.000660, -0.001100,0.001040,0.000570,0.000570,0.001070,0.001320,0.000860,0.001160,0.000600,0.000760, -0.000680,0.000760,0.000630,0.000600,0.000440,0.000810,0.000740,0.000670,0.000900,0.000550, -0.000520,0.001460,0.000890,0.001560,0.000580,0.001640,0.001170,0.000510,0.000960,0.000510, -0.000920,0.000710,0.000900,0.000510,0.001050,0.000970,0.000880,0.000440,0.000740,0.000680, -0.000820,0.000800,0.001040,0.000760,0.000540,0.000890,0.001080,0.001120,0.000730,0.001390, -0.001410,0.001640,0.000740,0.000570,0.000930,0.000890,0.000610,0.000510,0.001450,0.001450, -0.001130,0.000460,0.000510,0.001190,0.001190,0.000600,0.000650,0.001070,0.001090,0.000490, -0.000610,0.000520,0.000430,0.000580,0.000460,0.000340,0.000890,0.000890,0.001310,0.000650, -0.000860,0.000770,0.000550,0.001350,0.000820,0.000550,0.000840,0.000530,0.000940,0.000720, -0.000370,0.000520,0.000670,0.001030,0.000630,0.000330,0.000990,0.000550,0.000620,0.000720, -0.000940,0.000580,0.000350,0.000570,0.001380,0.000640,0.001180,0.000700,0.000700,0.000710, -0.000550,0.000580,0.000680,0.001030,0.000860,0.000850,0.000300,0.000760,0.000680,0.000820, -0.000610,0.001450,0.000730,0.000700,0.001350,0.000820,0.000670,0.000940,0.000490,0.001130, -0.000540,0.000540,0.000700,0.000790,0.000840,0.000600,0.000520,0.000730,0.000640,0.001020, -0.000780,0.000610,0.001330,0.000770,0.000610,0.000520,0.001130,0.001130,0.000530,0.000780, -0.000420,0.001250,0.000380,0.000720,0.000860,0.000470,0.000950,0.000540]) + 0.002,0.0021,0.0065,0.0036,0.0071,0.0031,0.0044,0.0021,0.0028,0.0023, + 0.0039,0.0071,0.0015,0.0055,0.00092,0.00092, + + #ExoClock error array + 0.0002,0.00013,0.00041,0.00047,0.0007,0.000417,0.0009,0.00042,0.001,0.0013,0.00057, + 0.0013,0.00095,0.00061,0.0007,0.000324,0.00044,0.0005,0.000463,0.000405,0.00082,0.00011, + 0.00066,0.00014,0.0016,0.00048,0.00063,0.00037,0.000313,0.00077,0.00029,0.00024,0.0011, + 0.00022,0.00063,0.00079,0.0008,0.00094,0.001,0.00036,0.00083,0.00026,0.00078,0.000289, + 0.00017,0.00039,0.001,0.00088,0.0009,0.00054,0.00086,0.0014,0.00063,0.00093,0.0013, + 0.0003,0.00038,0.00062,0.00075,0.00059,0.0003,0.00038,0.0009,0.00023,0.00083,0.00093, + 0.000324,0.00046,0.0002,0.00092,0.0019,0.00079,0.00036,0.00034,0.00022,0.00028,0.00011, + 0.0001,0.00018,0.0004,0.00029,0.00029,0.00094,0.00047,0.00029,0.000324,0.000417,0.00037, + 0.0004,0.00038,0.0004,0.00023,0.00033,0.00033,0.000394,0.000301,0.0003,0.000301,0.000301, + 0.00046,0.00026,0.000382,0.00027,0.00029,0.0002,0.0003,0.00034,0.000706,0.00019,0.00043, + 0.000336,0.00034,0.00019,0.00019,0.00032,0.00028,0.000324,0.00041,0.00029,0.00029,0.00026, + 0.00034,0.00034,0.00046,0.00043,0.00039,0.000486,0.0005,0.00049,0.00049,0.000347,0.000359, + 0.00022,0.00021,0.0003,0.00042,0.0004,0.0013,0.00034,0.00033,0.00055,0.0006,0.00023,0.00021, + 0.0007,0.0013,0.00035,0.00025,0.00034,0.00037,0.00028,0.00023,0.0006,0.00028,0.00039, + 0.00024,0.00022,0.00029,0.00026,0.00048,0.00032,0.0004,0.00018,0.0009,0.00021,0.0006, + 0.0006,0.00056,0.00023,0.0003,0.0003,0.00022,0.00034,0.00028,0.00027,0.00035,0.00031, + 0.00032,0.00033,0.0005,0.00031,0.00032,0.00091,0.00034,0.00038,0.0017,0.0004,0.0005, + 0.00026,0.0006,0.0006,0.0008,0.0003,0.0009,0.0003,0.00044,0.0008,0.0007,0.0009, + 0.0003,0.0007,0.0005,0.0003,0.0013,0.0007,0.0003,0.0004,0.0003,0.0012,0.0006, + 0.0005,0.0013,0.0004, + + #ETD error array + 0.000200,0.000600,0.000470,0.001000,0.001000,0.000980,0.001490,0.001290,0.000440,0.000140, + 0.001410,0.000110,0.000140,0.000590,0.001290,0.001120,0.000870,0.000700,0.000350,0.000500, + 0.001350,0.000380,0.000690,0.000850,0.000820,0.000470,0.000420,0.001090,0.000940,0.000830, + 0.000780,0.000540,0.001730,0.000710,0.000710,0.000750,0.001260,0.000920,0.001290,0.000730, + 0.000630,0.000570,0.000380,0.001030,0.001350,0.000570,0.000570,0.000780,0.000470,0.000900, + 0.000730,0.000910,0.001230,0.001190,0.000720,0.000770,0.001020,0.000590,0.000590,0.000660, + 0.001100,0.001040,0.000570,0.000570,0.001070,0.001320,0.000860,0.001160,0.000600,0.000760, + 0.000680,0.000760,0.000630,0.000600,0.000440,0.000810,0.000740,0.000670,0.000900,0.000550, + 0.000520,0.001460,0.000890,0.001560,0.000580,0.001640,0.001170,0.000510,0.000960,0.000510, + 0.000920,0.000710,0.000900,0.000510,0.001050,0.000970,0.000880,0.000440,0.000740,0.000680, + 0.000820,0.000800,0.001040,0.000760,0.000540,0.000890,0.001080,0.001120,0.000730,0.001390, + 0.001410,0.001640,0.000740,0.000570,0.000930,0.000890,0.000610,0.000510,0.001450,0.001450, + 0.001130,0.000460,0.000510,0.001190,0.001190,0.000600,0.000650,0.001070,0.001090,0.000490, + 0.000610,0.000520,0.000430,0.000580,0.000460,0.000340,0.000890,0.000890,0.001310,0.000650, + 0.000860,0.000770,0.000550,0.001350,0.000820,0.000550,0.000840,0.000530,0.000940,0.000720, + 0.000370,0.000520,0.000670,0.001030,0.000630,0.000330,0.000990,0.000550,0.000620,0.000720, + 0.000940,0.000580,0.000350,0.000570,0.001380,0.000640,0.001180,0.000700,0.000700,0.000710, + 0.000550,0.000580,0.000680,0.001030,0.000860,0.000850,0.000300,0.000760,0.000680,0.000820, + 0.000610,0.001450,0.000730,0.000700,0.001350,0.000820,0.000670,0.000940,0.000490,0.001130, + 0.000540,0.000540,0.000700,0.000790,0.000840,0.000600,0.000520,0.000730,0.000640,0.001020, + 0.000780,0.000610,0.001330,0.000770,0.000610,0.000520,0.001130,0.001130,0.000530,0.000780, + 0.000420,0.001250,0.000380,0.000720,0.000860,0.000470,0.000950,0.000540]) labels = np.full(len(Tc), 'EpW') - # # Tc = np.array([ # measured mid-transit times - # # 2459150.837905, 2459524.851045, - # # 2459546.613126, 2459565.643663, 2459584.690470, 2459584.686476, - # # 2459909.739104, 2459957.337739, 2459169.880602, 2458416.424861, - # # 2458428.664794, 2459145.400124, 2458430.025044, 2459164.440695, - # # 2458415.064846, 2458435.465269, 2458412.344658, 2459150.840070, - # # 2459160.360691, 2458431.384897, 2459146.760284, 2459154.920516, - # # 2458417.784945, 2459161.720466, 2459167.160761, 2458427.304939, - # # 2458413.705181, 2459163.080605, 2459153.560229, 2459168.520861, - # # 2458425.945068, 2459148.120269, 2458434.105274, 2458432.745423, - # # 2459152.200558, 2459165.800918, 2459159.000645, 2459149.480289, - # # 2455870.450027, 2456663.347570, 2457420.884390, 2457658.888900]) - - # # Tc_error = np.array([ - # # 0.001674, 0.000715, - # # 0.000758, 0.001560, 0.000371, 0.001651, - # # 0.000955, 0.000176, 0.000154, 0.000357, - # # 0.000381, 0.000141, 0.000362, 0.000149, - # # 0.000336, 0.000368, 0.000379, 0.000153, - # # 0.000153, 0.000349, 0.000149, 0.000146, - # # 0.000385, 0.000146, 0.000153, 0.000360, - # # 0.000356, 0.000147, 0.000147, 0.000146, - # # 0.000363, 0.000142, 0.000357, 0.000368, - # # 0.000160, 0.000160, 0.000151, 0.000160, - # # 0.000140, 0.000120, 0.000800, 0.000140]) - - # # labels for a legend - # labels = np.array(['TESS', 'TESS', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar']) P = 1.0914203 # orbital period for your target Tc_norm = Tc - Tc.min() # normalize the data to the first observation - # print(Tc_norm) orbit = np.rint(Tc_norm / P) # number of orbits since first observation (rounded to nearest integer) - # print(orbit) # make a n x 2 matrix with 1's in the first column and values of orbit in the second A = np.vstack([np.ones(len(Tc)), orbit]).T @@ -1302,17 +1071,17 @@ def prior_transform(upars): # min and max values to search between for fitting bounds = { - 'm': [P - 0.1, P + 0.1], # orbital period - 'b': [intercept - 0.1, intercept + 0.1] # mid-transit time + 'P': [slope - 0.1, slope + 0.1], # orbital period + 'T0': [intercept - 0.1, intercept + 0.1] # mid-transit time } # used to plot red overlay in O-C figure prior = { - 'm': [slope, slope_std_dev], # value from WLS (replace with literature value) - 'b': [intercept, intercept_std_dev] # value from WLS (replace with literature value) + 'P': [slope, slope_std_dev], # value from WLS (replace with literature value) + 'T0': [intercept, intercept_std_dev] # value from WLS (replace with literature value) } - lf = linear_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels) + lf = ephemeris_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels) lf.plot_triangle() plt.subplots_adjust(top=0.9, hspace=0.2, wspace=0.2) @@ -1348,7 +1117,7 @@ def prior_transform(upars): 'dPdN': [0, 0.0001] # best guess } - nlf = non_linear_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels) + nlf = decay_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels) nlf.plot_triangle() plt.subplots_adjust(top=0.9, hspace=0.2, wspace=0.2) @@ -1367,64 +1136,4 @@ def prior_transform(upars): for key in nlf.parameters: print(f"Parameter {key} = {nlf.parameters[key]:.2e} +- {nlf.errors[key]:.2e}") - # TODO extend plot, BIC Values - - dude() - # mearth = u.M_earth.to(u.kg).value - msun = u.M_sun.to(u.kg) - mjup = u.M_jup.to(u.kg) - - # Parameters for N-body Retrieval - # read more about adding particles: https://rebound.readthedocs.io/en/latest/addingparticles/ - nbody_prior = [ - {'m':1.12}, # star - - {'m':0.28*mjup/msun, - 'P':lf.parameters['P'], # inner planet - 'inc':3.14159/2, - 'e':0, - 'omega':0}, - - {'m':0.988*mjup/msun, - 'P':lf.parameters['P']*3, # outer planet - 'inc':3.14159/2, - 'e':0, - 'omega':0 }, - ] - - - # specify data to fit - data = [ - {}, # data for star (e.g. RV) - {'Tc':Tc, 'Tc_err':Tc_error}, # data for inner planet (e.g. Mid-transit times) - {} # data for outer planet (e.g. Mid-transit times) - ] - - - # set up where to look for a solution - bounds = [ - - {}, # no bounds on star parameters - - { # bounds for inner planet - 'P': [nbody_prior[1]['P']-0.025, nbody_prior[1]['P']+0.025], # based on solution from linear fit - # mid-transit is automatically adjusted in the fitter class - }, - - { # bounds for outer planet - 'P':[lf.parameters['P']*1.75, lf.parameters['P']*4.25], # orbital period [day] - 'm':[0.1*mjup/msun,1.5*mjup/msun], # mass [msun] - #'e':[0,0.05], # eccentricity - #'omega':[-20*np.pi/180, 20*np.pi/180] # arg of periastron [radians] - } - ] - - # run the fitter - nfit = nbody_fitter(data, nbody_prior, bounds) - # todo how could I phase fold the data to just shorter nbody simulations? - - print(nfit.parameters) - print(nfit.errors) - - -# main() + # TODO BIC Values \ No newline at end of file diff --git a/exotic/api/ew.py b/exotic/api/ew.py index c179e2cb..bf1071ab 100644 --- a/exotic/api/ew.py +++ b/exotic/api/ew.py @@ -1,10 +1,48 @@ -import os +# ########################################################################### # +# Copyright (c) 2019-2020, California Institute of Technology. +# All rights reserved. Based on Government Sponsored Research under +# contracts NNN12AA01C, NAS7-1407 and/or NAS7-03001. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# 3. Neither the name of the California Institute of +# Technology (Caltech), its operating division the Jet Propulsion +# Laboratory (JPL), the National Aeronautics and Space +# Administration (NASA), nor the names of its contributors may be +# used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE CALIFORNIA +# INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +# ########################################################################### # +# EXOplanet Transit Interpretation Code (EXOTIC) +# # NOTE: See companion file version.py for version info. +# ########################################################################### # import json -import urllib import numpy as np +import os +import urllib base_uri = "https://exoplanets.nasa.gov/watch_results/" + class ExoplanetWatch(): url = 'https://exoplanets.nasa.gov/api/v1/exoplanet_watch_results/?per_page=-1&order=planet_name+asc' def __init__(self) -> None: @@ -50,6 +88,7 @@ def get(self, target): print(f"{target} not in list") return {} + class ExoplanetWatchTarget(): def __init__(self,result) -> None: self.raw_result = result @@ -92,6 +131,7 @@ def check4floats(dictionary): pass return dictionary + class ExoplanetWatchObservation(): def __init__(self, name, observation) -> None: self.raw_observation = observation @@ -130,6 +170,7 @@ def get_data(self): jdata = json.loads(r.read().decode(r.info().get_param('charset') or 'utf-8')) return np.array(jdata[1:],dtype=float).T + def translate_keys(rdict): """ Translates the keys to a compatible format for EXOTIC/ELCA @@ -165,11 +206,11 @@ def translate_keys(rdict): lc_pars[k] = np.nan return lc_pars -if __name__ == "__main__": +if __name__ == "__main__": EW = ExoplanetWatch() print(EW.target_list) result = EW.get('WASP-33 b') print(result) time, flux, fluxerr, airmass, airmasscorr = result.observations[0].get_data() - print(result.observations[0].parameters) \ No newline at end of file + print(result.observations[0].parameters) diff --git a/exotic/api/gael_ld.py b/exotic/api/gael_ld.py index 2d53b162..1b42aa10 100644 --- a/exotic/api/gael_ld.py +++ b/exotic/api/gael_ld.py @@ -35,14 +35,12 @@ # EXOplanet Transit Interpretation Code (EXOTIC) # # NOTE: See companion file version.py for version info. # ########################################################################### # - -import logging -import numpy as np -import lmfit as lm -import matplotlib.pyplot as plt import ldtk -from ldtk import LDPSetCreator, BoxcarFilter from ldtk.ldmodel import LinearModel, QuadraticModel, NonlinearModel +import lmfit as lm +import logging +import matplotlib.pyplot as plt +import numpy as np log = logging.getLogger(__name__) @@ -92,9 +90,9 @@ def createldgrid(minmu, maxmu, orbp, munm = 1e3*np.array(avmu[loweri:upperi]) munmmin = 1e3*np.array(minmu[loweri:upperi]) munmmax = 1e3*np.array(maxmu[loweri:upperi]) - filters = [BoxcarFilter(str(mue), mun, mux) + filters = [ldtk.BoxcarFilter(str(mue), mun, mux) for mue, mun, mux in zip(munm, munmmin, munmmax)] - sc = LDPSetCreator(teff=(tstar, terr), logg=(loggstar, loggerr), + sc = ldtk.LDPSetCreator(teff=(tstar, terr), logg=(loggstar, loggerr), z=(fehstar, feherr), filters=filters) ps = sc.create_profiles(nsamples=int(1e4)) itpfail = False @@ -105,7 +103,7 @@ def createldgrid(minmu, maxmu, orbp, nfail = 1e0 while itpfail: nfail *= 2 - sc = LDPSetCreator(teff=(tstar, nfail*terr), logg=(loggstar, loggerr), + sc = ldtk.LDPSetCreator(teff=(tstar, nfail*terr), logg=(loggstar, loggerr), z=(fehstar, feherr), filters=filters) ps = sc.create_profiles(nsamples=int(1e4)) itpfail = False @@ -254,4 +252,3 @@ def nl_ldx(params, x, data=None, weights=None): if weights is None: return data - model return (data - model)/weights - diff --git a/exotic/api/joint_fitter.py b/exotic/api/joint_fitter.py index a746ca57..db2df69f 100644 --- a/exotic/api/joint_fitter.py +++ b/exotic/api/joint_fitter.py @@ -1,14 +1,60 @@ -import numpy as np -from scipy import stats +# ########################################################################### # +# Copyright (c) 2019-2020, California Institute of Technology. +# All rights reserved. Based on Government Sponsored Research under +# contracts NNN12AA01C, NAS7-1407 and/or NAS7-03001. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# 3. Neither the name of the California Institute of +# Technology (Caltech), its operating division the Jet Propulsion +# Laboratory (JPL), the National Aeronautics and Space +# Administration (NASA), nor the names of its contributors may be +# used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE CALIFORNIA +# INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +# ########################################################################### # +# EXOplanet Transit Interpretation Code (EXOTIC) +# # NOTE: See companion file version.py for version info. +# ########################################################################### # +from astropy import constants as const +from astropy import units as u from copy import deepcopy from itertools import cycle import matplotlib.pyplot as plt -from astropy import units as u -from astropy import constants as const +import numpy as np from pylightcurve.models.exoplanet_lc import eclipse_mid_time, transit_flux_drop -from ultranest import ReactiveNestedSampler - -from exotic.api.elca import glc_fitter, lc_fitter +from scipy import stats +try: + from ultranest import ReactiveNestedSampler +except ImportError: + import dynesty + import dynesty.plotting + from dynesty.utils import resample_equal + from scipy.stats import gaussian_kde + +try: + from elca import glc_fitter, lc_fitter +except ImportError: + from .elca import glc_fitter, lc_fitter AU = const.au.to(u.m).value Mjup = const.M_jup.to(u.kg).value @@ -16,8 +62,9 @@ Rsun = const.R_sun.to(u.m).value Grav = const.G.to(u.m**3/u.kg/u.day**2).value + def planet_orbit(period, sma_over_rs, eccentricity, inclination, periastron, mid_time, time_array, ww=0, mu=1, W=0): - # please see original: https://github.com/ucl-exoplanets/pylightcurve/blob/master/pylightcurve/models/exoplanet_lc.py + # see original @ https://github.com/ucl-exoplanets/pylightcurve/blob/master/pylightcurve/models/exoplanet_lc.py inclination = inclination * np.pi / 180.0 periastron = periastron * np.pi / 180.0 ww = ww * np.pi / 180.0 @@ -107,7 +154,7 @@ def phasecurve(times, values): values['ecc'], values['inc'], values['omega'], values['tmid'], times, method='quad', precision=3) c0 = edepth - values['c1'] - values['c3'] - brightness = 1 + c0 + values['c1']*np.cos(2*pi*(times-tme)/values['per']) + values['c2']*np.sin(2*pi*(times-tme)/values['per']) + values['c3']*np.cos(4*pi*(times-tme)/values['per']) + values['c4']*np.sin(4*pi*(times-tme)/values['per']) + brightness = 1 + c0 + values['c1']*np.cos(2*np.pi*(times-tme)/values['per']) + values['c2']*np.sin(2*np.pi*(times-tme)/values['per']) + values['c3']*np.cos(4*np.pi*(times-tme)/values['per']) + values['c4']*np.sin(4*np.pi*(times-tme)/values['per']) emask = np.floor(emodel) return (brightness*(emask) + (edepth+emodel)*(1-emask))*tmodel diff --git a/exotic/api/ld.py b/exotic/api/ld.py index cff8a5a7..5e939f1e 100644 --- a/exotic/api/ld.py +++ b/exotic/api/ld.py @@ -1,15 +1,52 @@ +# ########################################################################### # +# Copyright (c) 2019-2020, California Institute of Technology. +# All rights reserved. Based on Government Sponsored Research under +# contracts NNN12AA01C, NAS7-1407 and/or NAS7-03001. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# 3. Neither the name of the California Institute of +# Technology (Caltech), its operating division the Jet Propulsion +# Laboratory (JPL), the National Aeronautics and Space +# Administration (NASA), nor the names of its contributors may be +# used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE CALIFORNIA +# INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +# ########################################################################### # +# EXOplanet Transit Interpretation Code (EXOTIC) +# # NOTE: See companion file version.py for version info. +# ########################################################################### # import logging import numpy as np import re try: - from .gael_ld import createldgrid -except ImportError: from gael_ld import createldgrid -try: - from .filters import fwhm, fwhm_alias, fwhm_names_nonspecific except ImportError: + from .gael_ld import createldgrid +try: from filters import fwhm, fwhm_alias, fwhm_names_nonspecific +except ImportError: + from .filters import fwhm, fwhm_alias, fwhm_names_nonspecific log = logging.getLogger(__name__) ld_re_punct = r'[\W]' diff --git a/exotic/api/nea.py b/exotic/api/nea.py index 3cc35b5f..d78bb406 100644 --- a/exotic/api/nea.py +++ b/exotic/api/nea.py @@ -35,7 +35,9 @@ # EXOplanet Transit Interpretation Code (EXOTIC) # # NOTE: See companion file version.py for version info. # ########################################################################### # - +import astropy.constants as const +from astropy.coordinates import SkyCoord +import astropy.units as u from io import StringIO import json import numpy as np @@ -43,9 +45,6 @@ import pandas import requests import time -import astropy.units as u -import astropy.constants as const -from astropy.coordinates import SkyCoord from tenacity import retry, retry_if_exception_type, stop_after_attempt, \ wait_exponential diff --git a/exotic/api/nested_linear_fitter.py b/exotic/api/nested_linear_fitter.py index 1152dab1..229d4c8a 100644 --- a/exotic/api/nested_linear_fitter.py +++ b/exotic/api/nested_linear_fitter.py @@ -35,25 +35,32 @@ # EXOplanet Transit Interpretation Code (EXOTIC) # # NOTE: See companion file version.py for version info. # ########################################################################### # -# ########################################################################### # -# Various functions for fitting a linear model to data, including nested +# Various functions for fitting a linear model to data, including nested # sampling, linear least squares. Includes residual plotting, posteriors and # a periodogram analysis. # # ########################################################################### # - -import numpy as np -from itertools import cycle -import statsmodels.api as sm -import matplotlib.pyplot as plt -from exotic.api.plotting import corner -from ultranest import ReactiveNestedSampler +from astropy import constants as const from astropy.timeseries import LombScargle from astropy import units as u -from astropy import constants as const from collections import OrderedDict -import rebound import copy +from itertools import cycle +import matplotlib.pyplot as plt +import numpy as np +import statsmodels.api as sm +try: + from ultranest import ReactiveNestedSampler +except ImportError: + import dynesty + import dynesty.plotting + from dynesty.utils import resample_equal + from scipy.stats import gaussian_kde + +try: + from plotting import corner +except ImportError: + from .plotting import corner class linear_fitter(object): @@ -906,210 +913,6 @@ def plot_triangle(self): return fig -def empty_data(N): # orbit parameters in Nbody simulation - return { - 'x':np.zeros(N), - 'y':np.zeros(N), - 'z':np.zeros(N), - 'P':np.zeros(N), - 'a':np.zeros(N), - 'e':np.zeros(N), - 'inc':np.zeros(N), - 'Omega':np.zeros(N), - 'omega':np.zeros(N), - 'M':np.zeros(N), - } - -def integrate(sim, objects, Ndays, Noutputs): - # ps is now an array of pointers and will change as the simulation runs - ps = sim.particles - - times = np.linspace(0., Ndays, Noutputs) # units of days - - pdata = [empty_data(Noutputs) for i in range(len(objects)-1)] - star = {'x':np.zeros(Noutputs),'y':np.zeros(Noutputs),'z':np.zeros(Noutputs) } - - # run simulation time steps - for i,time in enumerate(times): - sim.integrate(time) - - # record data - for k in star.keys(): - star[k][i] = getattr(ps[0], k) - - for j in range(1,len(objects)): - for k in pdata[j-1].keys(): - pdata[j-1][k][i] = getattr(ps[j],k) - - sim_data = ({ - 'pdata':pdata, - 'star':star, - 'times':times, - 'objects':objects, - 'dt': Noutputs/(Ndays*24*60*60) # conversion factor to get seconds for RV semi-amp - }) - - return sim_data - -def generate(objects, Ndays=None, Noutputs=None): - # create rebound simulation - # for object parameters see: - # https://rebound.readthedocs.io/en/latest/_modules/rebound/particle.html - sim = rebound.Simulation() - - # sim.integrator = "tes" - # sim.ri_tes.dq_max = 1e-2 - # sim.ri_tes.recti_per_orbit = 1.61803398875 - # sim.ri_tes.epsilon = 1e-6 - - # sim.integrator = "whfast" - # sim.ri_whfast.corrector = 5 - # sim.ri_whfast.safe_mode = 0 - sim.units = ('day', 'AU', 'Msun') - for i in range(len(objects)): - sim.add( **objects[i] ) - sim.move_to_com() - - if Ndays and Noutputs: - return integrate(sim, objects, Ndays, Noutputs) - else: - return sim - -def find_zero(t1,dx1, t2,dx2): - # find zero with linear interpolation - m = (dx2-dx1)/(t2-t1) - t0 = -dx1/m + t1 - return t0 - -def transit_times(xp,xs,times): - # check for sign change in position difference - dx = xp-xs - tt = [] - for i in range(1,len(dx)): - if dx[i-1] >= 0 and dx[i] <= 0: - tt.append( find_zero(times[i-1],dx[i-1], times[i],dx[i]) ) - return np.array(tt) - -def TTV(epochs, tt): - N = len(epochs) - A = np.vstack([np.ones(N), epochs]).T - b, m = np.linalg.lstsq(A, tt, rcond=None)[0] - ttv = (tt-m*np.array(epochs)-b) - return [ttv,m,b] - -class nbody_fitter(): - - def __init__(self, data, prior=None, bounds=None, verbose=True): - self.data = data - self.bounds = bounds - self.prior = prior - self.verbose = verbose - self.fit_nested() - - def fit_nested(self): - - # set up some arrays for mapping sampler output - freekeys = [] - boundarray = [] - for i,planet in enumerate(self.bounds): - for bound in planet: - freekeys.append(f"{i}_{bound}") - boundarray.append(planet[bound]) - - # find min and max time for simulation - min_time = np.min(self.data[1]['Tc']) - max_time = np.max(self.data[1]['Tc']) - # todo extend for multiplanet systems - - sim_time = (max_time-min_time)*1.05 - # TODO extend for multiplanet systems - Tc_norm = self.data[1]['Tc'] - min_time # normalize the data to the first observation - self.orbit = np.rint(Tc_norm / self.prior[1]['P']).astype(int) # number of orbits since first observation (rounded to nearest integer) - - # numpify - boundarray = np.array(boundarray) - bounddiff = np.diff(boundarray,1).reshape(-1) - - # create queue and save simulation to - def loglike(pars): - chi2 = 0 - - # set parameters - for i,par in enumerate(pars): - idx,key = freekeys[i].split('_') - idx = int(idx) - if key == 'tmid': - continue - # this dict goes to REBOUND and needs specific keys - self.prior[idx][key] = par - - # take difference between data and simulation - # run N-body simulation - sim_data = generate(self.prior, sim_time, int(sim_time*24)) # uses linspace behind the scenes - - # json structure with analytics of interest from the simulation - # ttv_data = analyze(sim_data) # slow - - sim_shift = 0 - # loop over planets, check for data - for i,planet in enumerate(self.prior): - if self.data[i]: - # compute transit times from N-body simulation - Tc_sim = transit_times( sim_data['pdata'][i-1]['x'], sim_data['star']['x'], sim_data['times'] ) - - # derive an offset in time from the first planet - if i-1==0: - sim_shift = Tc_sim.min() - - # shift the first mid-transit in simulation to observation - Tc_sim -= sim_shift # first orbit has tmid at 0 - - # scale Tc_sim to data - try: - residual = self.data[i]['Tc'] - Tc_sim[self.orbit] - except: - #import pdb; pdb.set_trace - #ttv,m,b = TTV(np.arange(Tc_sim.shape[0]), Tc_sim) - #ttv1,m1,b1 = TTV(self.orbit, self.data[i]['Tc']) - # could be unstable orbit or not enough data - # switch to average error and clip by max epoch? - print(self.prior) - chi2 += -1e6 - continue - - Tc_sim += residual.mean() - - # take difference between data and simulation - try: - chi2 += -0.5*np.sum(((self.data[i]['Tc'] - Tc_sim[self.orbit])/self.data[i]['Tc_err'])**2) - except: - chi2 += -1e6 - print(self.prior) - # usually unstable orbit - - return chi2 - - - def prior_transform(upars): - return (boundarray[:,0] + bounddiff*upars) - - if self.verbose: - self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=1e5) - else: - self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=1e5, show_status=self.verbose, viz_callback=self.verbose) - - self.errors = {} - self.quantiles = {} - self.parameters = copy.deepcopy(self.prior) - - #for i, key in enumerate(freekeys): - # self.parameters[key] = self.results['maximum_likelihood']['point'][i] - # self.errors[key] = self.results['posterior']['stdev'][i] - # self.quantiles[key] = [ - # self.results['posterior']['errlo'][i], - # self.results['posterior']['errup'][i]] - - #def main(): if __name__ == "__main__": @@ -1360,66 +1163,4 @@ def prior_transform(upars): # for each free key print parameter and error for key in nlf.parameters: - print(f"Parameter {key} = {nlf.parameters[key]:.2e} +- {nlf.errors[key]:.2e}") - - # TODO extend plot, BIC Values - - dude() - # mearth = u.M_earth.to(u.kg).value - msun = u.M_sun.to(u.kg) - mjup = u.M_jup.to(u.kg) - - # Parameters for N-body Retrieval - # read more about adding particles: https://rebound.readthedocs.io/en/latest/addingparticles/ - nbody_prior = [ - {'m':1.12}, # star - - {'m':0.28*mjup/msun, - 'P':lf.parameters['P'], # inner planet - 'inc':3.14159/2, - 'e':0, - 'omega':0}, - - {'m':0.988*mjup/msun, - 'P':lf.parameters['P']*3, # outer planet - 'inc':3.14159/2, - 'e':0, - 'omega':0 }, - ] - - - # specify data to fit - data = [ - {}, # data for star (e.g. RV) - {'Tc':Tc, 'Tc_err':Tc_error}, # data for inner planet (e.g. Mid-transit times) - {} # data for outer planet (e.g. Mid-transit times) - ] - - - # set up where to look for a solution - bounds = [ - - {}, # no bounds on star parameters - - { # bounds for inner planet - 'P': [nbody_prior[1]['P']-0.025, nbody_prior[1]['P']+0.025], # based on solution from linear fit - # mid-transit is automatically adjusted in the fitter class - }, - - { # bounds for outer planet - 'P':[lf.parameters['P']*1.75, lf.parameters['P']*4.25], # orbital period [day] - 'm':[0.1*mjup/msun,1.5*mjup/msun], # mass [msun] - #'e':[0,0.05], # eccentricity - #'omega':[-20*np.pi/180, 20*np.pi/180] # arg of periastron [radians] - } - ] - - # run the fitter - nfit = nbody_fitter(data, nbody_prior, bounds) - # todo how could I phase fold the data to just shorter nbody simulations? - - print(nfit.parameters) - print(nfit.errors) - - -# main() + print(f"Parameter {key} = {nlf.parameters[key]:.2e} +- {nlf.errors[key]:.2e}") \ No newline at end of file diff --git a/exotic/api/output_aavso.py b/exotic/api/output_aavso.py index 4d953752..a24ff580 100644 --- a/exotic/api/output_aavso.py +++ b/exotic/api/output_aavso.py @@ -1,19 +1,57 @@ -#:-) +# ########################################################################### # +# Copyright (c) 2019-2020, California Institute of Technology. +# All rights reserved. Based on Government Sponsored Research under +# contracts NNN12AA01C, NAS7-1407 and/or NAS7-03001. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# 3. Neither the name of the California Institute of +# Technology (Caltech), its operating division the Jet Propulsion +# Laboratory (JPL), the National Aeronautics and Space +# Administration (NASA), nor the names of its contributors may be +# used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE CALIFORNIA +# INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +# ########################################################################### # +# EXOplanet Transit Interpretation Code (EXOTIC) +# # NOTE: See companion file version.py for version info. +# ########################################################################### # +import hashlib from json import dump, dumps -from tkinter import NONE from numpy import mean, median, std -from pathlib import Path import os +from pathlib import Path import re +from tkinter import NONE + try: - from utils import round_to_2 -except ImportError: from .utils import round_to_2 -try: - from version import __version__ except ImportError: + from utils import round_to_2 +try: from .version import __version__ -import hashlib +except ImportError: + from version import __version__ + class OutputFiles: def __init__(self, fit, p_dict, i_dict, planetdir): @@ -133,8 +171,7 @@ def aavso(self, airmasses, ld0, ld1, ld2, ld3, tmidstr): f"{round(self.fit.dataerr[aavsoC], 7)},{0.0}," f"{1.000}\n") -#*********** FOR CSV ************ - +# *********** FOR CSV ************ def aavso_csv(self, airmasses, ld0, ld1, ld2, ld3,tmidstr): priors_dict, filter_dict, results_dict = aavso_dicts(self.p_dict, self.fit, self.i_dict, self.durs, ld0, ld1, ld2, ld3) @@ -192,7 +229,6 @@ def aavso_csv(self, airmasses, ld0, ld1, ld2, ld3,tmidstr): f"{round(self.fit.airmass_model[aavsoC], 7)}\n") - def aavso_dicts(planet_dict, fit, i_dict, durs, ld0, ld1, ld2, ld3): priors = { 'Period': { @@ -277,4 +313,4 @@ def aavso_dicts(planet_dict, fit, i_dict, durs, ld0, ld1, ld2, ld3): 'units': "degrees" } - return priors, filter_type, results \ No newline at end of file + return priors, filter_type, results diff --git a/exotic/api/plate_solution.py b/exotic/api/plate_solution.py index 48a45c24..93642b01 100644 --- a/exotic/api/plate_solution.py +++ b/exotic/api/plate_solution.py @@ -35,11 +35,10 @@ # EXOplanet Transit Interpretation Code (EXOTIC) # # NOTE: See companion file version.py for version info. # ########################################################################### # - -import requests +from astropy.io.fits import PrimaryHDU, getdata, getheader from json import dumps from pathlib import Path -from astropy.io.fits import PrimaryHDU, getdata, getheader +import requests from tenacity import retry, retry_if_exception_type, retry_if_result, \ stop_after_attempt, wait_exponential diff --git a/exotic/api/plotting.py b/exotic/api/plotting.py index b464d0a9..4d26ab67 100644 --- a/exotic/api/plotting.py +++ b/exotic/api/plotting.py @@ -1,21 +1,60 @@ -import json -import numpy as np -from io import BytesIO +# ########################################################################### # +# Copyright (c) 2019-2020, California Institute of Technology. +# All rights reserved. Based on Government Sponsored Research under +# contracts NNN12AA01C, NAS7-1407 and/or NAS7-03001. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# 3. Neither the name of the California Institute of +# Technology (Caltech), its operating division the Jet Propulsion +# Laboratory (JPL), the National Aeronautics and Space +# Administration (NASA), nor the names of its contributors may be +# used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE CALIFORNIA +# INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +# ########################################################################### # +# EXOplanet Transit Interpretation Code (EXOTIC) +# # NOTE: See companion file version.py for version info. +# ########################################################################### # from astropy.io import fits # from astroscrappy import detect_cosmics -from scipy.ndimage import label -from skimage.transform import downscale_local_mean -from bokeh.plotting import figure, output_file, show -from bokeh.palettes import Viridis256 -from bokeh.models import ColorBar, LinearColorMapper, LogColorMapper, LogTicker -from bokeh.models import BoxZoomTool,WheelZoomTool,ResetTool,HoverTool,PanTool,FreehandDrawTool from bokeh.io import output_notebook -from pprint import pprint +from bokeh.models import BoxZoomTool, ColorBar, FreehandDrawTool, HoverTool, LinearColorMapper, LogColorMapper, \ + LogTicker, PanTool, ResetTool, WheelZoomTool +from bokeh.palettes import Viridis256 +from bokeh.plotting import figure, output_file, show +from io import BytesIO +import json +import logging import matplotlib.pyplot as plt -from matplotlib.ticker import MaxNLocator, NullLocator -from matplotlib.ticker import ScalarFormatter +from matplotlib.ticker import MaxNLocator, NullLocator, ScalarFormatter +import numpy as np +from pprint import pprint from scipy.interpolate import griddata -from scipy.ndimage import gaussian_filter +from scipy.ndimage import label, gaussian_filter +from skimage.transform import downscale_local_mean + +log = logging.getLogger(__name__) + def plot_image(filename, save=False, bg_min=60, bg_max=99): @@ -159,8 +198,7 @@ def corner(xs, bins=20, range=None, weights=None, color="k", hist_bin_factor=1, # Parse the parameter ranges. if range is None: if "extents" in hist2d_kwargs: - logging.warn("Deprecated keyword argument 'extents'. " - "Use 'range' instead.") + log.warn("Deprecated keyword argument 'extents'. Use 'range' instead.") range = hist2d_kwargs.pop("extents") else: range = [[x.min(), x.max()] for x in xs] diff --git a/exotic/api/rv_fitter.py b/exotic/api/rv_fitter.py index cd424a50..f4fdc6ee 100644 --- a/exotic/api/rv_fitter.py +++ b/exotic/api/rv_fitter.py @@ -1,12 +1,60 @@ +# ########################################################################### # +# Copyright (c) 2019-2020, California Institute of Technology. +# All rights reserved. Based on Government Sponsored Research under +# contracts NNN12AA01C, NAS7-1407 and/or NAS7-03001. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in +# the documentation and/or other materials provided with the +# distribution. +# 3. Neither the name of the California Institute of +# Technology (Caltech), its operating division the Jet Propulsion +# Laboratory (JPL), the National Aeronautics and Space +# Administration (NASA), nor the names of its contributors may be +# used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE CALIFORNIA +# INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +# ########################################################################### # +# EXOplanet Transit Interpretation Code (EXOTIC) +# # NOTE: See companion file version.py for version info. +# ########################################################################### # +from astropy import constants as const +from astropy import units as u import copy from itertools import cycle -import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt -from astropy import units as u -from astropy import constants as const -from ultranest import ReactiveNestedSampler -from exotic.api.elca import lc_fitter +import numpy as np +try: + from ultranest import ReactiveNestedSampler +except ImportError: + import dynesty + import dynesty.plotting + from dynesty.utils import resample_equal + from scipy.stats import gaussian_kde + +try: + from elca import lc_fitter +except ImportError: + from .elca import lc_fitter + Mjup = const.M_jup.to(u.kg).value Msun = const.M_sun.to(u.kg).value diff --git a/exotic/version.py b/exotic/version.py index 32206102..6a157dcb 100644 --- a/exotic/version.py +++ b/exotic/version.py @@ -1 +1 @@ -__version__ = '3.2.3' +__version__ = '3.3.0' diff --git a/requirements.txt b/requirements.txt index f2ca7781..7fa4feac 100644 --- a/requirements.txt +++ b/requirements.txt @@ -18,6 +18,7 @@ pylightcurve~=4.0.1 python_dateutil~=2.8 python_version pyvo~=1.4 +rebound~=3.28.3 requests~=2.31 scipy~=1.11.1 scikit-image~=0.21