% Richard Rieber
% September 26, 2006
% rrieber@gmail.com
%
% function [a,ecc,inc,Omega,w,nu,w_true,u_true,lambda_true] = elorb(R,V)
%
% function [a,ecc,inc,Omega,w,nu] = elorb(R,V,U)
%
% function [...] = elorb(R,V,U,tol)
%
% Purpose: This function converts ECI Cartesian coordinates (R & V)
% to Kepler orbital elements (a, e, i, O, w, nu).
% This function is derived from algorithm 9 on pg. 120 of
% David A. Vallado's Fundamentals of Astrodynamics and
% Applications (2nd Edition)
%
% Inputs: R: n x 3 radius vector from the center of the Earth to the
% spacecraft (ECI coordinates) in km. n rows should be used
% for multiple inputs
% V: A n x 3 velocity vector of the space craft in km/s.
% n rows should be used for multiple inputs
% U: Gravitational constant of body being orbited (km^3/s^2). Default is Earth
% at 398600.4415 km^3/s^2. [OPTIONAL]
% tol: A tolerance value for checking equality. Default is 10^-8.
% [OPTIONAL]
%
% Outputs: The user has two options for outputs. If 6 outputs are asked for, the
% normal 6 Kepler elements will be returned (a,ecc,inc,Omega,w,nu).
% If 9 outputs are asked for, w_true, u_true, and lambda_true will also
% be returned. These parameters are important for near equatorial elliptical
% orbits, inclined near circular orbits, and near equatorial, near circular
% orbits
%
% a: Semi-major axis of orbit in km
% ecc: Eccentricity of orbit
% inc: Inclination of orbit in radians
% Omega: Right ascension of the ascending node in radians
% w: Argument of perigee in radians
% nu: True anomaly in radians
% w_true: True longitude of periapse in radians if equatorial
% elliptical orbit [OPTIONAL]
% u_true: Argument of Latitude in radians if orbiit is circular
% inclined [OPTIONAL]
% lambda_true: True longitude in radians if orbit is equatorial
% circular [OPTIONAL]
function [a,ecc,inc,Omega,w,nu,w_true,u_true,lambda_true] = elorb(R,V,U)
if nargin < 2 || nargin > 4
error('Incorrect number of inputs: see help elorb')
elseif nargin == 2
U = 398600.4415; %km^3/s^2 Gravitational Constant of Earth
tol = 10^-8; %A tolerance for checking equality
elseif nargin == 3
tol = 10^-8; %A tolerance for checking equality
end
[x,y] = size(R);
if y ~= 3
error('Radius vector is not the right size: see help elorb')
end
if size(R) ~= size(V)
error('Velocity and Radius vectors are not the same size. Check inputs')
end
ecc = zeros(1,x);
a = zeros(1,x);
w = zeros(1,x);
nu = zeros(1,x);
Omega = zeros(1,x);
inc = zeros(1,x);
w_true = zeros(1,x);
u_true = zeros(1,x);
lambda_true = zeros(1,x);
for j = 1:x
h = cross(R(j,:),V(j,:)); %Specfic angular momentum vector
N = cross([0,0,1],h);
%Eccenticity vector
e_vec = ((norm(V(j,:))^2 - U/norm(R(j,:)))*R(j,:) - dot(R(j,:),V(j,:))*V(j,:))/U;
ecc(j) = norm(e_vec); %Magnitude of eccentricty vector
zeta = norm(V(j,:))^2/2 - U/norm(R(j,:)); %Specific mechanical energy of orbit
if (1 - abs(ecc(j))) > tol %Checking to see if orbit is parabolic
a(j) = -U/(2*zeta); %Semi-major axis
P = a(j)*(1-ecc(j)^2); %semi-parameter
else
P = norm(h)^2/U; %Semi-parameter
a(j) = inf;
end
inc(j) = acos(h(3)/norm(h)); %inclination of orbit in radians
Omega(j) = acos(N(1)/norm(N)); %Right ascension of ascending node in radians
if N(2) < 0 %Checking for quadrant
Omega(j) = 2*pi - Omega;
end
w(j) = acos(dot(N,e_vec)/(norm(N)*ecc(j))); %Argument of perigee in radians
if e_vec(3) < 0 %Checking for quadrant
w(j) = 2*pi - w(j);
end
nu(j) = acos(dot(e_vec,R(j,:))/(ecc(j)*norm(R(j,:)))); %True anomaly in radians
if dot(R(j,:),V(j,:)) < 0 %Checking for quadrant
nu(j) = 2*pi - nu(j);
end
%%%%%%%%%%%% Special Cases %%%%%%%%%%%%%%%
% Elliptical Equatorial
w_true(j) = acos(e_vec(1)/ecc(j)); %True longitude of periapse
if e_vec(2) < 0 %Checking for quadrant
w_true(j) = 2*pi - w_true(j);
end
% Circular Inclined
u_true(j) = acos(dot(N,R(j,:))/(norm(N)*norm(R(j,:)))); %Argument of Latitude
if R(j,3) < 0 %Checking for quadrant
u_true(j) = 2*pi - u_true(j);
end
% Circular Equatorial
lambda_true(j) = acos(R(j,1)/norm(R(j,:))); %True Longitude
if R(j,2) < 0 %Checking for quadrant
lambda_true(j) = 2*pi - lambda_true(j);
end
end
- 1
- 2
- 3
- 4
前往页