Skip to content

Commit

Permalink
Merge pull request #86 from darioizzo/master
Browse files Browse the repository at this point in the history 
Equinoctial elements work.
  • Loading branch information
@darioizzo
darioizzo authored Jan 29, 2019
2 parents bae1006 + 512f0ee commit 04638d5
Show file tree
Hide file tree
Showing 5 changed files with 63 additions and 63 deletions.
8 changes: 4 additions & 4 deletions pykep/core/core.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -590,14 +590,14 @@ BOOST_PYTHON_MODULE(_core)
"- v: velocity (cartesian)\n"
"- mu: gravity parameter\n\n"
"- retrogade: uses the retrograde parameters. Default value is False.\n\n"
"Returns the modified equinoctial elements a(1-e^2),h,k,p,q,L\n"
"Returns the modified equinoctial elements a(1-e^2),f,g,h,k,L\n"
"L is the true mean longitude\n"
"Example:: \n\n"
" E = ic2eq(r = [1,0,0], v = [0,1,0], mu =1.0)",
(arg("r"), arg("v"), arg("mu") = 1.0, arg("retrograde") = false));

def("eq2ic", &eq2ic_wrapper, "pykep.eq2ic(EQ,mu = 1.0, retrograde = False)\n\n"
"- EQ: modified equinoctial elements a(1-e^2),h,k,p,q,L\n"
"- EQ: modified equinoctial elements a(1-e^2),f,g,h,k,L\n"
"- mu: gravity parameter (l^3/s^2)\n\n"
"Returns cartesian elements from Keplerian elements\n"
"L is the true longitude\n"
Expand All @@ -606,7 +606,7 @@ BOOST_PYTHON_MODULE(_core)
(arg("eq"), arg("mu") = 1.0, arg("retrograde") = false));

