Merge pull request #1228 from rzellem/nbody

Use-case document + fix to multi-lc fitter for citisens

examples/README.md

@@ -0,0 +1,68 @@
+# EXOTIC Use Cases
+
+This repository contains examples of how to use the EXOTIC software to perform a variety of tasks related to exoplanet transit science. The package is designed to be used with FITS images, photometric data, radial velocity data, and ephemeris data. The examples below are organized by the type of data used in the analysis.
+
+
+## [Programmatic Access to Exoplanet Watch Results](Exoplanet_Watch_API.ipynb) 
+
+Download parameters that are derived from photometric data and light curves for a given target. Over 350 targets and 5000 publicly available light curves from ground-based and space-based telescopes.
+
+![](epw_results.png)
+
+## Acknowledgements
+
+If you use any Exoplanet Watch data or in your publication, you are required to include the observers of those data as co-authors on your paper. To get in touch with your anonymous observer, contact the [AAVSO](https://www.aavso.org/) with their observer code.
+
+If you make use of Exoplanet Watch in your work, please cite the papers [Zellem et al. 2020](https://ui.adsabs.harvard.edu/abs/2020PASP..132e4401Z/abstract) and [Pearson et al. 2022](https://ui.adsabs.harvard.edu/abs/2022AJ....164..178P/abstract) and include the following standard acknowledgment in any published material that makes use of Exoplanet Watch data: **“This publication makes use of data products from Exoplanet Watch, a citizen science project managed by NASA's Jet Propulsion Laboratory on behalf of NASA's Universe of Learning. This work is supported by NASA under award number NNX16AC65A to the Space Telescope Science Institute, in partnership with Caltech/IPAC, Center for Astrophysics|Harvard & Smithsonian, and NASA Jet Propulsion Laboratory."**
+
+## [Fit a transit light curve](single_transit/transit_fit_example.py)
+
+Start from a photometric timeseries and derive transit parameters like planetary radius, inclination and mid-transit time. Optimized in a Bayesian framework ([Ultranest](https://johannesbuchner.github.io/UltraNest/index.html)) with posteriors to assess degeneracies and uncertainties.
+
+![](single_transit/bestfit.png)
+![](single_transit/triangle.png)
+
+
+## Fit multiple light curves simultaneously with shared and individual parameters
+
+- [Simultaneous airmass detrending](multiple_transit/Multiple_Lightcurve_fit.ipynb) (more robust but takes much longer)
+
+- [Airmass detrending prior to simultaneous fit](multiple_transit/Multiple_Lightcurve_Fit_Detrended.ipynb)
+
+The notebooks above are also compatible with TESS data! Just don't include an `a2` parameter for airmass detrending in the local bounds.
+
+![](multiple_transit/glc_fit.png)
+
+![](multiple_transit/glc_mosaic.png)
+
+
+## [TESS light curve generation](tess.py)
+
+![](tess/2458415_06_wasp-77ab_lightcurve.png)
+
+## [Ephemeris fitting](ephemeris/fit_ephemeris.py)
+- Observed - Calculated plot with colors coded to data source
+
+![](ephemeris/oc.png)
+
+- [Orbital Decay](ephemeris/fit_decay.py)
+
+![](ephemeris/decay_oc.png)
+![](ephemeris/decay_posterior.png)
+
+- Periodogram for transit timing search with up to two orders
+
+![](ephemeris/periodogram.png)
+
+## [N-body interpretation of periodogram](nbody/README.md)
+
+![](nbody/report.png)
+
+## [Radial Velocity](radial_velocity/rv_example.py)
+
+![](radial_velocity/RV_bestfit.png)
+
+
+## [Joint Fit of Transit Photometry, Radial Velocity, and Ephemeris data (transit/eclipse times)](joint_rv_transit/joint_fit.py)
+
+![](joint_rv_transit/joint_posterior.png)

examples/ephemeris/fit_decay.py

@@ -0,0 +1,181 @@
+import numpy as np
+import statsmodels.api as sm
+from exotic.api.ephemeris import decay_fitter
+import matplotlib.pyplot as plt
+
+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
+    ])
+
+    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
+    ])
+
+    # labels for a legend
+    labels = np.array(['Data']*len(Tc)) 
+    # TODO give individual labels to each dataset
+
+    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
+
+    # perform the weighted least squares regression
+    res = sm.WLS(Tc, A, weights=1.0 / Tc_error ** 2).fit()
+    # use sm.WLS for weighted LS, sm.OLS for ordinary LS, or sm.GLS for general LS
+
+    params = res.params  # retrieve the slope and intercept of the fit from res
+    std_dev = np.sqrt(np.diagonal(res.normalized_cov_params))
+
+    slope = params[1]
+    slope_std_dev = std_dev[1]
+    intercept = params[0]
+    intercept_std_dev = std_dev[0]
+
+    # 3 sigma clip based on residuals
+    calculated = orbit * slope + intercept
+    residuals = (Tc - calculated) / Tc_error
+    mask = np.abs(residuals) < 3
+    Tc = Tc[mask]
+    Tc_error = Tc_error[mask]
+    labels = labels[mask]
+
+    # print(res.summary())
+    # print("Params =",params)
+    # print("Error matrix =",res.normalized_cov_params)
+    # print("Standard Deviations =",std_dev)
+
+    print("Weighted Linear Least Squares Solution")
+    print("T0 =", intercept, "+-", intercept_std_dev)
+    print("P =", slope, "+-", slope_std_dev)
+
+    # min and max values to search between for fitting
+    bounds = {
+        'P': [P - 0.1, P + 0.1],                   # orbital period [day]
+        'T0': [intercept - 0.1, intercept + 0.1],  # mid-transit time [BJD]
+        'dPdN': [-0.00001, 0.00001],               # orbital period decay rate [day/epoch]
+    }
+
+    # used to plot red overlay in O-C figure
+    prior = {
+        '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 = decay_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels)
+
+    lf.plot_triangle()
+    plt.subplots_adjust(top=0.9, hspace=0.2, wspace=0.2)
+    plt.savefig("decay_posterior.png")
+    plt.close()
+
+    fig, ax = lf.plot_oc()
+    plt.tight_layout()
+    plt.savefig("decay_oc.png")
+    plt.show()
+    plt.close()