Skip to content

Commit

Permalink
planekeplerian udpla now constructed from various element types
Browse files Browse the repository at this point in the history
  • Loading branch information
darioizzo committed Sep 20, 2023
1 parent 936a158 commit a6f3f34
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 12 deletions.
6 changes: 2 additions & 4 deletions include/kep3/planet.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -231,10 +231,8 @@ struct kep3_DLL_PUBLIC_INLINE_CLASS planet_inner final : planet_inner_base {
double v2 = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
double en = v2 / 2. - mu / R;
if (en > 0) {
// If the energy is positive we have an hyperbolae
throw std::logic_error(
"Hyperbolic conditions detected no period available for '" +
get_name() + "'");
// If the energy is positive we have an hyperbolae and we return nan
return std::numeric_limits<double>::quiet_NaN();
} else {
double a = -mu / 2. / en;
return kep3::pi * 2. * std::sqrt(a * a * a / mu);
Expand Down
3 changes: 1 addition & 2 deletions include/kep3/planets/keplerian.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -49,11 +49,10 @@ class kep3_DLL_PUBLIC keplerian {
}

public:
// NOTE: in here par is a,e,i,W,w,M (Mean anomaly, not true anomaly)
// NOTE: added_param contains mu_self, radius and safe_radius
explicit keplerian(const epoch &ref_epoch, const std::array<double, 6> &par,
double mu_central_body = 1., std::string name = "Unknown",
std::array<double, 3> added_params = {-1., -1., -1.});
std::array<double, 3> added_params = {-1., -1., -1.}, kep3::elements_type el_t = kep3::elements_type::KEP_F);
// Constructor from pos_vel
explicit keplerian(const epoch &ref_epoch = kep3::epoch(0),
const std::array<std::array<double, 3>, 2> &pos_vel =
Expand Down
46 changes: 41 additions & 5 deletions src/planets/keplerian.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,9 @@
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "kep3/core_astro/eq2par2eq.hpp"
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>

Expand Down Expand Up @@ -41,27 +44,60 @@ keplerian::keplerian(const epoch &ref_epoch,
mu_central_body / R;
(en > 0) ? m_ellipse = false : m_ellipse = true;
double a = -m_mu_central_body / 2. / en;
m_period = kep3::pi * 2. * std::sqrt(a * a * a / m_mu_central_body);
if (m_ellipse) {
m_period = kep3::pi * 2. * std::sqrt(a * a * a / m_mu_central_body);
} else {
m_period = std::numeric_limits<double>::quiet_NaN();
}
}

keplerian::keplerian(const epoch &ref_epoch, const std::array<double, 6> &par,
keplerian::keplerian(const epoch &ref_epoch,
const std::array<double, 6> &par_in,
double mu_central_body, std::string name,
std::array<double, 3> added_params)
std::array<double, 3> added_params,
kep3::elements_type el_type)
: m_ref_epoch(ref_epoch), m_pos_vel_0(), m_name(std::move(name)),
m_mu_central_body(mu_central_body), m_mu_self(added_params[0]),
m_radius(added_params[1]), m_safe_radius(added_params[2]), m_period(),
m_ellipse() {
// orbital parameters a,e,i,W,w,f will be stored here
std::array<double, 6> par(par_in);
// we convert according to the chosen input
switch (el_type) {
case kep3::elements_type::KEP_F:
break;
case kep3::elements_type::KEP_M:
if (par[0] < 0) {
throw std::logic_error("Mean anomaly is only available for ellipses.");
}
par[5] = kep3::m2f(par[5], par[1]);
break;
case kep3::elements_type::MEQ:
par = kep3::eq2par(par, false);
break;
case kep3::elements_type::MEQ_R:
par = kep3::eq2par(par, true);
break;
default:
throw std::logic_error("You should not go here!");
}

if (par[0] * (1 - par[1]) <= 0) {
throw std::domain_error(
"A Keplerian planet constructor was called with with non compatible "
"a,e:"
"The following must hold: a<0 -> e>1 [a>0 -> e<1].");
}

m_pos_vel_0 = kep3::par2ic(par, mu_central_body);

(par[0] < 0) ? m_ellipse = false : m_ellipse = true;
m_period =
kep3::pi * 2. * std::sqrt(par[0] * par[0] * par[0] / m_mu_central_body);
if (m_ellipse) {
m_period =
kep3::pi * 2. * std::sqrt(par[0] * par[0] * par[0] / m_mu_central_body);
} else {
m_period = std::numeric_limits<double>::quiet_NaN();
}
}

std::array<std::array<double, 3>, 2>
Expand Down
2 changes: 1 addition & 1 deletion test/planet_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -189,7 +189,7 @@ TEST_CASE("construction") {
simple_udpla_mu_h udpla{};
REQUIRE_NOTHROW(planet(udpla));
planet pla(udpla);
REQUIRE_THROWS_AS((pla.period()), std::logic_error);
REQUIRE(!std::isfinite(pla.period()));
}
{
// Constructor from a simple udpla with mu and elliptical orbit
Expand Down

0 comments on commit a6f3f34

Please sign in to comment.