def("eq2par", &eq2par_wrapper, "pykep.eq2par(EQ, retrograde = False)\n\n"
"- EQ: modified equinoctial elements a(1-e^2),h,k,p,q,L\n"
"- EQ: modified equinoctial elements a(1-e^2),f,g,h,k,L\n"
"- retrogade: uses the retrograde parameters. Default value is False.\n\n"
"Returns the osculating Keplerian elements a,e,i,W,w,E \n"
"E is the eccentric anomaly for e<1, the Gudermannian for e>1\n"
Expand All @@ -619,7 +619,7 @@ BOOST_PYTHON_MODULE(_core)
def("par2eq", &par2eq_wrapper, "pykep.par2eq(E, retrograde = False)\n\n"
"- E: osculating keplerian elements a,e,i,W,w,E ( l, ND, rad, rad, rad, rad)\n"
"- retrogade: uses the retrograde parameters. Default value is False.\n\n"
"Returns the modified equinoctial elements a(1-e^2),h,k,p,q,L\n"
"Returns the modified equinoctial elements a(1-e^2),f,g,h,k,L\n"
"L is the true mean longitude\n"
"Example:: \n\n"
" E = pk.par2eq(E = [1.0,0.1,0.2,0.3,0.4,0.5], retrograde = False)",
Expand Down
40 changes: 20 additions & 20 deletions src/core_functions/eq2ic.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,10 +34,10 @@ namespace kep_toolbox

/// From equinoctial to cartesian
/**
* Transforms equinoctial elements (a,h,k,p,q,L) to cartesian.
* Transforms equinoctial elements (p,f,g,h,k,L) to cartesian.
* Note that we use the true longitude LF (i.e. L + om + Om)
*
* @param[in] EQ the the equinoctial elements (p,h,k,p,q,L) in SI.
* @param[in] EQ the the equinoctial elements (p,f,g,h,k,L) in SI.
* @param[in] mu gravitational parameter, (m^3/s^2)
*
* @param[out] r the cartesian position corresponding to the input elements
Expand All @@ -58,34 +58,34 @@ void eq2ic(const vettore6D &EQ, const double &mu, vettore3D &r, vettore3D &v, co
// computations
// to make sense
auto par = std::abs(EQ[0]);
auto h = EQ[1];
auto k = EQ[2];
auto p = EQ[3];
auto q = EQ[4];
auto f = EQ[1];
auto g = EQ[2];
auto h = EQ[3];
auto k = EQ[4];
auto L = EQ[5];

// We compute the equinoctial reference frame
auto den = p * p + q * q + 1;
auto fx = (1 - p * p + q * q) / den;
auto fy = (2 * p * q) / den;
auto fz = (-2 * I * p) / den;
auto den = k * k + h * h + 1;
auto fx = (1 - k * k + h * h) / den;
auto fy = (2 * k * h) / den;
auto fz = (-2 * I * k) / den;

auto gx = (2 * I * p * q) / den;
auto gy = (1 + p * p - q * q) * I / den;
auto gz = (2 * q) / den;
auto gx = (2 * I * k * h) / den;
auto gy = (1 + k * k - h * h) * I / den;
auto gz = (2 * h) / den;

auto wx = (2 * p) / den;
auto wy = (-2 * q) / den;
auto wz = (1 - p * p - q * q) * I / den;
auto wx = (2 * k) / den;
auto wy = (-2 * h) / den;
auto wz = (1 - k * k - h * h) * I / den;

// Auxiliary
auto b = std::sqrt(1 - h * h - k * k);
auto radius = par / (1 + h * std::sin(L) + k * std::cos(L));
auto b = std::sqrt(1 - g * g - f * f);
auto radius = par / (1 + g * std::sin(L) + f * std::cos(L));
// In the equinoctial reference frame
auto X = radius * std::cos(L);
auto Y = radius * std::sin(L);
auto VX = -std::sqrt(mu / par) * (h + std::sin(L));
auto VY = std::sqrt(mu / par) * (k + std::cos(L));
auto VX = -std::sqrt(mu / par) * (g + std::sin(L));
auto VY = std::sqrt(mu / par) * (f + std::cos(L));

/*print("r is: ", radius, "\n");
print("f is: [", fx, ", ", fy, ", ", fz, "]\n");
Expand Down
8 changes: 4 additions & 4 deletions src/core_functions/eq2par.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,12 +36,12 @@ namespace kep_toolbox

/// From modified equinoctial elements to osculating parameters
/**
* Transforms modified equinoctial elements (p,h,k,p,q,L) to osculating Keplerian parameters (a,e,i,W,w,E).
* Transforms modified equinoctial elements (p,f,g,h,k,L) to osculating Keplerian parameters (a,e,i,W,w,E).
* Note that we use the eccentric anomaly E (or gudermannian) and the true longitude L (i.e. f + om + Om).
*
* Note that nan will be generated if the osculating parameters are not defined (e=0, i=0, etc..)
*
* @param[in] EQ the equinoctial elements (p,h,k,p,q,L), anomaly in rad.
* @param[in] EQ the equinoctial elements (p,f,g,h,k,L), anomaly in rad.
* @param[in] retrogade forces retrograde elements to be used
*
* @param[out] E the osculating Keplerian parameters (a,e,i,W,w,E).
Expand All @@ -56,13 +56,13 @@ void eq2par(vettore6D &E, const vettore6D &EQ, const bool retrogade = false)
}
auto ecc = std::sqrt(EQ[1] * EQ[1] + EQ[2] * EQ[2]);
auto tmp = std::sqrt(EQ[3] * EQ[3] + EQ[4] * EQ[4]);
auto zita = std::atan2(EQ[1] / ecc, EQ[2] / ecc);
auto zita = std::atan2(EQ[2] / ecc, EQ[1] / ecc);

E[1] = ecc;
E[0] = EQ[0] / (1 - ecc * ecc);

E[2] = boost::math::constants::pi<double>() / 2 * (1 - I) + 2 * I * std::atan(tmp);
E[3] = std::atan2(EQ[3] / tmp, EQ[4] / tmp);
E[3] = std::atan2(EQ[4] / tmp, EQ[3] / tmp);
E[4] = zita - I * E[3];
E[5] = EQ[5] - zita;
if (E[1] < 1) {
Expand Down
56 changes: 28 additions & 28 deletions src/core_functions/ic2eq.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,15 +35,15 @@ namespace kep_toolbox
{
/// From cartesian to modified equinoctial
/**
* Transforms cartesian coordinates (r,v) to modified equinoctial elements (p,h,k,p,q,L).
* Transforms cartesian coordinates (r,v) to modified equinoctial elements (p,f,g,h,k,L).
* Note that we use the true longitude L (i.e. f + om + Om).
*
* @param[in] r0 cartesian position vector.
* @param[in] v0 cartesian velocity vector.
* @param[in] mu gravitational parameter.
* @param[in] retrogade forces retrograde elements to be used
*
* @param[out] EQ the modified equinoctial elements (p,h,k,p,q,L).
* @param[out] EQ the modified equinoctial elements (p,f,g,h,k,L).
*
*/
template <class vettore3D, class vettore6D>
Expand All @@ -56,8 +56,8 @@ void ic2eq(const vettore3D &r0, const vettore3D &v0, const double &mu, vettore6D
I = 1;
}
// The equinoctial reference frame
vettore3D f = {0.0, 0.0, 0.0};
vettore3D g = {0.0, 0.0, 0.0};
vettore3D fv = {0.0, 0.0, 0.0};
vettore3D gv = {0.0, 0.0, 0.0};
vettore3D w = {0.0, 0.0, 0.0};
// The eccentricity vector
vettore3D e = {0.0, 0.0, 0.0};
Expand All @@ -74,55 +74,55 @@ void ic2eq(const vettore3D &r0, const vettore3D &v0, const double &mu, vettore6D
cross(w, r0, v0);
vers(w, w);

auto p = w[0] / (1 + I * w[2]);
auto q = -w[1] / (1 + I * w[2]);
auto den = p * p + q * q + 1;
f[0] = (1. - p * p + q * q) / den;
f[1] = (2. * p * q) / den;
f[2] = (-2. * I * p) / den;
auto k = w[0] / (1 + I * w[2]);
auto h = -w[1] / (1 + I * w[2]);
auto den = k * k + h * h + 1;
fv[0] = (1. - k * k + h * h) / den;
fv[1] = (2. * k * h) / den;
fv[2] = (-2. * I * k) / den;

g[0] = (2. * I * p * q) / den;
g[1] = (1. + p * p - q * q) * I / den;
g[2] = (2. * q) / den;
gv[0] = (2. * I * k * h) / den;
gv[1] = (1. + k * k - h * h) * I / den;
gv[2] = (2. * h) / den;

// 2 - We compute evett: the eccentricity vector
vettore3D Dum_Vec = {0.0, 0.0, 0.0};
cross(Dum_Vec, v0, ang);
for (auto i = 0u; i < 3; ++i) {
e[i] = Dum_Vec[i] / mu - r0[i] / R0;
}
auto h = dot(e, g);
auto k = dot(e, f);
auto g = dot(e, gv);
auto f = dot(e, fv);
auto ecc = norm(e);

// 3 - We compute the true longitude L
// This solution is certainly not the most elegant, but it works and will never be singular.

auto det1 = (g[1]*f[0]-f[1]*g[0]); // xy
auto det2 = (g[2]*f[0]-f[2]*g[0]); // xz
auto det3 = (g[2]*f[1]-f[2]*g[1]); // yz
auto det1 = (gv[1]*fv[0]-fv[1]*gv[0]); // xy
auto det2 = (gv[2]*fv[0]-fv[2]*gv[0]); // xz
auto det3 = (gv[2]*fv[1]-fv[2]*gv[1]); // yz
auto max = std::max({std::abs(det1), std::abs(det2), std::abs(det3)});

double X, Y;
if (std::abs(det1) == max) {
X = (g[1]*r0[0] - g[0] * r0[1]) / det1;
Y = (-f[1]*r0[0]+f[0]*r0[1]) / det1;
X = (gv[1]*r0[0] - gv[0] * r0[1]) / det1;
Y = (-fv[1]*r0[0]+fv[0]*r0[1]) / det1;
} else if (std::abs(det2) == max) {
X = (g[2]*r0[0] - g[0] * r0[2]) / det2;
Y = (-f[2]*r0[0]+f[0]*r0[2]) / det2;
X = (gv[2]*r0[0] - gv[0] * r0[2]) / det2;
Y = (-fv[2]*r0[0]+fv[0]*r0[2]) / det2;
} else {
X = (g[2]*r0[1] - g[1] * r0[2]) / det3;
Y = (-f[2]*r0[1]+f[1]*r0[2]) / det3;
X = (gv[2]*r0[1] - gv[1] * r0[2]) / det3;
Y = (-fv[2]*r0[1]+fv[1]*r0[2]) / det3;
}

auto L = std::atan2(Y/R0, X/R0);

// 5 - We assign the results
EQ[0] = a * (1. - ecc * ecc);
EQ[1] = h;
EQ[2] = k;
EQ[3] = p;
EQ[4] = q;
EQ[1] = f;
EQ[2] = g;
EQ[3] = h;
EQ[4] = k;
EQ[5] = L;
}
} // namespace kep_toolbox end
Expand Down
14 changes: 7 additions & 7 deletions src/core_functions/par2eq.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ namespace kep_toolbox
/// From osculating Keplerian to equinoctial
/**
* Transforms osculating Keplerian parameters (a,e,i,W,w,E) to
* modified equinoctial elements (p,h,k,p,q,L).
* modified equinoctial elements (p,f,g,h,k,L).
* Note that we use the eccentric anomaly E (or Gudermannian) and the true longitude L (i.e. f + om + Om)
*
* @param[in] E the osculating Keplerian parameters (a,e,i,W,w,E), anomalies in rad.
Expand All @@ -51,16 +51,16 @@ void par2eq(vettore6D &EQ, const vettore6D &E, const bool retrogade = false)
int I;
if (retrogade) {
I = -1;
EQ[3] = 1. / std::tan(E[2] / 2) * std::sin(E[3]);
EQ[4] = 1. / std::tan(E[2] / 2) * std::cos(E[3]);
EQ[3] = 1. / std::tan(E[2] / 2) * std::cos(E[3]);
EQ[4] = 1. / std::tan(E[2] / 2) * std::sin(E[3]);
} else {
I = 1;
EQ[3] = std::tan(E[2] / 2) * std::sin(E[3]);
EQ[4] = std::tan(E[2] / 2) * std::cos(E[3]);
EQ[3] = std::tan(E[2] / 2) * std::cos(E[3]);
EQ[4] = std::tan(E[2] / 2) * std::sin(E[3]);
}
EQ[0] = E[0] * (1 - E[1] * E[1]);
EQ[1] = E[1] * std::sin(E[4] + I * E[3]);
EQ[2] = E[1] * std::cos(E[4] + I * E[3]);
EQ[1] = E[1] * std::cos(E[4] + I * E[3]);
EQ[2] = E[1] * std::sin(E[4] + I * E[3]);
double f;
if (E[1] < 1) {
f = e2f(E[5], E[1]);
Expand Down

0 comments on commit 04638d5

Please sign in to comment.