From 3069841f6f39ac2783b1429d370f5dcd007a6507 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 7 Jul 2026 13:47:25 -0700 Subject: [PATCH 1/8] remove unnecessary code --- orbitize/plot.py | 164 +++++++---------------------------------------- 1 file changed, 22 insertions(+), 142 deletions(-) diff --git a/orbitize/plot.py b/orbitize/plot.py index 26373519..1628e400 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -612,39 +612,10 @@ def plot_orbits( ) # Calculate ra/dec offsets for all epochs of this orbit - if (rv_time_series == True) or (rv_time_series2 == True): - raoff0, deoff0, vz = kepler.calc_orbit( - epochs_seppa[i, :], - sma[i], - ecc[i], - inc[i], - aop[i], - pan[i], - tau[i], - plx[i], - mtot[i], - tau_ref_epoch=results.tau_ref_epoch, - mass_for_Kamp=m0[i], - ) - - raoff[i, :] = raoff0 - deoff[i, :] = deoff0 - else: - raoff0, deoff0, _ = kepler.calc_orbit( - epochs_seppa[i, :], - sma[i], - ecc[i], - inc[i], - aop[i], - pan[i], - tau[i], - plx[i], - mtot[i], - tau_ref_epoch=results.tau_ref_epoch, - ) + raoff0, deoff0, vz0 = results.system.compute_all_orbits(results.post[i,:], epochs=epochs_seppa[i, :]) - raoff[i, :] = raoff0 - deoff[i, :] = deoff0 + raoff[i, :] = raoff0[:,object_to_plot,0] + deoff[i, :] = deoff0[:,object_to_plot,0] yr_epochs = Time(epochs_seppa[i, :], format="mjd").decimalyear @@ -658,13 +629,10 @@ def plot_orbits( plt.sca(ax2) plt.plot(yr_epochs, pas, color=sep_pa_color) - # plot RV orbits here - if rv_time_series == True: + # plot stellar RV orbit posterior draws here + if rv_time_series: plt.sca(ax3) - # scale back to primary RV semi amplitude - vz0 = vz * (-(mtot[i] - m0[i]) / np.median(m0[i])) - epochs_rv = np.linspace( rv_data["epoch"][0] - 3 * 365, epochs_seppa[0, -1], @@ -673,44 +641,27 @@ def plot_orbits( plt.plot( Time(epochs_rv, format="mjd").decimalyear, - vz0 + gamma3[i], + vz0[:,0,0] + gamma3[i], color=sep_pa_color, ) - if rv_time_series2 == True: - if rv_time_series == True: - plt.sca(ax4) - else: - plt.sca(ax3) - - # scale back to primary RV semi amplitude - if rv_time_series: - epochs_rv = np.linspace( - rv_data["epoch"][0] - 3 * 365, - epochs_seppa[0, -1], - num_epochs_to_plot, - ) + # plot relative RV orbit posterior draws here + if rv_time_series2: - plt.plot( - Time(epochs_rv, format="mjd").decimalyear, - vz, - color=sep_pa_color, - ) - else: - rv_data2 = results.data[results.data["object"] == 1] - rv_data2 = rv_data2[rv_data2["quant_type"] == "rv"] + rv_data2 = results.data[results.data["object"] == 1] + rv_data2 = rv_data2[rv_data2["quant_type"] == "rv"] - epochs_rv2 = np.linspace( - rv_data2["epoch"][0] - 3 * 365, - epochs_seppa[0, -1], - num_epochs_to_plot, - ) + epochs_rv2 = np.linspace( + rv_data2["epoch"][0] - 3 * 365, + epochs_seppa[0, -1], + num_epochs_to_plot, + ) - plt.plot( - Time(epochs_rv2, format="mjd").decimalyear, - vz, - color=sep_pa_color, - ) + plt.plot( + Time(epochs_rv2, format="mjd").decimalyear, + vz0[:,object_to_plot,0], + color=sep_pa_color, + ) # Plot sep/pa instruments if plot_astrometry_insts: @@ -835,21 +786,6 @@ def plot_orbits( (rv_data["instrument"] == insts[i].encode()) | (rv_data["instrument"] == insts[i]) )[0] - # choose the orbit with the best log probability - best_like = np.where(results.lnlike == np.amax(results.lnlike))[0][0] - - med_ga = [results.post[best_like, i] for i in gam_idx] - - # Get the posteriors for this index and convert to standard basis - best_post = results.basis.to_standard_basis(results.post[best_like].copy()) - - # Get the masses for the best posteriors: - best_m0 = best_post[results.standard_param_idx["m0"]] - best_m1 = best_post[ - results.standard_param_idx["m{}".format(object_to_plot)] - ] - best_mtot = best_m0 + best_m1 - # colour/shape scheme scheme for rv data points clrs = ("#0496FF", "#372554", "#FF1053", "#3A7CA5", "#143109") symbols = ("o", "^", "v", "s") @@ -863,9 +799,7 @@ def plot_orbits( rvs = inst_data["quant1"] epochs = inst_data["epoch"] epochs = Time(epochs, format="mjd").decimalyear - # don't include this so we can plot more orbits - # rvs -= med_ga[i] - # rvs -= best_post[results.param_idx[gams[i]]] + plt.scatter( epochs, rvs, @@ -888,26 +822,6 @@ def plot_orbits( else: plt.legend(fontsize=20, loc=1) - ## calculate the predicted rv trend using the best orbit - # _, _, vz = kepler.calc_orbit( - # epochs_seppa[0, :], - # best_post[results.standard_param_idx['sma{}'.format(object_to_plot)]], - # best_post[results.standard_param_idx['ecc{}'.format(object_to_plot)]], - # best_post[results.standard_param_idx['inc{}'.format(object_to_plot)]], - # best_post[results.standard_param_idx['aop{}'.format(object_to_plot)]], - # best_post[results.standard_param_idx['pan{}'.format(object_to_plot)]], - # best_post[results.standard_param_idx['tau{}'.format(object_to_plot)]], - # best_post[results.standard_param_idx['plx']], best_mtot, - # tau_ref_epoch=results.tau_ref_epoch, mass_for_Kamp=best_m0 - # ) - # - # - ## scale to the RV semiampltude of primary - # vz=vz*-(best_m1)/np.median(best_m0) - # - ## plot rv trend - # plt.plot(Time(epochs_seppa[0, :],format='mjd').decimalyear, vz, color=sep_pa_color) - if rv_time_series2 == True: if rv_time_series == False: # get list of rv instruments @@ -933,22 +847,6 @@ def plot_orbits( (rv_data["instrument"] == insts[i].encode()) | (rv_data["instrument"] == insts[i]) )[0] - # choose the orbit with the best log probability - best_like = np.where(results.lnlike == np.amax(results.lnlike))[0][0] - med_ga = [results.post[best_like, i] for i in gam_idx] - - # Get the posteriors for this index and convert to standard basis - best_post = results.basis.to_standard_basis( - results.post[best_like].copy() - ) - - # Get the masses for the best posteriors: - best_m0 = best_post[results.standard_param_idx["m0"]] - best_m1 = best_post[ - results.standard_param_idx["m{}".format(object_to_plot)] - ] - best_mtot = best_m0 + best_m1 - # colour/shape scheme scheme for rv data points clrs = ("#0496FF", "#372554", "#FF1053", "#3A7CA5", "#143109") symbols = ("o", "^", "v", "s") @@ -979,9 +877,7 @@ def plot_orbits( rvs2 = inst_data2["quant1"] epochs2 = inst_data2["epoch"] epochs2 = Time(epochs2, format="mjd").decimalyear - # don't include this so we can plot more orbits - # rvs -= med_ga[i] - # rvs -= best_post[results.param_idx[gams[i]]] + plt.scatter( epochs2, rvs2, @@ -1040,13 +936,6 @@ def plot_orbits( cbar.ax.tick_params(labelsize=15) cbar.set_label(label=cbar_param, size=20) - # hard code custom things - # ax2.set_xlim(2000, 2025) - # if rv_time_series: - # ax3.set_xlim(1995, 2025) - # if rv_time_series2: - # ax4.set_xlim(1995, 2025) - ax1.locator_params(axis="x", nbins=6) ax1.locator_params(axis="y", nbins=6) ax2.locator_params(axis="x", nbins=6) @@ -1062,15 +951,6 @@ def plot_orbits( fig.tight_layout() - # if (rv_time_series == True) and (rv_time_series2 == True): - # return fig, ax1, ax2, ax3, ax4 - # elif (rv_time_series == True) and (rv_time_series2 == False): - # return fig, ax1, ax2, ax3 - # elif (rv_time_series == False) and (rv_time_series2 == True): - # return fig, ax1, ax2, ax3 - # else: - # return fig, ax1, ax2 - return fig From ec731e179583cbc89fde3fb83480d0a751504f1f Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 7 Jul 2026 13:47:50 -0700 Subject: [PATCH 2/8] add rv + astrom e2e test --- tests/end-to-end-tests/hd4747.py | 40 ++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 tests/end-to-end-tests/hd4747.py diff --git a/tests/end-to-end-tests/hd4747.py b/tests/end-to-end-tests/hd4747.py new file mode 100644 index 00000000..29bbd367 --- /dev/null +++ b/tests/end-to-end-tests/hd4747.py @@ -0,0 +1,40 @@ +""" +Stellar RVs + astrometry, using data for HD 4747. +""" + +# Initialize Driver to Run MCMC +filename = "{}/HD4747.csv".format(DATADIR) + +num_secondary_bodies = 1 +system_mass = 0.84 # [Msol] +plx = 53.18 # [mas] +mass_err = 0.04 # [Msol] +plx_err = 0.12 # [mas] + +num_temps = 5 +num_walkers = 50 +num_threads = 50 # or a different number if you prefer + +my_driver = driver.Driver( + filename, + "MCMC", + num_secondary_bodies, + system_mass, + plx, + mass_err=mass_err, + plx_err=plx_err, + system_kwargs={"fit_secondary_mass": True, "tau_ref_epoch": 0}, + mcmc_kwargs={ + "num_temps": num_temps, + "num_walkers": num_walkers, + "num_threads": num_threads, + }, +) + +total_orbits = 100_000 +burn_steps = 1000 +thin = 2 + +# Run Sampler +m = my_driver.sampler +m.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) \ No newline at end of file From d8ce782c2f0c880971ab189a09d7a7df00dac806 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 7 Jul 2026 13:48:04 -0700 Subject: [PATCH 3/8] in-progress unit debugging script --- tests/test_rv_plotting.py | 29 +++++++++++++++++++++++------ 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/tests/test_rv_plotting.py b/tests/test_rv_plotting.py index b628a136..7ff548f4 100644 --- a/tests/test_rv_plotting.py +++ b/tests/test_rv_plotting.py @@ -1,8 +1,10 @@ -from orbitize import driver, DATADIR, plot +from orbitize import driver, DATADIR, plot, system import multiprocessing as mp +import numpy as np +import astropy.units as u, astropy.constants as cst -def test_rv_default_inst(): +def test_rv_default_inst(save_plot=False): # Initialize Driver to Run MCMC filename = "{}/HD4747.csv".format(DATADIR) @@ -13,7 +15,7 @@ def test_rv_default_inst(): plx_err = 0.12 # [mas] num_temps = 5 - num_walkers = 30 + num_walkers = 50 num_threads = mp.cpu_count() # or a different number if you prefer my_driver = driver.Driver( @@ -32,7 +34,7 @@ def test_rv_default_inst(): }, ) - total_orbits = 1000 + total_orbits = 100 burn_steps = 10 thin = 2 @@ -41,6 +43,17 @@ def test_rv_default_inst(): m.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) epochs = my_driver.system.data_table["epoch"] + # Set all posterior values equal to ones from Xuan+ 22 (to visually inspect what the plot should look like) + m.results.post[:, m.results.system.param_idx['m1']] = (67.2 * u.M_jup / u.M_sun).to('').value + m.results.post[:, m.results.system.param_idx['sma1']] = 10.0 + m.results.post[:, m.results.system.param_idx['ecc1']] = 0.7317 + m.results.post[:, m.results.system.param_idx['inc1']] = np.radians(48.0) + m.results.post[:, m.results.system.param_idx['aop1']] = np.radians(267.2) + m.results.post[:, m.results.system.param_idx['pan1']] = np.radians(89.4) + m.results.post[:, m.results.system.param_idx['tau1']] = ((2462615 - my_driver.system.tau_ref_epoch) / (33.2 * 365.25)) % 1 + m.results.post[:, m.results.system.param_idx['gamma_defrv']] = 0.0 + + # Test plotting with single orbit _ = m.results.plot_orbits( object_to_plot=1, # Plots orbits for the first (and only) companion @@ -60,6 +73,10 @@ def test_rv_default_inst(): ) + if save_plot: + import matplotlib.pyplot as plt + plt.savefig('HD4747_orbit.png') + def test_rv_multiple_inst(): filename = "{}/HR7672_joint.csv".format(DATADIR) @@ -235,7 +252,7 @@ def test_secondary_and_primary_rv(): if __name__ == "__main__": - # test_rv_default_inst() + test_rv_default_inst(save_plot=True) # test_rv_multiple_inst() - test_secondary_rv() + # test_secondary_rv() # test_secondary_and_primary_rv() From c230bd965b53f6738805b1c934687632736e92b0 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 7 Jul 2026 17:23:02 -0700 Subject: [PATCH 4/8] working e2e test --- tests/end-to-end-tests/hd4747.py | 69 ++++++++++++++++++++++++++++---- 1 file changed, 62 insertions(+), 7 deletions(-) diff --git a/tests/end-to-end-tests/hd4747.py b/tests/end-to-end-tests/hd4747.py index 29bbd367..dd244866 100644 --- a/tests/end-to-end-tests/hd4747.py +++ b/tests/end-to-end-tests/hd4747.py @@ -1,6 +1,10 @@ """ Stellar RVs + astrometry, using data for HD 4747. """ +import matplotlib.pyplot as plt +from orbitize import DATADIR, driver, plot, results +import numpy as np +from astropy import units as u, constants as cst # Initialize Driver to Run MCMC filename = "{}/HD4747.csv".format(DATADIR) @@ -13,7 +17,7 @@ num_temps = 5 num_walkers = 50 -num_threads = 50 # or a different number if you prefer +num_threads = 1 # or a different number if you prefer my_driver = driver.Driver( filename, @@ -31,10 +35,61 @@ }, ) -total_orbits = 100_000 -burn_steps = 1000 -thin = 2 +# values from Xuan et al 2022 +sma_cen = 10.0 +sma_sig = 0.2 +ecc_cen = .7317 +ecc_sig = .0014 +inc_cen = np.radians(48) +inc_sig = np.radians(1.1) +aop_cen = np.radians(267.2) +aop_sig = np.radians(0.5) +pan_cen = np.radians(89.4) +pan_sig = np.radians(1.1) +tau_cen = ((2462615 - my_driver.system.tau_ref_epoch) / (33.2 * 365.25)) % 1 +tau_sig = 0.01 + +assert my_driver.system.param_idx['sma1'] == 0 +assert my_driver.system.param_idx['ecc1'] == 1 +assert my_driver.system.param_idx['inc1'] == 2 +assert my_driver.system.param_idx['aop1'] == 3 +assert my_driver.system.param_idx['pan1'] == 4 +assert my_driver.system.param_idx['tau1'] == 5 + +walker_centres = np.array( + [sma_cen, ecc_cen, inc_cen, aop_cen, pan_cen, tau_cen] +) +walker_1sigmas = np.array( + [sma_sig, ecc_sig, inc_sig, aop_sig, pan_sig, tau_sig] +) + +new_pos = np.random.standard_normal((num_temps, num_walkers, 6)) * walker_1sigmas + walker_centres + +my_driver.sampler.curr_pos[:,:,:6] = np.copy(new_pos) +my_driver.sampler.curr_pos[:,:,my_driver.system.param_idx['m1']] = np.random.normal(loc=(67.2*u.M_jup/u.M_sun).to('').value, scale=0.01, size=(num_temps, num_walkers)) + +my_driver.sampler.check_prior_support() + +assert len(my_driver.system.rv[0]) == 56 + +total_orbits = 10_00_000 +burn_steps = 10_000 +thin = 10 + +if __name__ == '__main__': + + # Run Sampler + # m = my_driver.sampler + # m.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + # my_driver.sampler.results.save_results('hd4747_orbits.hdf5') + + my_results = results.Results() + my_results.load_results('hd4747_orbits.hdf5') + + # plot.plot_corner(my_results) + # plt.savefig('hd4747_corner.png') + + plot.plot_orbits(my_results, rv_time_series=True, start_mjd=51044.0) + plt.savefig('hd4747_orbit.png') + -# Run Sampler -m = my_driver.sampler -m.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) \ No newline at end of file From 961946741a1929a3412951581ff2a3c5dda6032b Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 7 Jul 2026 17:23:10 -0700 Subject: [PATCH 5/8] Revert "remove unnecessary code" This reverts commit 3069841f6f39ac2783b1429d370f5dcd007a6507. --- orbitize/plot.py | 164 ++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 142 insertions(+), 22 deletions(-) diff --git a/orbitize/plot.py b/orbitize/plot.py index 1628e400..26373519 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -612,10 +612,39 @@ def plot_orbits( ) # Calculate ra/dec offsets for all epochs of this orbit - raoff0, deoff0, vz0 = results.system.compute_all_orbits(results.post[i,:], epochs=epochs_seppa[i, :]) + if (rv_time_series == True) or (rv_time_series2 == True): + raoff0, deoff0, vz = kepler.calc_orbit( + epochs_seppa[i, :], + sma[i], + ecc[i], + inc[i], + aop[i], + pan[i], + tau[i], + plx[i], + mtot[i], + tau_ref_epoch=results.tau_ref_epoch, + mass_for_Kamp=m0[i], + ) + + raoff[i, :] = raoff0 + deoff[i, :] = deoff0 + else: + raoff0, deoff0, _ = kepler.calc_orbit( + epochs_seppa[i, :], + sma[i], + ecc[i], + inc[i], + aop[i], + pan[i], + tau[i], + plx[i], + mtot[i], + tau_ref_epoch=results.tau_ref_epoch, + ) - raoff[i, :] = raoff0[:,object_to_plot,0] - deoff[i, :] = deoff0[:,object_to_plot,0] + raoff[i, :] = raoff0 + deoff[i, :] = deoff0 yr_epochs = Time(epochs_seppa[i, :], format="mjd").decimalyear @@ -629,10 +658,13 @@ def plot_orbits( plt.sca(ax2) plt.plot(yr_epochs, pas, color=sep_pa_color) - # plot stellar RV orbit posterior draws here - if rv_time_series: + # plot RV orbits here + if rv_time_series == True: plt.sca(ax3) + # scale back to primary RV semi amplitude + vz0 = vz * (-(mtot[i] - m0[i]) / np.median(m0[i])) + epochs_rv = np.linspace( rv_data["epoch"][0] - 3 * 365, epochs_seppa[0, -1], @@ -641,27 +673,44 @@ def plot_orbits( plt.plot( Time(epochs_rv, format="mjd").decimalyear, - vz0[:,0,0] + gamma3[i], + vz0 + gamma3[i], color=sep_pa_color, ) - # plot relative RV orbit posterior draws here - if rv_time_series2: + if rv_time_series2 == True: + if rv_time_series == True: + plt.sca(ax4) + else: + plt.sca(ax3) + + # scale back to primary RV semi amplitude + if rv_time_series: + epochs_rv = np.linspace( + rv_data["epoch"][0] - 3 * 365, + epochs_seppa[0, -1], + num_epochs_to_plot, + ) - rv_data2 = results.data[results.data["object"] == 1] - rv_data2 = rv_data2[rv_data2["quant_type"] == "rv"] + plt.plot( + Time(epochs_rv, format="mjd").decimalyear, + vz, + color=sep_pa_color, + ) + else: + rv_data2 = results.data[results.data["object"] == 1] + rv_data2 = rv_data2[rv_data2["quant_type"] == "rv"] - epochs_rv2 = np.linspace( - rv_data2["epoch"][0] - 3 * 365, - epochs_seppa[0, -1], - num_epochs_to_plot, - ) + epochs_rv2 = np.linspace( + rv_data2["epoch"][0] - 3 * 365, + epochs_seppa[0, -1], + num_epochs_to_plot, + ) - plt.plot( - Time(epochs_rv2, format="mjd").decimalyear, - vz0[:,object_to_plot,0], - color=sep_pa_color, - ) + plt.plot( + Time(epochs_rv2, format="mjd").decimalyear, + vz, + color=sep_pa_color, + ) # Plot sep/pa instruments if plot_astrometry_insts: @@ -786,6 +835,21 @@ def plot_orbits( (rv_data["instrument"] == insts[i].encode()) | (rv_data["instrument"] == insts[i]) )[0] + # choose the orbit with the best log probability + best_like = np.where(results.lnlike == np.amax(results.lnlike))[0][0] + + med_ga = [results.post[best_like, i] for i in gam_idx] + + # Get the posteriors for this index and convert to standard basis + best_post = results.basis.to_standard_basis(results.post[best_like].copy()) + + # Get the masses for the best posteriors: + best_m0 = best_post[results.standard_param_idx["m0"]] + best_m1 = best_post[ + results.standard_param_idx["m{}".format(object_to_plot)] + ] + best_mtot = best_m0 + best_m1 + # colour/shape scheme scheme for rv data points clrs = ("#0496FF", "#372554", "#FF1053", "#3A7CA5", "#143109") symbols = ("o", "^", "v", "s") @@ -799,7 +863,9 @@ def plot_orbits( rvs = inst_data["quant1"] epochs = inst_data["epoch"] epochs = Time(epochs, format="mjd").decimalyear - + # don't include this so we can plot more orbits + # rvs -= med_ga[i] + # rvs -= best_post[results.param_idx[gams[i]]] plt.scatter( epochs, rvs, @@ -822,6 +888,26 @@ def plot_orbits( else: plt.legend(fontsize=20, loc=1) + ## calculate the predicted rv trend using the best orbit + # _, _, vz = kepler.calc_orbit( + # epochs_seppa[0, :], + # best_post[results.standard_param_idx['sma{}'.format(object_to_plot)]], + # best_post[results.standard_param_idx['ecc{}'.format(object_to_plot)]], + # best_post[results.standard_param_idx['inc{}'.format(object_to_plot)]], + # best_post[results.standard_param_idx['aop{}'.format(object_to_plot)]], + # best_post[results.standard_param_idx['pan{}'.format(object_to_plot)]], + # best_post[results.standard_param_idx['tau{}'.format(object_to_plot)]], + # best_post[results.standard_param_idx['plx']], best_mtot, + # tau_ref_epoch=results.tau_ref_epoch, mass_for_Kamp=best_m0 + # ) + # + # + ## scale to the RV semiampltude of primary + # vz=vz*-(best_m1)/np.median(best_m0) + # + ## plot rv trend + # plt.plot(Time(epochs_seppa[0, :],format='mjd').decimalyear, vz, color=sep_pa_color) + if rv_time_series2 == True: if rv_time_series == False: # get list of rv instruments @@ -847,6 +933,22 @@ def plot_orbits( (rv_data["instrument"] == insts[i].encode()) | (rv_data["instrument"] == insts[i]) )[0] + # choose the orbit with the best log probability + best_like = np.where(results.lnlike == np.amax(results.lnlike))[0][0] + med_ga = [results.post[best_like, i] for i in gam_idx] + + # Get the posteriors for this index and convert to standard basis + best_post = results.basis.to_standard_basis( + results.post[best_like].copy() + ) + + # Get the masses for the best posteriors: + best_m0 = best_post[results.standard_param_idx["m0"]] + best_m1 = best_post[ + results.standard_param_idx["m{}".format(object_to_plot)] + ] + best_mtot = best_m0 + best_m1 + # colour/shape scheme scheme for rv data points clrs = ("#0496FF", "#372554", "#FF1053", "#3A7CA5", "#143109") symbols = ("o", "^", "v", "s") @@ -877,7 +979,9 @@ def plot_orbits( rvs2 = inst_data2["quant1"] epochs2 = inst_data2["epoch"] epochs2 = Time(epochs2, format="mjd").decimalyear - + # don't include this so we can plot more orbits + # rvs -= med_ga[i] + # rvs -= best_post[results.param_idx[gams[i]]] plt.scatter( epochs2, rvs2, @@ -936,6 +1040,13 @@ def plot_orbits( cbar.ax.tick_params(labelsize=15) cbar.set_label(label=cbar_param, size=20) + # hard code custom things + # ax2.set_xlim(2000, 2025) + # if rv_time_series: + # ax3.set_xlim(1995, 2025) + # if rv_time_series2: + # ax4.set_xlim(1995, 2025) + ax1.locator_params(axis="x", nbins=6) ax1.locator_params(axis="y", nbins=6) ax2.locator_params(axis="x", nbins=6) @@ -951,6 +1062,15 @@ def plot_orbits( fig.tight_layout() + # if (rv_time_series == True) and (rv_time_series2 == True): + # return fig, ax1, ax2, ax3, ax4 + # elif (rv_time_series == True) and (rv_time_series2 == False): + # return fig, ax1, ax2, ax3 + # elif (rv_time_series == False) and (rv_time_series2 == True): + # return fig, ax1, ax2, ax3 + # else: + # return fig, ax1, ax2 + return fig From cee01e477832470a5b83132bdcb098c8840d19c3 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 7 Jul 2026 17:23:56 -0700 Subject: [PATCH 6/8] Revert "in-progress unit debugging script" This reverts commit d8ce782c2f0c880971ab189a09d7a7df00dac806. --- tests/test_rv_plotting.py | 29 ++++++----------------------- 1 file changed, 6 insertions(+), 23 deletions(-) diff --git a/tests/test_rv_plotting.py b/tests/test_rv_plotting.py index 7ff548f4..b628a136 100644 --- a/tests/test_rv_plotting.py +++ b/tests/test_rv_plotting.py @@ -1,10 +1,8 @@ -from orbitize import driver, DATADIR, plot, system +from orbitize import driver, DATADIR, plot import multiprocessing as mp -import numpy as np -import astropy.units as u, astropy.constants as cst -def test_rv_default_inst(save_plot=False): +def test_rv_default_inst(): # Initialize Driver to Run MCMC filename = "{}/HD4747.csv".format(DATADIR) @@ -15,7 +13,7 @@ def test_rv_default_inst(save_plot=False): plx_err = 0.12 # [mas] num_temps = 5 - num_walkers = 50 + num_walkers = 30 num_threads = mp.cpu_count() # or a different number if you prefer my_driver = driver.Driver( @@ -34,7 +32,7 @@ def test_rv_default_inst(save_plot=False): }, ) - total_orbits = 100 + total_orbits = 1000 burn_steps = 10 thin = 2 @@ -43,17 +41,6 @@ def test_rv_default_inst(save_plot=False): m.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) epochs = my_driver.system.data_table["epoch"] - # Set all posterior values equal to ones from Xuan+ 22 (to visually inspect what the plot should look like) - m.results.post[:, m.results.system.param_idx['m1']] = (67.2 * u.M_jup / u.M_sun).to('').value - m.results.post[:, m.results.system.param_idx['sma1']] = 10.0 - m.results.post[:, m.results.system.param_idx['ecc1']] = 0.7317 - m.results.post[:, m.results.system.param_idx['inc1']] = np.radians(48.0) - m.results.post[:, m.results.system.param_idx['aop1']] = np.radians(267.2) - m.results.post[:, m.results.system.param_idx['pan1']] = np.radians(89.4) - m.results.post[:, m.results.system.param_idx['tau1']] = ((2462615 - my_driver.system.tau_ref_epoch) / (33.2 * 365.25)) % 1 - m.results.post[:, m.results.system.param_idx['gamma_defrv']] = 0.0 - - # Test plotting with single orbit _ = m.results.plot_orbits( object_to_plot=1, # Plots orbits for the first (and only) companion @@ -73,10 +60,6 @@ def test_rv_default_inst(save_plot=False): ) - if save_plot: - import matplotlib.pyplot as plt - plt.savefig('HD4747_orbit.png') - def test_rv_multiple_inst(): filename = "{}/HR7672_joint.csv".format(DATADIR) @@ -252,7 +235,7 @@ def test_secondary_and_primary_rv(): if __name__ == "__main__": - test_rv_default_inst(save_plot=True) + # test_rv_default_inst() # test_rv_multiple_inst() - # test_secondary_rv() + test_secondary_rv() # test_secondary_and_primary_rv() From 4a1af9046ff15781a3b47aa2dc6a73f42d3c1378 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 7 Jul 2026 17:25:02 -0700 Subject: [PATCH 7/8] comment in full test --- tests/end-to-end-tests/hd4747.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/end-to-end-tests/hd4747.py b/tests/end-to-end-tests/hd4747.py index dd244866..458ea650 100644 --- a/tests/end-to-end-tests/hd4747.py +++ b/tests/end-to-end-tests/hd4747.py @@ -79,15 +79,15 @@ if __name__ == '__main__': # Run Sampler - # m = my_driver.sampler - # m.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) - # my_driver.sampler.results.save_results('hd4747_orbits.hdf5') + m = my_driver.sampler + m.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin) + my_driver.sampler.results.save_results('hd4747_orbits.hdf5') my_results = results.Results() my_results.load_results('hd4747_orbits.hdf5') - # plot.plot_corner(my_results) - # plt.savefig('hd4747_corner.png') + plot.plot_corner(my_results) + plt.savefig('hd4747_corner.png') plot.plot_orbits(my_results, rv_time_series=True, start_mjd=51044.0) plt.savefig('hd4747_orbit.png') From e6bbd4e4ab1acecd544a870d9d43024fc265ae04 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 7 Jul 2026 17:26:52 -0700 Subject: [PATCH 8/8] fix rv plotting bug --- orbitize/plot.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/orbitize/plot.py b/orbitize/plot.py index 26373519..afaf85c1 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -672,7 +672,7 @@ def plot_orbits( ) plt.plot( - Time(epochs_rv, format="mjd").decimalyear, + Time(epochs_seppa[i, :], format="mjd").decimalyear, vz0 + gamma3[i], color=sep_pa_color, ) @@ -692,7 +692,7 @@ def plot_orbits( ) plt.plot( - Time(epochs_rv, format="mjd").decimalyear, + Time(epochs_seppa[i, :], format="mjd").decimalyear, vz, color=sep_pa_color, ) @@ -707,7 +707,7 @@ def plot_orbits( ) plt.plot( - Time(epochs_rv2, format="mjd").decimalyear, + Time(epochs_seppa[i, :], format="mjd").decimalyear, vz, color=sep_pa_color, )