Computing_Spacecraft_Position_and_Velocity_from_Orbital_Elements

Problem Statement

Required Knowledge

This problem can be solved with the following knowledge:

• Converting Classical Orbital Elements (CoE) to spacecraft's position and velocity vectors (rv2coe)
• Converting Satellite coordinate system (RSW) to Geocentric Equatorial System (IJK) (rsw2ijk)
• Solving Kepler's Equation with Newton's Iteration method

Solution Steps

1. Obtain the semi-major axis ($a$), eccentricity ($e$), and angular momentum ($h$) from the initial position and velocity vector with rv2coe.m

   Computing-Spacecraft-Position-and-Velocity-from-Orbital-Elements/rv2coe.m

   Lines 1 to 37 in 74a83d5

   	function [a, e_norm, i, omega, w, f, h_norm] = rv2coe(r,v,mu)
   	h = cross(r,v); % Angular momentum vector
   	h_norm = norm(h);
   	k = [0; 0; 1]; % Unit vector pointing to the normal vector of the plane
   	n = cross(k, h); % Nodal vector
   	e = (-cross(h,v)/mu) -(r/norm(r));
   	e_norm = norm(e); % eccentricity
   	a = (norm(h)^2)/(mu*(1-norm(e)^2)); % semi-major axis
   	i = acos(dot(k,h)/norm(h)); % inclination
   	%% Right ascension of the ascending node, Omega
   	if (n(2) < 0)
   	omega = 2*pi - acos(n(1)/norm(n));
   	else
   	omega = acos(n(1)/norm(n));
   	end
   	omega_deg = omega * 180/pi;
   	%% Argument of periaposis, w
   	if (dot(e,k)<0)
   	w = 2*pi - acos(dot(n,e)/(norm(n)*norm(e)));
   	else
   	w = acos(dot(n,e)/(norm(n)*norm(e)));
   	end
   	w_deg = w*180/pi;
   	%% True anomaly, f
   	if (dot(r,v)<0)
   	f = 2*pi - acos(dot(r,e)/(norm(r)*norm(e)));
   	else
   	f = acos(dot(r,e)/(norm(r)*norm(e)));
   	end
   	f_deg = f*180/pi;
   	end

2. For each time step, obtain the eccentric anomaly ($E$) from the given mean anomaly ($M$) with Newton's iteration. This has been coded in eccanomaly_newt.m

   Computing-Spacecraft-Position-and-Velocity-from-Orbital-Elements/eccanomaly_newt.m

   Lines 1 to 15 in b3d1e54

   	function E = eccanomaly_newt(E0,M,e)
   	% Calculate Eccentric Anomaly (E) from Newton's iteration
   	E = E0;
   	while (true)
   	% Newton's Iteration
   	Enew = E - ((E - e*sin(E) - M)/(1 - e*cos(E)));
   	% Converge condition
   	if abs(Enew - E) < 1e-8
   	break
   	end
   	E = Enew;
   	end
   	end

3. Calculate True Anomaly ($f$) from the following equation: $${f = 2\arctan\left(\sqrt{ 1+e \over 1-e} \times \tan{E \over 2}\right)}$$

4. Position and velocity vector in RSW coordinate frame can be obtained as follow:

$${ { {\vec{r}_{rsw}} = \left\lbrack \matrix{r \cr 0 \cr 0} \right\rbrack } ,{{\vec{v}_{rsw}} = \left\lbrack \matrix{{h \over a(1-e^2)} e\sin f \cr h/r \cr 0} \right\rbrack} }$$

5. Transform position and velocity vectors from RSW to IJK coordinate with directional cosine matrices. This is coded in rsw2ijk.m

Computing-Spacecraft-Position-and-Velocity-from-Orbital-Elements/rsw2ijk.m

Lines 1 to 21 in 366f8e2

	function [r, v, u] = rsw2ijk(r_rsw, v_rsw, omega, i, w, f)
	function mat1 = rot1(x)
	mat1 = [1 0 0;
	0 cos(x) sin(x);
	0 -sin(x) cos(x)];
	end
	function mat3 = rot3(x)
	mat3 = [cos(x) sin(x) 0;
	-sin(x) cos(x) 0;
	0 0 1];
	end
	u = f + w;
	rot313 = rot3(-omega) * rot1(-i) * rot3(-u);
	r = rot313 * r_rsw;
	v = rot313 * v_rsw;
	end

The main loop for finding spacecraft's position and velocity is shown below:

Computing-Spacecraft-Position-and-Velocity-from-Orbital-Elements/main.m

Lines 31 to 53 in 82d4dfc

	% Calculate the spacecraf's future position
	for j = 1:length(t)
	% Eccentric anomaly (rad)
	E0 = M(j) + e*sin(M(j)) + e*e*0.5*sin(2*M(j));
	E = eccanomaly_newt(E0, M(j), e);
	% True anomaly (rad)
	f = 2*atan(sqrt((1+e)/(1-e))*tan(E/2));
	% Get r and v vector in RSW coordinate (coe2rv)
	r_rsw_norm = a*(1-e*cos(E));
	r_rsw = [r_rsw_norm; 0; 0];
	v_rsw = [h * e * sin(f)/(a*(1 - e^2)); h/r_rsw_norm; 0];
	% Convert RSW to IJK coordinate (rsw2ijk)
	[r, v, u] = rsw2ijk(r_rsw, v_rsw, omega, i, w, f)
	% Extract x, y, z components
	rx(j) = r(1); ry(j) = r(2); rz(j) = r(3);
	vx(j) = v(1); vy(j) = v(2); vz(j) = v(3);
	end

Results

More time steps

Animated Simulation

Molniya_Orbit_Compressed.mp4

Answer

References

[1] H. Curtis, Orbital mechanics for engineering students, Butterworth-Heinemann, 2013.

[2] J. E. Prussing, B. A. Conway, Orbital mechanics, Oxford University Press, 2012.