From b3e5819525ac73f397705c33abacd7e3a736c686 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Thu, 22 Nov 2018 16:04:59 +0100 Subject: [PATCH 1/2] Equinoctial elements names changed to follow Walker original paper. Read Tle changed to py3 syntax --- pykep/util/__init__.py | 4 +-- src/core_functions/eq2ic.h | 40 +++++++++++++------------- src/core_functions/eq2par.h | 8 +++--- src/core_functions/ic2eq.h | 56 ++++++++++++++++++------------------- src/core_functions/par2eq.h | 14 +++++----- 5 files changed, 61 insertions(+), 61 deletions(-) diff --git a/pykep/util/__init__.py b/pykep/util/__init__.py index cd950e3d..bbbee4db 100644 --- a/pykep/util/__init__.py +++ b/pykep/util/__init__.py @@ -69,9 +69,9 @@ def read_tle(tle_file, verbose=False, with_name=True): with open(tle_file, 'r') as f: for line in f: if with_name: - line1 = f.next()[:69] + line1 = f.readline()[:69] else: line1 = line[:69] - line2 = f.next()[:69] + line2 = f.readline()[:69] planet_list.append(tle(line1, line2)) return planet_list diff --git a/src/core_functions/eq2ic.h b/src/core_functions/eq2ic.h index 69615728..280a7f65 100644 --- a/src/core_functions/eq2ic.h +++ b/src/core_functions/eq2ic.h @@ -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 @@ -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"); diff --git a/src/core_functions/eq2par.h b/src/core_functions/eq2par.h index 8660050c..73ee5c5a 100644 --- a/src/core_functions/eq2par.h +++ b/src/core_functions/eq2par.h @@ -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). @@ -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() / 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) { diff --git a/src/core_functions/ic2eq.h b/src/core_functions/ic2eq.h index 8034d61e..7fb3af97 100644 --- a/src/core_functions/ic2eq.h +++ b/src/core_functions/ic2eq.h @@ -35,7 +35,7 @@ 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. @@ -43,7 +43,7 @@ namespace kep_toolbox * @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 @@ -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}; @@ -74,16 +74,16 @@ 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}; @@ -91,38 +91,38 @@ void ic2eq(const vettore3D &r0, const vettore3D &v0, const double &mu, vettore6D 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 diff --git a/src/core_functions/par2eq.h b/src/core_functions/par2eq.h index 98bf3a61..96234c85 100644 --- a/src/core_functions/par2eq.h +++ b/src/core_functions/par2eq.h @@ -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. @@ -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]); From 590d71643c1ed42433089863ba31cf300f2e41ae Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Thu, 22 Nov 2018 16:06:55 +0100 Subject: [PATCH 2/2] Python docs changed to new equinoctial elements convention --- pykep/core/core.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pykep/core/core.cpp b/pykep/core/core.cpp index c26adcba..596d2768 100644 --- a/pykep/core/core.cpp +++ b/pykep/core/core.cpp @@ -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" @@ -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" @@ -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)",