Convert Position/Velocity Vectors to Keplerian Elements in MATLAB

vec2orbElem(rs, vs, mus) converts the position vector (rs) and velocity vector (vs) of an object into Keplerian orbital elements using the gravitational parameter (mus).

Input:

- rs: 3n x 1 stacked initial position vector:

[r1(1); r1(2); r1(3); r2(1); r2(2); r2(3); ...; rn(1); rn(2); rn(3)]

or 3 x n position matrix.

- vs: 3n x 1 stacked initial velocity vector or 3 x n matrix.

- mus: gravitational parameter (G*m_i) where G is the gravitational constant and m_i is the mass of the i-th object. If all vectors represent the same object, mus can be a scalar.

Output:

- a: semi-major axis

- e: eccentricity

- E: eccentric anomaly

- i: inclination

- omega: argument of periapsis

- Omega: longitude of ascending node

- P: orbital period

- tau: time of periapsis passage

- A, B: direction matrices (see Vinti, 1998)

All units must be consistent; for example, if position is in AU, time should be in days.