m_planet_units variable not changing planet m_planet units in exoplanet.orbits.KeplerianOrbit call

[dyahalomi]

Describe the bug
I updated my exoplanet version from <=0.4.x to 0.5.1 and noticed a several order of magnitudes change to the RV signal amplitude produced using exoplanet.orbits.KeplerianOrbit. Going into my code, I was able to fix the issue by changing my units for m_planet to M_sun rather that M_earth. However, I set in my exoplanet.orbits.KeplerianOrbit function the variable m_planet_units=M_earth (where M_earth is imported by "from astropy.constants import M_earth"), but this doesn't seem to be working.

To Reproduce
Run exoplanet.orbits.KeplerianOrbit with an m_planet in Earth units and m_planet_units=M_earth. The units of m_planet still appear to be in M_sun units.

Expected behavior
I expected that calling m_planet_units=M_earth would change the planet units to m_earth.

Your setup (please complete the following information):

• Version of exoplanet: 0.5.1
• Operating system: macOS version 11.4
• Python version & installation method (pip, conda, etc.): python version 3.9.5. exoplanet and other packages installed with pip.

Additional context

[dfm]

Thanks! Can you write a small block of code that executes and raises an exception if the results aren't what you expect? Then we can turn it into a unit test.

[dyahalomi]

Sure, this should work -- thanks! I also added a plotting bit below the unit test, if helpful to see visually.

---

import numpy as np
import matplotlib.pyplot as plt
import exoplanet as xo
from astropy.constants import M_earth


P_earth = 365.256 #days
e_earth = 0.0167 
Tper_earth= 2454115.5208333 #BJD
omega_earth = np.radians(102.9)
Omega_earth = np.radians(0.0)
inclination_earth = np.radians(45.0)
m_earth = 1 #units m_earth


orbit = xo.orbits.KeplerianOrbit(
    period = P_earth,
    ecc = e_earth, 
    t_periastron = Tper_earth, 
    omega = omega_earth, 
    Omega = Omega_earth,  
    incl = inclination_earth,
    m_planet = m_earth, 
    m_planet_units = M_earth)


times_rv = np.linspace(Tper_earth, Tper_earth+1000, 10000) 
rv = orbit.get_radial_velocity(times_rv)
rv_orbit = rv.eval()


#Earth RV signal around sun should only be around 0.1 m/s, to give some flexibility here let's increase this
#by an order of magnitude and say that the max amplitude in rv_orbit minus the min amplitude in rv_orbit
#in units of m/s cannot excede 1 m/s
rv_diff = np.max(rv_orbit)-np.min(rv_orbit)
assert(rv_diff < 1.0 ), "ERROR! THIS IS NOT EARTH LIKE! max amplitude in rv_orbit minus the min amplitude in rv_orbit is greater than 1.0 m/s."

### if you get past the assertion error, then plot RV plot
fig, ax = plt.subplots(1, figsize = [18,10])
fig.suptitle("RV Signal", fontsize = 45)


ax.plot(times_rv, rv_orbit, color = 'k')

ax.set_xlabel("time [BJD]", fontsize = 27)
ax.set_ylabel("RV [m/s]", fontsize = 27)

fig.tight_layout()
plt.show()