Skip to content

Commit

Permalink
Merge pull request #280 from sblunt/rebound_naming
Browse files Browse the repository at this point in the history 
updated function names in test_rebound
  • Loading branch information
@sblunt
sblunt committed Oct 13, 2021
2 parents fb7a912 + d62dc87 commit 90f1f55
Showing 1 changed file with 53 additions and 67 deletions.
120 changes: 53 additions & 67 deletions tests/test_rebound.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@
data_table = read_file(input_file)
num_secondary_bodies = 4

epochs = Time(np.linspace(2020, 2025, num=int(1000)), format='decimalyear').mjd
epochs = Time(np.linspace(2020, 2022, num=int(1000)), format='decimalyear').mjd

sma1 = sma[0]
ecc1 = ecc[0]
Expand Down Expand Up @@ -79,96 +79,82 @@
m1, m2, m3, m4, m_st
])

# these arrays have shape (n_epochs x n_bodies x 1)
ra, dec, _ = hr8799_sys.compute_all_orbits(params_arr, epochs, comp_rebound=True)
def test_8799_rebound_vs_kepler(plotname=None):
"""
Test that for the HR 8799 system between 2020 and 2022, the rebound and
keplerian (with epicycle approximation) backends produce similar results.
"""


def calc_diff():
import matplotlib.pyplot as plt
from astropy.time import Time
assert hr8799_sys.track_planet_perturbs

rra, rde, rvz = hr8799_sys.compute_all_orbits(params_arr, epochs=epochs, comp_rebound=True)
kra, kde, kvz = hr8799_sys.compute_all_orbits(params_arr, epochs=epochs, comp_rebound=False)
#ora = ra
#odec = dec
rra, rde, _ = hr8799_sys.compute_all_orbits(params_arr, epochs=epochs, comp_rebound=True)
kra, kde, _ = hr8799_sys.compute_all_orbits(params_arr, epochs=epochs, comp_rebound=False)

delta_ra = abs(rra-kra[:,:,0])
delta_de = abs(rde-kde[:,:,0])
#delta_vz = abs(rvz-_)
yepochs = Time(epochs, format='mjd').decimalyear

if len(sma)==1:
plt.plot(yepochs, delta_ra, label = 'Planet X: RA offsets')
plt.plot(yepochs, delta_de, label = 'Planet X: Dec offsets')
#plt.plot(yepochs, delta_vz, label = 'Planet X: RV offsets')
# check that the difference between these two solvers is smaller than
# GRAVITY precision over the next two years
gravity_precision_mas = 100 * 1e-3
assert np.all(delta_ra < gravity_precision_mas)
assert np.all(delta_de < gravity_precision_mas)

elif len(sma)==4:

# check that the difference between these two solvers is larger than
# 10 uas, though
assert np.max(delta_ra) > 10 * 1e-3
assert np.max(delta_de) > 10 * 1e-3

if plotname is not None:

# plot 1: RA & Dec differences over time
fig, (ax1, ax2) = plt.subplots(2)
fig.suptitle('Massive Orbits in Rebound vs. Orbitize approx.')

ax1.plot(yepochs, delta_ra[:,0], 'black', label = 'Star') #fourth planet
ax1.plot(yepochs, delta_ra[:,0], 'black', label = 'Star')
ax2.plot(yepochs, delta_de[:,0], 'dimgray', label = 'Star')
#plt.plot(yepochs, delta_vz[:,0], 'silver', label = 'Star')

ax1.plot(yepochs, delta_ra[:,1], 'brown', label = 'Planet E: RA offsets') #first planet
ax2.plot(yepochs, delta_de[:,1], 'red', label = 'Planet E: Dec offsets')
#plt.plot(yepochs, delta_vz[:,1], 'pink', label = 'Planet B: RV offsets')

ax1.plot(yepochs, delta_ra[:,2], 'coral', label = 'Planet D: RA offsets') #second planet
ax2.plot(yepochs, delta_de[:,2], 'orange', label = 'Planet D: Dec offsets')
#plt.plot(yepochs, delta_vz[:,2], 'gold', label = 'Planet C: RV offsets')

