keplerianFrequency formula
From a pure Celestial Mechanics point of view, omega = std::sqrt(coupleMass/std::pow(std::abs(coupleMass/(2*energy)),3))
does not seem correct to me. The Keplerian Energy is E=-\mathcal{G} M_* m_{pl}/(2a)
so a=|\mathcal{G} M_* m_{pl}/(2E)|
(assuming E<0
as is the case for elliptical orbits). The Keplerian frequency is \Omega_K = \sqrt{\mathcal{G}(M_* + m_{pl})/a^3}
. So I think the problem is that it's not omega = std::sqrt(coupleMass/std::pow(std::abs(coupleMass/(2*energy)),3))
but rather omega = std::sqrt(coupleMass/std::pow(std::abs(mass()*XMSTAR/(2*energy)),3))
(I'm guessing here mass()
is the planet's mass; also \mathcal{G}=1
so it's dropped).
(If I could have it my way, I'd write:
Planet::keplerianFrequency() const {
real mu = XMSTAR+mass();
real kKepl= XMSTAR*mass();
real distance = position().norm();
real energy = (kKepl/mu)*velocity().norm2()/2 - kKepl/distance;
real omega = std::sqrt(mu/std::pow(std::abs(kKepl/(2*energy)),3));
return omega;
}
This notation would be natural for celestial mechanicians: the name mu
is from the Standard gravitational parameter; the name kKepl
is from the Keplerian potential in the Kepler problem
)
EDIT: I modified the code snippet to make it correct. Like this the energy
really is the mechanical energy. However, this turns out to be equivalent to the legacy code, which used the name energy
for a rescaled energy. Except for the slight abuse of nomenclature, the legacy code is correct.