diff --git a/HDF5_Test.py b/HDF5_Test.py new file mode 100644 index 00000000..4d108539 --- /dev/null +++ b/HDF5_Test.py @@ -0,0 +1,19 @@ +from orbitize import results +import matplotlib.pyplot as plt +hdf5_filename = "my_posterior_brightness.hdf5" + + +loaded_results = results.Results() # Create blank results object for loading +loaded_results.load_results(hdf5_filename) + +print(loaded_results.post) + +orbit_plot_fig = loaded_results.plot_orbits( + object_to_plot=1, + num_orbits_to_plot=100, + start_mjd=0, + sep_pa_end_year=1861 +) +# Save or show the figure +orbit_plot_fig.savefig("reloaded_short_bright_orbit_plot.png") +plt.show() \ No newline at end of file diff --git a/HDF5_shortpos.py b/HDF5_shortpos.py new file mode 100644 index 00000000..b78dca41 --- /dev/null +++ b/HDF5_shortpos.py @@ -0,0 +1,20 @@ +from orbitize import results +import matplotlib.pyplot as plt +hdf5_filename = "/home/fmolina/orbitize/tests/my_posterior_shortpos.hdf5" + + +loaded_results = results.Results() # Create blank results object for loading +loaded_results.load_results(hdf5_filename) + +print(loaded_results.labels) + + +orbit_plot_fig = loaded_results.plot_orbits( + object_to_plot=1, + num_orbits_to_plot=100, + start_mjd=0, + sep_pa_end_year=1859 +) +# save fig +orbit_plot_fig.savefig("/home/fmolina/orbitize/reloaded_shortpos_plot.png") +plt.show() \ No newline at end of file diff --git a/mcmc_test.py b/mcmc_test.py new file mode 100644 index 00000000..79c096dc --- /dev/null +++ b/mcmc_test.py @@ -0,0 +1,82 @@ +import numpy as np + +import orbitize +from orbitize import driver +import multiprocessing as mp + +filename = "{}/simulated_ra_dec_data.csv".format(orbitize.DATADIR) + +# system parameters +num_secondary_bodies = 1 +total_mass = 1.75 # [Msol] +plx = 30 # [mas] +mass_err = 0.05 # [Msol] +plx_err = 0.12 # [mas] + +# MCMC parameters +num_temps = 20 +num_walkers = 1000 +num_threads = 20 # or a different number if you prefer, mp.cpu_count() for example + + +my_driver = driver.Driver( + filename, + "MCMC", + num_secondary_bodies, + total_mass, + plx, + mass_err=mass_err, + plx_err=plx_err, + mcmc_kwargs={ + "num_temps": num_temps, + "num_walkers": num_walkers, + "num_threads": num_threads, + }, system_kwargs={"restrict_angle_ranges": True}, +) + + +if __name__ == '__main__': + + total_orbits = 50000000 # number of steps x number of walkers (at lowest temperature) + burn_steps = 500000 # steps to burn in per walker + thin = 10 # only save every 2nd step + + my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + + corner_plot_fig = ( + my_driver.sampler.results.plot_corner() + ) # Creates a corner plot and returns Figure object + corner_plot_fig.savefig( + "my_corner_plot.png" + ) # This is matplotlib.figure.Figure.savefig() + + epochs = my_driver.system.data_table["epoch"] + + orbit_plot_fig = my_driver.sampler.results.plot_orbits( + object_to_plot=1, # Plot orbits for the first (and only, in this case) companion + num_orbits_to_plot=100, # Will plot 100 randomly selected orbits of this companion + start_mjd=epochs[0], # Minimum MJD for colorbar (here we choose first data epoch) + sep_pa_end_year=1870.0 + ) + orbit_plot_fig.savefig( + "my_orbit_plot.png" + ) # This is matplotlib.figure.Figure.savefig() + + sma_chains, ecc_chains = my_driver.sampler.examine_chains( + param_list=["sma1", "ecc1"], n_walkers=5 + ) + + from orbitize import results + + + hdf5_filename = "my_posterior.hdf5" + import os + +# To avoid weird behaviours, delete saved file if it already exists from a previous run of this notebook + if os.path.isfile(hdf5_filename): + os.remove(hdf5_filename) + + my_driver.sampler.results.save_results(hdf5_filename) + loaded_results = results.Results() # Create blank results object for loading + loaded_results.load_results(hdf5_filename) + diff --git a/orbitize/example_data/orbital_data_with_id.csv b/orbitize/example_data/orbital_data_with_id.csv new file mode 100644 index 00000000..c8f15a9d --- /dev/null +++ b/orbitize/example_data/orbital_data_with_id.csv @@ -0,0 +1,51 @@ +epoch,object,raoff,decoff,raoff_err,decoff_err,brightness,brightness_err +0.0,1,0.811231477773281,0.7594337387611091,0.01,0.01,0.09358490399640425,0.01 +51.02040816326531,1,0.8086976534668815,0.7884882174763149,0.01,0.01,0.09862733963920862,0.01 +102.04081632653062,1,0.8018051182543783,0.8135752998448389,0.01,0.01,0.10499582237388896,0.01 +153.06122448979593,1,0.7904122778426745,0.8345316311433326,0.01,0.01,0.11302378885599486,0.01 +204.08163265306123,1,0.7743165421613408,0.8511262472797302,0.01,0.01,0.12318805546636806,0.01 +255.10204081632654,1,0.753249548414535,0.8630503873758559,0.01,0.01,0.1361800248346027,0.01 +306.12244897959187,1,0.7268712100315663,0.8699046823746542,0.01,0.01,0.1530226526335016,0.01 +357.14285714285717,1,0.694762808612201,0.871183254587578,0.01,0.01,0.17526767302647034,0.01 +408.16326530612247,1,0.6564196199390899,0.8662543021471459,0.01,0.01,0.2053375622554458,0.01 +459.18367346938777,1,0.6112440678303533,0.8543369207025031,0.01,0.01,0.2471342724034155,0.01 +510.2040816326531,1,0.5585412933875817,0.8344743818340995,0.01,0.01,0.30714302736709537,0.01 +561.2244897959184,1,0.49752057854752274,0.8055051180185213,0.01,0.01,0.39642605066463393,0.01 +612.2448979591837,1,0.42730869340142713,0.7660347500990113,0.01,0.01,0.5339941531274536,0.01 +663.265306122449,1,0.3469855704957162,0.7144164850051454,0.01,0.01,0.7510921949147653,0.01 +714.2857142857143,1,0.25565954973242283,0.6487544841765712,0.01,0.01,1.089659065489047,0.01 +765.3061224489796,1,0.15260947588849744,0.5669573266103906,0.01,0.01,1.5668452007999316,0.01 +816.3265306122449,1,0.037533685832437026,0.4668885487044668,0.01,0.01,2.062626571880702,0.01 +867.3469387755102,1,-0.08904312701616794,0.34668838046440426,0.01,0.01,2.2566005333932724,0.01 +918.3673469387755,1,-0.22515974344986556,0.2053664293775624,0.01,0.01,1.9742762557365308,0.01 +969.3877551020408,1,-0.36661560366594714,0.04375763921362707,0.01,0.01,1.4724201265196528,0.01 +1020.4081632653061,1,-0.5061474469043882,-0.13417842573759436,0.01,0.01,1.0378267478370446,0.01 +1071.4285714285716,1,-0.6330566925723007,-0.32005017911521655,0.01,0.01,0.7457820123781923,0.01 +1122.4489795918369,1,-0.7340642827698997,-0.5005631649102175,0.01,0.01,0.5692403465256193,0.01 +1173.4693877551022,1,-0.7960233988208258,-0.6590354524244593,0.01,0.01,0.3677023181237994,0.01 +1224.4897959183675,1,-0.8100093753970008,-0.779341500288983,0.01,0.01,0.310805460038203,0.01 +1275.5102040816328,1,-0.7747106962319243,-0.8508197707379554,0.01,0.01,0.41184810403842464,0.01 +1326.530612244898,1,-0.6968201319288353,-0.8712731010397949,0.01,0.01,0.4267037403343691,0.01 +1377.5510204081634,1,-0.5881886796558359,-0.8463291454704698,0.01,0.01,0.45692118233449847,0.01 +1428.5714285714287,1,-0.46179505907320606,-0.7861524685145167,0.01,0.01,0.4835731029852539,0.01 +1479.591836734694,1,-0.32874319217715686,-0.7018568896193422,0.01,0.01,0.47696040679432183,0.01 +1530.6122448979593,1,-0.19704231837725736,-0.603169261080095,0.01,0.01,0.4173135946760028,0.01 +1581.6326530612246,1,-0.0717085052125464,-0.49753210156043054,0.01,0.01,0.32665907139175554,0.01 +1632.6530612244899,1,0.04455725697773338,-0.39013739669169883,0.01,0.01,0.24572306795533427,0.01 +1683.6734693877552,1,0.15061964729048624,-0.2843546293265123,0.01,0.01,0.189864891641152,0.01 +1734.6938775510205,1,0.246302320136533,-0.18222582266793214,0.01,0.01,0.1541293228877271,0.01 +1785.7142857142858,1,0.3319636092307561,-0.08488826595156958,0.01,0.01,0.1309447698576423,0.01 +1836.734693877551,1,0.4082247422781462,0.007108706138521259,0.01,0.01,0.11532666162626576,0.01 +1887.7551020408164,1,0.47580638683140664,0.09358245909783232,0.01,0.01,0.10445257538616691,0.01 +1938.7755102040817,1,0.5354357374912007,0.17456889152858981,0.01,0.01,0.09672357063471888,0.01 +1989.795918367347,1,0.5877970769939762,0.25022587899384413,0.01,0.01,0.09119647427470219,0.01 +2040.8163265306123,1,0.6335080358663165,0.3207715550593558,0.01,0.01,0.08729085845318117,0.01 +2091.8367346938776,1,0.673110446052766,0.3864454566166559,0.01,0.01,0.08463812457943871,0.01 +2142.857142857143,1,0.707069071798058,0.44748437788656903,0.01,0.01,0.0830008016213394,0.01 +2193.877551020408,1,0.7357742511645264,0.5041075168486121,0.01,0.01,0.08222762558027762,0.01 +2244.8979591836737,1,0.75954615418438,0.5565073629841016,0.01,0.01,0.08222788911616931,0.01 +2295.918367346939,1,0.7786393591852832,0.6048440073048987,0.01,0.01,0.08295689503656585,0.01 +2346.9387755102043,1,0.7932470300017178,0.6492413535736495,0.01,0.01,0.0844083805697474,0.01 +2397.9591836734694,1,0.8035043081850315,0.6897842167972191,0.01,0.01,0.08661187268900492,0.01 +2448.979591836735,1,0.8094907181927744,0.7265156142343624,0.01,0.01,0.08963415615760459,0.01 +2500.0,1,0.811231477773281,0.7594337387611084,0.01,0.01,0.09358490399640412,0.01 diff --git a/orbitize/example_data/reflected_light_example.csv b/orbitize/example_data/reflected_light_example.csv new file mode 100644 index 00000000..1d7cd32a --- /dev/null +++ b/orbitize/example_data/reflected_light_example.csv @@ -0,0 +1,10 @@ +epoch,object,raoff,decoff,raoff_err,decoff_err,brightness,brightness_err +57298,1,253.72,92.35,2.98,2.85,0.5,0.1 +57606,1,236.63,127.94,9.77,9.18,0.5,0.1 +57645,1,234.52,123.39,1.79,1.03,0.5,0.1 +57946,1,210.76,152.09,1.94,1.88,0.5,0.1 +58276,1,167.49,180.87,1.61,16.97,, +58287,1,177.67,174.6,1.67,1.67,0.5,0.1 +58365,1,165.7,185.33,3.28,3.66,0.5,0.1 +58368,1,170.38,185.94,2.52,2.74,0.5,0.1 +58414,1,161.64,176.21,13.6,14.31,, \ No newline at end of file diff --git a/orbitize/example_data/simulated_ra_dec_data.csv b/orbitize/example_data/simulated_ra_dec_data.csv new file mode 100644 index 00000000..afc737cf --- /dev/null +++ b/orbitize/example_data/simulated_ra_dec_data.csv @@ -0,0 +1,51 @@ +epoch,object,raoff,decoff,raoff_err,decoff_err +0.0,1,-44.97595264191645,-40.40063509461096,0.01,0.01 +0.02040816326530612,1,-44.30820072074837,-45.71952453300766,0.01,0.01 +0.04081632653061224,1,-41.04700178820174,-48.201058902486146,0.01,0.01 +0.061224489795918366,1,-35.74054958085938,-47.965832008749146,0.01,0.01 +0.08163265306122448,1,-29.106223608047067,-45.48435880817074,0.01,0.01 +0.1020408163265306,1,-21.8253835127802,-41.370259974357374,0.01,0.01 +0.12244897959183673,1,-14.428317330209142,-36.211598246977516,0.01,0.01 +0.14285714285714285,1,-7.2699135408338,-30.48491432830629,0.01,0.01 +0.16326530612244897,1,-0.5565700744014579,-24.536472548444245,0.01,0.01 +0.18367346938775508,1,5.611711944213546,-18.598583371815057,0.01,0.01 +0.2040816326530612,1,11.203113265151215,-12.816822522196352,0.01,0.01 +0.22448979591836732,1,16.22563205150931,-7.27593756563277,0.01,0.01 +0.24489795918367346,1,20.708218518339898,-2.020291901969007,0.01,0.01 +0.26530612244897955,1,24.68910142063662,2.93143637707865,0.01,0.01 +0.2857142857142857,1,28.208969483025083,7.576234944513116,0.01,0.01 +0.3061224489795918,1,31.307217197702858,11.920150575624033,0.01,0.01 +0.32653061224489793,1,34.02003706808991,15.974029767980111,0.01,0.01 +0.3469387755102041,1,36.37958133646282,19.75081855840348,0.01,0.01 +0.36734693877551017,1,38.41371660681627,23.263874472541467,0.01,0.01 +0.3877551020408163,1,40.14608684277597,26.525923557689968,0.01,0.01 +0.4081632653061224,1,41.59631855303417,29.548420660047412,0.01,0.01 +0.42857142857142855,1,42.7802730988534,32.34115487431027,0.01,0.01 +0.44897959183673464,1,43.71029297512206,34.91199683952046,0.01,0.01 +0.4693877551020408,1,44.3954130505847,37.266719491329106,0.01,0.01 +0.4897959183673469,1,44.8415215664951,39.408846678587985,0.01,0.01 +0.5102040816326531,1,45.05146258834253,41.339496315275746,0.01,0.01 +0.5306122448979591,1,45.02507581579308,43.057194255088284,0.01,0.01 +0.5510204081632653,1,44.75917087078464,44.55763784084132,0.01,0.01 +0.5714285714285714,1,44.247434007528476,45.83339021912297,0.01,0.01 +0.5918367346938775,1,43.480265674963206,46.87348554962607,0.01,0.01 +0.6122448979591836,1,42.444548173623666,47.66292429386228,0.01,0.01 +0.6326530612244897,1,41.12334547730126,48.182034574487595,0.01,0.01 +0.6530612244897959,1,39.495542722455234,48.40567440725761,0.01,0.01 +0.673469387755102,1,37.53544417786146,48.302249377305294,0.01,0.01 +0.6938775510204082,1,35.212369445793755,47.832525882942925,0.01,0.01 +0.7142857142857142,1,32.49032544307804,46.948236929022904,0.01,0.01 +0.7346938775510203,1,29.32789804090339,45.5905171549874,0.01,0.01 +0.7551020408163265,1,25.678621015774937,43.68828706664445,0.01,0.01 +0.7755102040816326,1,21.49227022844291,41.156869864693746,0.01,0.01 +0.7959183673469387,1,16.717838261070998,37.89742910959155,0.01,0.01 +0.8163265306122448,1,11.309414286826271,33.79835583006353,0.01,0.01 +0.836734693877551,1,5.236842173178732,28.74062952787787,0.01,0.01 +0.8571428571428571,1,-1.4962572079392693,22.610509345316203,0.01,0.01 +0.8775510204081632,1,-8.824224576507511,15.324492907047599,0.01,0.01 +0.8979591836734693,1,-16.579434040144825,6.872290977814089,0.01,0.01 +0.9183673469387754,1,-24.445019923631108,-2.619219253000618,0.01,0.01 +0.9387755102040816,1,-31.917621523654393,-12.809602893650707,0.01,0.01 +0.9591836734693877,1,-38.3162586778345,-23.088729032293035,0.01,0.01 +0.9795918367346939,1,-42.88164071438466,-32.600661274349676,0.01,0.01 +1.0,1,-44.97595264191646,-40.40063509461095,0.01,0.01 diff --git a/orbitize/kepler.py b/orbitize/kepler.py index 447bbcc7..60e4d5af 100644 --- a/orbitize/kepler.py +++ b/orbitize/kepler.py @@ -1,6 +1,7 @@ """ This module solves for the orbit of the planet given Keplerian parameters. """ + import numpy as np import astropy.units as u import astropy.constants as consts @@ -13,30 +14,29 @@ if cuda_ext: # Configure GPU context for CUDA accelerated compute from orbitize import gpu_context + kep_gpu_ctx = gpu_context.gpu_context() + def tau_to_manom(date, sma, mtot, tau, tau_ref_epoch): """ Gets the mean anomlay - + Args: date (float or np.array): MJD sma (float): semi major axis (AU) mtot (float): total mass (M_sun) tau (float): epoch of periastron, in units of the orbital period tau_ref_epoch (float): reference epoch for tau - + Returns: float or np.array: mean anomaly on that date [0, 2pi) """ - period = np.sqrt( - 4 * np.pi**2.0 * (sma * u.AU)**3 / - (consts.G * (mtot * u.Msun)) - ) + period = np.sqrt(4 * np.pi**2.0 * (sma * u.AU) ** 3 / (consts.G * (mtot * u.Msun))) period = period.to(u.day).value - frac_date = (date - tau_ref_epoch)/period + frac_date = (date - tau_ref_epoch) / period frac_date %= 1 mean_anom = (frac_date - tau) * 2 * np.pi @@ -45,11 +45,65 @@ def tau_to_manom(date, sma, mtot, tau, tau_ref_epoch): return mean_anom -def calc_orbit( - epochs, sma, ecc, inc, aop, pan, tau, plx, mtot, mass_for_Kamp=None, tau_ref_epoch=58849, tolerance=1e-9, - max_iter=100, use_c=True, use_gpu=False +def times2trueanom_and_eccanom( + sma, + epochs, + mtot, + ecc, + tau, + tau_ref_epoch=58849, + tolerance=1e-9, + max_iter=100, + use_c=True, + use_gpu=False, ): + + + n_orbs = np.size(sma) # num sets of input orbital parameters + n_dates = np.size(epochs) # number of dates to compute offsets and vz + + # Necessary for _calc_ecc_anom, for now + if np.isscalar(epochs): # just in case epochs is given as a scalar + epochs = np.array([epochs]) + ecc_arr = np.tile(ecc, (n_dates, 1)) + + # # compute mean anomaly (size: n_orbs x n_dates) + manom = tau_to_manom(epochs[:, None], sma, mtot, tau, tau_ref_epoch) + # compute eccentric anomalies (size: n_orbs x n_dates) + eanom = _calc_ecc_anom( + manom, + ecc_arr, + tolerance=tolerance, + max_iter=max_iter, + use_c=use_c, + use_gpu=use_gpu, + ) + + # compute the true anomalies (size: n_orbs x n_dates) + # Note: matrix multiplication makes the shapes work out here and below + tanom = 2.0 * np.arctan(np.sqrt((1.0 + ecc) / (1.0 - ecc)) * np.tan(0.5 * eanom)) + + return tanom, eanom + + +def calc_orbit( + epochs, + sma, + ecc, + inc, + aop, + pan, + tau, + plx, + mtot, + mass_for_Kamp=None, + tau_ref_epoch=58849, + tolerance=1e-9, + max_iter=100, + use_c=True, + use_gpu=False, +): """ Returns the separation and radial velocity of the body given array of orbital parameters (size n_orbs) at given epochs (array of size n_dates) @@ -89,33 +143,32 @@ def calc_orbit( Written: Jason Wang, Henry Ngo, 2018 """ - n_orbs = np.size(sma) # num sets of input orbital parameters - n_dates = np.size(epochs) # number of dates to compute offsets and vz # return planetary RV if `mass_for_Kamp` is not defined if mass_for_Kamp is None: mass_for_Kamp = mtot + ecc + + tanom, eanom = times2trueanom_and_eccanom( + sma, + epochs, + mtot, + ecc, + tau, + tau_ref_epoch=tau_ref_epoch, + tolerance=tolerance, + max_iter=max_iter, + use_c=use_c, + use_gpu=use_gpu, + ) - # Necessary for _calc_ecc_anom, for now - if np.isscalar(epochs): # just in case epochs is given as a scalar - epochs = np.array([epochs]) - ecc_arr = np.tile(ecc, (n_dates, 1)) - - # # compute mean anomaly (size: n_orbs x n_dates) - manom = tau_to_manom(epochs[:, None], sma, mtot, tau, tau_ref_epoch) - # compute eccentric anomalies (size: n_orbs x n_dates) - eanom = _calc_ecc_anom(manom, ecc_arr, tolerance=tolerance, max_iter=max_iter, use_c=use_c, use_gpu=use_gpu) - - # compute the true anomalies (size: n_orbs x n_dates) - # Note: matrix multiplication makes the shapes work out here and below - tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom)) # compute 3-D orbital radius of second body (size: n_orbs x n_dates) radius = sma * (1.0 - ecc * np.cos(eanom)) # compute ra/dec offsets (size: n_orbs x n_dates) # math from James Graham. Lots of trig - c2i2 = np.cos(0.5*inc)**2 - s2i2 = np.sin(0.5*inc)**2 + c2i2 = np.cos(0.5 * inc) ** 2 + s2i2 = np.sin(0.5 * inc) ** 2 arg1 = tanom + aop + pan arg2 = tanom + aop - pan c1 = np.cos(arg1) @@ -124,39 +177,46 @@ def calc_orbit( s2 = np.sin(arg2) # updated sign convention for Green Eq. 19.4-19.7 - raoff = radius * (c2i2*s1 - s2i2*s2) * plx - deoff = radius * (c2i2*c1 + s2i2*c2) * plx + raoff = radius * (c2i2 * s1 - s2i2 * s2) * plx + deoff = radius * (c2i2 * c1 + s2i2 * c2) * plx # compute the radial velocity (vz) of the body (size: n_orbs x n_dates) # first comptue the RV semi-amplitude (size: n_orbs x n_dates) - Kv = np.sqrt(consts.G / (1.0 - ecc**2)) * (mass_for_Kamp * u.Msun * - np.sin(inc)) / np.sqrt(mtot * u.Msun) / np.sqrt(sma * u.au) + Kv = ( + np.sqrt(consts.G / (1.0 - ecc**2)) + * (mass_for_Kamp * u.Msun * np.sin(inc)) + / np.sqrt(mtot * u.Msun) + / np.sqrt(sma * u.au) + ) # Convert to km/s - Kv = Kv.to(u.km/u.s) + Kv = Kv.to(u.km / u.s) # compute the vz - vz = Kv.value * (ecc*np.cos(aop) + np.cos(aop + tanom)) + vz = Kv.value * (ecc * np.cos(aop) + np.cos(aop + tanom)) # Squeeze out extra dimension (useful if n_orbs = 1, does nothing if n_orbs > 1) vz = np.squeeze(vz)[()] return raoff, deoff, vz -def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gpu=False): + +def _calc_ecc_anom( + manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gpu=False +): """ - Computes the eccentric anomaly from the mean anomlay. - Code from Rob De Rosa's orbit solver (e < 0.95 use Newton, e >= 0.95 use Mikkola) + Computes the eccentric anomaly from the mean anomlay. + Code from Rob De Rosa's orbit solver (e < 0.95 use Newton, e >= 0.95 use Mikkola) - Args: - manom (float/np.array): mean anomaly, either a scalar or np.array of any shape - ecc (float/np.array): eccentricity, either a scalar or np.array of the same shape as manom - tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. - max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. - use_c (bool, optional): Use the C solver if configured. Defaults to False - use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False + Args: + manom (float/np.array): mean anomaly, either a scalar or np.array of any shape + ecc (float/np.array): eccentricity, either a scalar or np.array of the same shape as manom + tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. + max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. + use_c (bool, optional): Use the C solver if configured. Defaults to False + use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False -Return: - eanom (float/np.array): eccentric anomalies, same shape as manom + Return: + eanom (float/np.array): eccentric anomalies, same shape as manom - Written: Jason Wang, 2018 + Written: Jason Wang, 2018 """ if np.isscalar(ecc) or (np.shape(manom) == np.shape(ecc)): @@ -166,7 +226,7 @@ def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gp # If manom is a scalar, make it into a one-element array if np.isscalar(manom): - manom = np.array((manom, )) + manom = np.array((manom,)) # If ecc is a scalar, make it the same shape as manom if np.isscalar(ecc): @@ -186,19 +246,26 @@ def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gp # Now low eccentricities ind_low = np.where(~ecc_zero & ecc_low) - if len(ind_low[0]) > 0: - eanom[ind_low] = _newton_solver_wrapper(manom[ind_low], ecc[ind_low], tolerance, max_iter, use_c, use_gpu) - + if len(ind_low[0]) > 0: + eanom[ind_low] = _newton_solver_wrapper( + manom[ind_low], ecc[ind_low], tolerance, max_iter, use_c, use_gpu + ) + # Now high eccentricities - ind_high = np.where(~ecc_zero & ~ecc_low | (eanom == -1)) # The C and CUDA solvers return the unphysical value -1 if they fail to converge - if len(ind_high[0]) > 0: - eanom[ind_high] = _mikkola_solver_wrapper(manom[ind_high], ecc[ind_high], use_c, use_gpu) + ind_high = np.where( + ~ecc_zero & ~ecc_low | (eanom == -1) + ) # The C and CUDA solvers return the unphysical value -1 if they fail to converge + if len(ind_high[0]) > 0: + eanom[ind_high] = _mikkola_solver_wrapper( + manom[ind_high], ecc[ind_high], use_c, use_gpu + ) return np.squeeze(eanom)[()] + def _newton_solver_wrapper(manom, ecc, tolerance, max_iter, use_c=False, use_gpu=False): """ - Wrapper for the various (Python, C, CUDA) implementations of the Newton-Raphson solver + Wrapper for the various (Python, C, CUDA) implementations of the Newton-Raphson solver for eccentric anomaly. Args: @@ -213,18 +280,21 @@ def _newton_solver_wrapper(manom, ecc, tolerance, max_iter, use_c=False, use_gpu Written: Devin Cody, 2021 """ eanom = np.empty_like(manom) - + if cuda_ext and use_gpu: # the CUDA solver returns eanom = -1 if it doesnt converge after max_iter iterations eanom = _CUDA_newton_solver(manom, ecc, tolerance=tolerance, max_iter=max_iter) elif cext and use_c: # the C solver returns eanom = -1 if it doesnt converge after max_iter iterations - eanom = _kepler._c_newton_solver(manom, ecc, tolerance=tolerance, max_iter=max_iter) + eanom = _kepler._c_newton_solver( + manom, ecc, tolerance=tolerance, max_iter=max_iter + ) else: eanom = _newton_solver(manom, ecc, tolerance=tolerance, max_iter=max_iter) return eanom + def _newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): """ Newton-Raphson solver for eccentric anomaly. @@ -232,11 +302,11 @@ def _newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): Args: manom (np.array): array of mean anomalies ecc (np.array): array of eccentricities - tolerance (float, optional): absolute tolerance of iterative computation. + tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. - max_iter (int, optional): maximum number of iterations before switching. + max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. - eanom0 (np.array): array of first guess for eccentric anomaly, same + eanom0 (np.array): array of first guess for eccentric anomaly, same shape as manom (optional) @@ -262,23 +332,29 @@ def _newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): abs_diff = np.abs(diff) ind = np.where(abs_diff > tolerance) niter = 0 - while ((ind[0].size > 0) and (niter <= max_iter)): + while (ind[0].size > 0) and (niter <= max_iter): eanom[ind] -= diff[ind] # If it hasn't converged after half the iterations are done, try starting from pi - if niter == (max_iter//2): + if niter == (max_iter // 2): eanom[ind] = np.pi - diff[ind] = (eanom[ind] - (ecc[ind] * np.sin(eanom[ind])) - manom[ind]) / \ - (1.0 - (ecc[ind] * np.cos(eanom[ind]))) + diff[ind] = (eanom[ind] - (ecc[ind] * np.sin(eanom[ind])) - manom[ind]) / ( + 1.0 - (ecc[ind] * np.cos(eanom[ind])) + ) abs_diff[ind] = np.abs(diff[ind]) ind = np.where(abs_diff > tolerance) niter += 1 if niter >= max_iter: - print(manom[ind], eanom[ind], diff[ind], ecc[ind], '> {} iter.'.format(max_iter)) - eanom[ind] = _mikkola_solver_wrapper(manom[ind], ecc[ind]) # Send remaining orbits to the analytical version, this has not happened yet... + print( + manom[ind], eanom[ind], diff[ind], ecc[ind], "> {} iter.".format(max_iter) + ) + eanom[ind] = _mikkola_solver_wrapper( + manom[ind], ecc[ind] + ) # Send remaining orbits to the analytical version, this has not happened yet... return eanom + def _CUDA_newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): """ Helper function for calling the CUDA implementation of the Newton-Raphson solver for eccentric anomaly. @@ -298,16 +374,17 @@ def _CUDA_newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): manom = np.asarray(manom) ecc = np.asarray(ecc) eanom = np.empty_like(manom) - tolerance = np.asarray(tolerance, dtype = np.float64) + tolerance = np.asarray(tolerance, dtype=np.float64) max_iter = np.asarray(max_iter) - + kep_gpu_ctx.newton(manom, ecc, eanom, eanom0, tolerance, max_iter) return eanom + def _mikkola_solver_wrapper(manom, ecc, use_c=False, use_gpu=False): """ - Wrapper for the various (Python, C, CUDA) implementations of Analtyical Mikkola solver + Wrapper for the various (Python, C, CUDA) implementations of Analtyical Mikkola solver Args: manom (np.array): array of mean anomalies between 0 and 2pi @@ -353,26 +430,29 @@ def _mikkola_solver(manom, ecc): beta = (0.5 * manom) / ((4.0 * ecc) + 0.5) aux = np.sqrt(beta**2.0 + alpha**3.0) - z = np.abs(beta + aux)**(1.0/3.0) + z = np.abs(beta + aux) ** (1.0 / 3.0) - s0 = z - (alpha/z) - s1 = s0 - (0.078*(s0**5.0)) / (1.0 + ecc) - e0 = manom + (ecc * (3.0*s1 - 4.0*(s1**3.0))) + s0 = z - (alpha / z) + s1 = s0 - (0.078 * (s0**5.0)) / (1.0 + ecc) + e0 = manom + (ecc * (3.0 * s1 - 4.0 * (s1**3.0))) se0 = np.sin(e0) ce0 = np.cos(e0) - f = e0-ecc*se0-manom - f1 = 1.0-ecc*ce0 - f2 = ecc*se0 - f3 = ecc*ce0 + f = e0 - ecc * se0 - manom + f1 = 1.0 - ecc * ce0 + f2 = ecc * se0 + f3 = ecc * ce0 f4 = -f2 - u1 = -f/f1 - u2 = -f/(f1+0.5*f2*u1) - u3 = -f/(f1+0.5*f2*u2+(1.0/6.0)*f3*u2*u2) - u4 = -f/(f1+0.5*f2*u3+(1.0/6.0)*f3*u3*u3+(1.0/24.0)*f4*(u3**3.0)) + u1 = -f / f1 + u2 = -f / (f1 + 0.5 * f2 * u1) + u3 = -f / (f1 + 0.5 * f2 * u2 + (1.0 / 6.0) * f3 * u2 * u2) + u4 = -f / ( + f1 + 0.5 * f2 * u3 + (1.0 / 6.0) * f3 * u3 * u3 + (1.0 / 24.0) * f4 * (u3**3.0) + ) + + return e0 + u4 - return (e0 + u4) def _CUDA_mikkola_solver(manom, ecc): """ diff --git a/orbitize/plot.py b/orbitize/plot.py index a02cfc1e..1177115f 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -173,7 +173,7 @@ def plot_orbits( show_colorbar=True, cmap=cmap, sep_pa_color="lightgrey", - sep_pa_end_year=2025.0, + sep_pa_end_year=1875.0, cbar_param="Epoch [year]", mod180=False, rv_time_series=False, @@ -1077,7 +1077,7 @@ def plot_residuals( num_orbits_to_plot=100, num_epochs_to_plot=100, sep_pa_color="lightgrey", - sep_pa_end_year=2025.0, + sep_pa_end_year=1925.0, cbar_param="Epoch [year]", mod180=False, ): diff --git a/orbitize/read_input.py b/orbitize/read_input.py index f37f5e41..cd6f6a9a 100644 --- a/orbitize/read_input.py +++ b/orbitize/read_input.py @@ -179,6 +179,10 @@ def read_file(filename): have_seppacorr = np.zeros( num_measurements, dtype=bool ) # zeros are False + if "brightness" in input_table.columns: + have_brightness = ~input_table["brightness"].mask + else: + have_brightness = np.zeros(num_measurements, dtype=bool) if "rv" in input_table.columns: have_rv = ~input_table["rv"].mask else: @@ -231,11 +235,14 @@ def read_file(filename): else: have_rv = np.zeros(num_measurements, dtype=bool) # zeros are False - # Rob: not sure if we need this but adding just in case if "instrument" in input_table.columns: have_inst = np.ones(num_measurements, dtype=bool) else: have_inst = np.zeros(num_measurements, dtype=bool) + if "brightness" in input_table.columns: + have_brightness = np.ones(num_measurements, dtype=bool) + else: + have_brightness = np.zeros(num_measurements, dtype=bool) # orbitize! backwards compatability since we added new columns, some old data formats may not have them # fill in with default values @@ -280,9 +287,9 @@ def read_file(filename): if not isinstance(row["object"], (int, np.int32, np.int64)): raise Exception("Invalid object ID. Object IDs must be integers.") - # determine input quantity type (RA/DEC, SEP/PA, or RV) + # determine input quantity type (RA/DEC, SEP/PA, RV, or BRIGHTNESS) if orbitize_style: - if row["quant_type"] == "rv": # special format for rv rows + if row["quant_type"] == "rv" or row["quant_type"] == 'brightness': # special format for rv rows output_table.add_row( [ MJD, @@ -332,6 +339,7 @@ def read_file(filename): ) else: # When not in orbitize style + if have_ra[index] and have_dec[index]: # check if there's a covariance term if have_radeccorr[index]: @@ -428,6 +436,20 @@ def read_file(filename): "defrv", ] ) + if have_brightness[index]: + output_table.add_row( + [ + MJD, + row["object"], + row["brightness"], + row["brightness_err"], + None, + None, + None, + "brightness", + "defbr", + ] + ) return output_table diff --git a/orbitize/system.py b/orbitize/system.py index 1da30d3f..55e8a97a 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -98,10 +98,15 @@ def __init__( # List of index arrays corresponding to each rv for each body self.rv = [] + # index arrays corresponding to brightness for each body + self.brightness = [] + self.fit_astrometry = True radec_indices = np.where(self.data_table["quant_type"] == "radec") seppa_indices = np.where(self.data_table["quant_type"] == "seppa") + brightness_indices = np.where(self.data_table["quant_type"] == "brightness") + if len(radec_indices[0]) == 0 and len(seppa_indices[0]) == 0: self.fit_astrometry = False rv_indices = np.where(self.data_table["quant_type"] == "rv") @@ -140,6 +145,7 @@ def __init__( np.intersect1d(self.body_indices[body_num], seppa_indices) ) self.rv.append(np.intersect1d(self.body_indices[body_num], rv_indices)) + self.brightness.append(np.intersect1d(self.body_indices[body_num], brightness_indices)) # we should track the influence of the planet(s) on each other/the star if: # we are not fitting massless planets and @@ -302,6 +308,7 @@ def __init__( self.param_idx = self.basis.param_idx + def save(self, hf): """ Saves the current object to an hdf5 file @@ -380,6 +387,7 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): dec_perturb = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits)) vz = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits)) + brightness_out = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits)) # mass/mtot used to compute each Keplerian orbit will be needed later to compute perturbations if self.track_planet_perturbs: @@ -483,16 +491,33 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): tau_ref_epoch=self.tau_ref_epoch, ) + tanom, eanom = kepler.times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=self.tau_ref_epoch) + + + R = (sma*(1-ecc**2))/(1+ecc*np.cos(tanom)) + + z = (R)*(-np.cos(argp)*np.sin(inc)*np.sin(tanom)-np.cos(tanom)*np.sin(inc)*np.sin(argp)) + + B = np.arctan2(-R, z)+ np.pi + + Alpha = (1/np.pi)*(np.sin(B)+(np.pi-B)*np.cos(B)) + + Albedo = 0.5 + brightness = Albedo*Alpha/R**2 + + # raoff, decoff, vz are scalers if the length of epochs is 1 if len(epochs) == 1: raoff = np.array([raoff]) decoff = np.array([decoff]) vz_i = np.array([vz_i]) + brightness_out = np.array([brightness]) # add Keplerian ra/deoff for this body to storage arrays ra_kepler[:, body_num, :] = np.reshape(raoff, (n_epochs, n_orbits)) dec_kepler[:, body_num, :] = np.reshape(decoff, (n_epochs, n_orbits)) vz[:, body_num, :] = np.reshape(vz_i, (n_epochs, n_orbits)) + brightness_out[:, body_num, :] = np.reshape(brightness, (n_epochs, n_orbits)) # vz_i is the ith companion radial velocity if self.fit_secondary_mass: @@ -564,7 +589,7 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): else: return raoff, deoff, vz else: - return raoff, deoff, vz + return raoff, deoff, vz, brightness_out def compute_model(self, params_arr, use_rebound=False): """ @@ -599,7 +624,7 @@ def compute_model(self, params_arr, use_rebound=False): standard_params_arr, comp_rebound=True ) else: - raoff, decoff, vz = self.compute_all_orbits(standard_params_arr) + raoff, decoff, vz, brightness = self.compute_all_orbits(standard_params_arr) if len(standard_params_arr.shape) == 1: n_orbits = 1 @@ -651,6 +676,14 @@ def compute_model(self, params_arr, use_rebound=False): model[self.rv[body_num], 0] = vz[self.rv[body_num], body_num, :] model[self.rv[body_num], 1] = np.nan + # TODO (farrah): add brightness to model here (use RV block above as template) + # Brightness + if len(self.brightness[body_num]) > 0: + model[self.brightness[body_num], 0] = brightness[self.brightness[body_num], body_num, :] + model[self.brightness[body_num], 1] = np.nan + # (assume self.brightness is array of indices of epochs with brightness measurements) + + # if we have abs astrometry measurements in the input file (i.e. not # from Hipparcos or Gaia), add the parallactic & proper motion here by # calling AbsAstrom compute_model diff --git a/tests/Visual.py b/tests/Visual.py new file mode 100644 index 00000000..6e418eb9 --- /dev/null +++ b/tests/Visual.py @@ -0,0 +1,50 @@ +from orbitize import system, read_input, DATADIR +import os +import numpy as np +import matplotlib.pyplot as plt + + +num_secondary_bodies = 1 + +input_file = os.path.join(DATADIR, "GJ504.csv") +data_table = read_input.read_file(input_file) + +system_mass = 1.47 +plx = 24.30 + +test_system = system.System(num_secondary_bodies, data_table, system_mass, plx) + +params_arr = np.array( + [ + 10.0, # sma + 0.9, # ecc + np.radians(90), # inc + np.radians(60), # aop + np.radians(90), # pan + 0.0, # tau + plx, + system_mass, + ] +) +epochs = np.linspace(0, 365 * 30, int(1e3)) +ra, dec, vz, brightness = test_system.compute_all_orbits(params_arr, epochs=epochs) + +fig, ax = plt.subplots(2, 1, figsize=(5, 10)) + +ax[0].scatter( + epochs, + brightness, + color=plt.cm.RdYlBu((epochs - epochs[0]) / (epochs[-1] - epochs[0])), +) + +ax[1].scatter( + ra[:, 1, :], + dec[:, 1, :], + color=plt.cm.RdYlBu((epochs - epochs[0]) / (epochs[-1] - epochs[0])), +) + +ax[1].scatter([0], [0], color="red") + +ax[1].axis("equal") + +plt.savefig("visual4farrah.png") \ No newline at end of file diff --git a/tests/brightness_mcmc_test.py b/tests/brightness_mcmc_test.py new file mode 100644 index 00000000..1e4dbeb4 --- /dev/null +++ b/tests/brightness_mcmc_test.py @@ -0,0 +1,82 @@ +import numpy as np + +import orbitize +from orbitize import driver +import multiprocessing as mp + +filename = "{}/orbital_data_with_id.csv".format(orbitize.DATADIR) + +# system parameters +num_secondary_bodies = 1 +total_mass = 1.75 # [Msol] +plx = 51.44 # [mas] +mass_err = 0.05 # [Msol] +plx_err = 0.12 # [mas] + +# MCMC parameters +num_temps = 20 +num_walkers = 1000 +num_threads = 20 # or a different number if you prefer, mp.cpu_count() for example + + +my_driver = driver.Driver( + filename, + "MCMC", + num_secondary_bodies, + total_mass, + plx, + mass_err=mass_err, + plx_err=plx_err, + mcmc_kwargs={ + "num_temps": num_temps, + "num_walkers": num_walkers, + "num_threads": num_threads, + }, system_kwargs={"restrict_angle_ranges": True}, +) + + +if __name__ == '__main__': + + total_orbits = 100000000 # number of steps x number of walkers (at lowest temperature) + burn_steps = 50000 # steps to burn in per walker + thin = 10 # only save every 2nd step + + my_driver.sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + + corner_plot_fig = ( + my_driver.sampler.results.plot_corner() + ) # Creates a corner plot and returns Figure object + corner_plot_fig.savefig( + "my_corner_plot.png" + ) # This is matplotlib.figure.Figure.savefig() + + epochs = my_driver.system.data_table["epoch"] + + orbit_plot_fig = my_driver.sampler.results.plot_orbits( + object_to_plot=1, # Plot orbits for the first (and only, in this case) companion + num_orbits_to_plot=100, # Will plot 100 randomly selected orbits of this companion + start_mjd=epochs[0], # Minimum MJD for colorbar (here we choose first data epoch) + sep_pa_end_year=1870.0 + ) + orbit_plot_fig.savefig( + "my_orbit_plot.png" + ) # This is matplotlib.figure.Figure.savefig() + + sma_chains, ecc_chains = my_driver.sampler.examine_chains( + param_list=["sma1", "ecc1"], n_walkers=5 + ) + + from orbitize import results + + + hdf5_filename = "my_posterior.hdf5" + import os + +# To avoid weird behaviours, delete saved file if it already exists from a previous run of this notebook + if os.path.isfile(hdf5_filename): + os.remove(hdf5_filename) + + my_driver.sampler.results.save_results(hdf5_filename) + loaded_results = results.Results() # Create blank results object for loading + loaded_results.load_results(hdf5_filename) + diff --git a/tests/save_chain_corner_plot.py b/tests/save_chain_corner_plot.py new file mode 100644 index 00000000..9d02bd55 --- /dev/null +++ b/tests/save_chain_corner_plot.py @@ -0,0 +1,9 @@ +# Save Chains Corner Plot +import orbitize.results + +import matplotlib.pyplot as plt + +myResults = orbitize.results.Results() # create empty Results object +myResults.load_results("my_posterior.hdf5") # load from file + +corner_figure = myResults.plot_corner() \ No newline at end of file diff --git a/tests/test_brightness.py b/tests/test_brightness.py new file mode 100644 index 00000000..7a3f3d94 --- /dev/null +++ b/tests/test_brightness.py @@ -0,0 +1,131 @@ +""" +Tests for the reflected-light calculation in system.compute_all_orbits +""" + +from orbitize import system, sampler, plot +from orbitize import DATADIR, read_input +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from astropy.time import Time + + +def test_brightness_calculation(make_plot=False): + + num_secondary_bodies = 1 + + # input_file = os.path.join(DATADIR, "GJ504.csv") + input_file = os.path.join(DATADIR, "betaPic.csv") + data_table = read_input.read_file(input_file) + + times = data_table["epoch"].value + + system_mass = 1.47 + plx = 24.30 + + test_system = system.System(num_secondary_bodies, data_table, system_mass, plx) + + params = np.array( + [ + 10.0, + 0.3, + np.radians(89), + np.radians(21), + np.radians(70), + 0.0, + 51.5, + 1.75, + ] + ) + + ra, dec, vz, brightness = test_system.compute_all_orbits(params) + + if make_plot: + plt.figure() + plt.scatter(times, brightness) + plt.xlabel("Time [dy]", fontsize=18) + plt.ylabel("Brightness", fontsize=18) + plt.savefig("Test_brightness2.png") + + +def test_read_input_with_brightness(): + + num_secondary_bodies = 1 + + input_file = os.path.join(DATADIR, "reflected_light_example.csv") + + data_table = read_input.read_file(input_file) + + times = data_table["epoch"].value + + print(data_table) + + # TODO (Farrah): add a test that asserts the brightness column of the data table is + # what you expect (hint: check in the reflected_light_example.csv to see what + # the brightness values should be + + +def test_assert_nan(): + num_secondary_bodies = 1 + + input_file = os.path.join(DATADIR, "reflected_light_example.csv") + + data_table = read_input.read_file(input_file) + + brightness_values = data_table["brightness"].value + assert pd.isna(data_table["brightness"][4]), "This is supposed to say NAN" + print(brightness_values) + +def test_compute_posteriors(): + + num_secondary_bodies = 1 + + input_file = os.path.join(DATADIR, "orbital_data_with_id.csv") + data_table = read_input.read_file(input_file) + + system_mass = 1.47 + plx = 24.30 + + test_system = system.System(num_secondary_bodies, data_table, system_mass, plx) + + params_arr = np.array( + [ + 10.0, # sma + 0.3, # ecc + np.radians(0), # inc + np.radians(45), # aop + np.radians(90), # pan + 0.0, # tau + 51.5, # plx + 1.75, # stellar mass + ] + ) + epochs = np.linspace(0, 365*30, int(1e3)) + ra, dec, vz, brightness = test_system.compute_all_orbits(params_arr, epochs=epochs) + + # fig, ax = plt.subplots(2, 1, figsize=(5,10)) + + + # ax[0].scatter(epochs, brightness, color=plt.cm.RdYlBu((epochs-epochs[0])/(epochs[-1] - epochs[0]))) + # ax[1].scatter(ra[:,1,:], dec[:,1,:], color=plt.cm.RdYlBu((epochs-epochs[0])/(epochs[-1] - epochs[0]))) + + # ax[1].axis('equal') + # plt.savefig('visual4farrah.png') + + model = test_system.compute_model(params_arr) + + test_mcmc = sampler.MCMC(test_system, 1, 50) + + test_mcmc.run_sampler(10) + + # fig = plot.plot_orbits(test_mcmc.results, start_mjd=Time(2000., format='decimalyear').mjd, sep_pa_end_year=2500.) + + # plt.savefig('foo.png') + + + +if __name__ == "__main__": + #test_brightness_calculation() + # test_read_input_with_brightness() + test_compute_posteriors()