% This example shows how to use the state vectors of spacecraft A and B
% to find the position, velocity and acceleration of B relative
% to A in the LVLH frame attached to A. For more details about the theory and
% algorithm look at Chapter 7 of  H. D. Curtis, Orbital Mechanics for
% Engineering Students,Second Edition, Elsevier 2010

Input

State vectors of Satellite A and B

RA = [-266.77,  3865.8, 5426.2];     % [km]
VA = [-6.4836, -3.6198, 2.4156];     % [km/s]
RB = [-5890.7, -2979.8, 1792.2];     % [km]
VB = [0.93583, -5.2403, -5.5009];    % [km/s]

% Earth gravitational parameter
mu = 398600;                        % [km^3/s^2]

Algorithm

hA = cross(RA, VA);                 % Angular momentum of A
% Unit vectors i, j,k of the co-moving frame
i = RA/norm(RA);  k = hA/norm(hA); j = cross(k,i);

% Transformation matrix Qxx:
QXx   = [i; j; k];
Om    = hA/norm(RA)^2;
Om_dt = -2*VA*RA'/norm(RA)^2.*Om;

% Accelerations of A and B,inertial frame
aA = -mu*RA/norm(RA)^3;
aB = -mu*RB/norm(RB)^3;

% Relative position,inertial frame
Rr = RB - RA;
% Relative position,LVLH frame attached to A
R_BA = QXx*Rr';

% Relative velocity,inertial frame
Vr = VB - VA - cross(Om,Rr);
% Relative velocity,LVLH frame attached to A
V_BA = QXx*Vr';

% Relative acceleration, inertial frame
ar = aB - aA - cross(Om_dt,Rr) - cross(Om,cross(Om,Rr))- 2*cross(Om,Vr);
% Relative acceleration,LVLH frame attached to A
a_BA = QXx*ar';

fprintf('Position of B relative to A in LVLH frame attached to A \n');
fprintf('R_BA = [%4.2f %4.2f %4.2f] km \n\n', R_BA);

fprintf('Velocity of B relative to A in LVLH frame attached to A \n');
fprintf('V_BA = [%6.4f %6.4f %6.4f] km/s \n\n', V_BA);

fprintf('Acceleration of B relative to A in LVLH frame attached to A \n');
fprintf('a_BA = [%8.8f %8.8f %8.8f] km/s^2 \n', a_BA);

Position of B relative to A in LVLH frame attached to A 
R_BA = [-6701.22 6828.28 -406.24] km 

Velocity of B relative to A in LVLH frame attached to A 
V_BA = [0.3168 0.1120 1.2470] km/s 

Acceleration of B relative to A in LVLH frame attached to A 
a_BA = [-0.00022213 -0.00018083 0.00050590] km/s^2