Source code for orbitdeterminator.util.state_kep

'''
Takes a state vector (x, y, z, vx, vy, vz) where v is the velocity of the satellite
and transforms it into a set of keplerian elements (a, e, i, ω, Ω, v)

'''

import numpy as np
import math

[docs]def state_kep(r, v): ''' Converts state vector to orbital elements. Args: r (numpy array): position vector v (numpy array): velocity vector Returns: numpy array: array of the computed keplerian elements kep(0): semimajor axis (kilometers) kep(1): orbital eccentricity (non-dimensional) (0 <= eccentricity < 1) kep(2): orbital inclination (degrees) kep(3): argument of perigee (degress) kep(4): right ascension of ascending node (degrees) kep(5): true anomaly (degrees) ''' mu = 398600.4405 mag_r = np.sqrt(r.dot(r)) mag_v = np.sqrt(v.dot(v)) h = np.cross(r, v) mag_h = np.sqrt(h.dot(h)) e = ((np.cross(v, h)) / mu) - (r / mag_r) mag_e = np.sqrt(e.dot(e)) n = np.array([-h[1], h[0], 0]) mag_n = np.sqrt(n.dot(n)) true_anom = math.acos(np.clip(np.dot(e,r)/(mag_r * mag_e), -1, 1)) if np.dot(r, v) < 0: true_anom = 2 * math.pi - true_anom true_anom = math.degrees(true_anom) i = math.acos(np.clip(h[2] / mag_h, -1, 1)) i = math.degrees(i) ecc = mag_e raan = math.acos(np.clip(n[0] / mag_n, -1, 1)) if n[1] < 0: raan = 2 * math.pi - raan raan = math.degrees(raan) per = math.acos(np.clip(np.dot(n, e) / (mag_n * mag_e), -1, 1)) if e[2] < 0: per = 2 * math.pi - per per = math.degrees(per) a = 1 / ((2 / mag_r) - (mag_v**2 / mu)) if i >= 360.0: i = i - 360 if raan >= 360.0: raan = raan - 360 if per >= 360.0: per = per - 360 kep = np.zeros(6) kep[0] = a kep[1] = ecc kep[2] = i kep[3] = per kep[4] = raan kep[5] = true_anom return kep
if __name__ == "__main__": r = np.array([5.0756899358316559e+03, -4.5590381308371752e+03, 1.9322228177731663e+03]) v = np.array([1.3360847905126974e+00, -1.5698574946888049e+00, -7.2117328822023676e+00]) kep = state_kep(r, v) print(kep)