ax1.plot(yepochs, delta_ra[:,3], 'greenyellow', label = 'Planet C: RA offsets') #third planet
ax2.plot(yepochs, delta_de[:,3], 'green', label = 'Planet C: Dec offsets')
#plt.plot(yepochs, delta_vz[:,3], 'darkgreen', label = 'Planet D: RV offsets')


ax1.plot(yepochs, delta_ra[:,4], 'dodgerblue', label = 'Planet B: RA offsets') #fourth planet
ax2.plot(yepochs, delta_de[:,4], 'blue', label = 'Planet B: Dec offsets')
#plt.plot(yepochs, delta_vz[:,4], 'indigo', label = 'Planet E: RV offsets')




else:
print('I dont feel like it')

plt.xlabel('year')
plt.ylabel('milliarcseconds')
ax1.legend()
ax2.legend()
plt.savefig('foo.png')

def plot_orbit():

kra, kde, kvz = kepler.calc_orbit(epochs,sma,ecc,inc,aop,pan,tau,plx,mtot,tau_ref_epoch=tau_ref_epoch)
rra, rdec, rvz = nbody.calc_orbit(epochs,sma,ecc,inc,aop,pan,tau,plx,mtot,tau_ref_epoch=tau_ref_epoch)

plt.plot(kra[:,1:5], kde[:,1:5], 'indigo', label = 'Orbitize approx.')
plt.plot(kra[-1,1:5], kde[-1,1:5],'o')

plt.plot(rra, rdec, 'r', label = 'Rebound', alpha = 0.25)
plt.plot(rra[-1], rdec[-1], 'o', alpha = 0.25)
plt.xlabel('year')
plt.ylabel('milliarcseconds')
ax1.legend()
ax2.legend()
plt.savefig('{}_differences.png'.format(plotname), dpi=250)

# plot 2: orbit tracks over time
plt.figure()
plt.plot(kra[:,1:5,0], kde[:,1:5,0], 'indigo', label = 'Orbitize approx.')
plt.plot(kra[-1,1:5,0], kde[-1,1:5,0],'o')

plt.plot(0, 0, '*')
plt.legend()
plt.savefig('foo.png')

def plot_orbit2():

rra, rde, rvz = hr8799_sys.compute_all_orbits(params_arr, epochs=epochs, comp_rebound=True)
kra, kde, kvz = hr8799_sys.compute_all_orbits(params_arr, epochs=epochs, comp_rebound=False)

plt.plot(kra[:,0], kde[:,0], 'indigo', label = 'Orbitize approx.')
plt.plot(kra[-1,0], kde[-1,0],'o')

plt.plot(rra[:,0], rde[:,0], 'r', label = 'Rebound', alpha = 0.25)
plt.plot(rra[-1,0], rde[-1,0], 'o', alpha = 0.25)
plt.plot(rra, rde, 'r', label = 'Rebound', alpha = 0.25)
plt.plot(rra[-1], rde[-1], 'o', alpha = 0.25)

plt.plot(0, 0, '*')
plt.legend()
plt.savefig('{}_orbittracks.png'.format(plotname), dpi=250)

# plot 3: primary orbit track over time
plt.figure()
plt.plot(kra[:,0,0], kde[:,0,0], 'indigo', label = 'Orbitize approx.')
plt.plot(kra[-1,0,0], kde[-1,0,0],'o')

plt.plot(0, 0, '*')
plt.legend()
plt.savefig('foo.png')

calc_diff()
# plot_orbit()
plt.plot(rra[:,0], rde[:,0], 'r', label = 'Rebound', alpha = 0.25)
plt.plot(rra[-1,0], rde[-1,0], 'o', alpha = 0.25)

plt.plot(0, 0, '*')
plt.legend()
plt.savefig('{}_primaryorbittrack.png'.format(plotname), dpi=250)

if __name__=='__main__':
test_8799_rebound_vs_kepler(plotname='hr8799_diffs')

0 comments on commit 90f1f55

Please sign in to comment.