Small Satellites

Home » Posts tagged 'Two Body Problem'

Tag Archives: Two Body Problem

Relative Motion of Satellites, Numerical Simulation

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. We numerically solve fundamental equation of relative two-body motion to obtain the trajectory of spacecraft B relative to spacecraft A and the distance between two satellites. 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

clc;
clear all;
% Initial State vectors of Satellite A and B
RA0 = [-266.77,  3865.8, 5426.2];     % [km]
VA0 = [-6.4836, -3.6198, 2.4156];     % [km/s]
RB0 = [-5890.7, -2979.8, 1792.2];     % [km]
VB0 = [0.93583, -5.2403, -5.5009];    % [km/s]
mu  = 398600;            % Earth’s gravitational parameter [km^3/s^2]
t   = 0;                 % initial time
dt  = 15;                % Simulation time step [s]
dT  = 4*24*3600;         % Simulation interval  [s]
% Using fourth-order Runge–Kutta method to solve fundamental equation
% of relative two-body motion
F_r = @(R) -mu/(norm(R)^3)*R;
VA = VA0; RA  = RA0;
VB = VB0; RB  = RB0;
ind = 1;
while (t <= dT)
    % Relative position
    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(ind,:) = QXx*Rr';

    % A Satellite
    k_1 = dt*F_r(RA);  k_2 = dt*F_r(RA+0.5*k_1);
    k_3 = dt*F_r(RA+0.5*k_2);  k_4 = dt*F_r(RA+k_3);
    VA  = VA + (1/6)*(k_1+2*k_2+2*k_3+k_4);
    RA  = RA + VA*dt;
    % B Satellite
    k_1 = dt*F_r(RB);  k_2 = dt*F_r(RB+0.5*k_1);
    k_3 = dt*F_r(RB+0.5*k_2);     k_4 = dt*F_r(RB+k_3);
    VB = VB + (1/6)*(k_1+2*k_2+2*k_3+k_4);
    RB  = RB + VB*dt;

    R_A(ind,:)  = RA;
    R_B(ind,:)  = RB;
    time(ind)   = t;
    t   = t+dt;
    ind = ind+1;
end
r_BA = (R_BA(:,1).^2+R_BA(:,2).^2+R_BA(:,3).^2).^0.5;
close all;
figure(1);
hold on;
plot3(R_A(:,1),R_A(:,2),R_A(:,3),'r');
plot3(R_B(:,1),R_B(:,2),R_B(:,3),'y');
title('Satellites orbits around earth');
legend('Satellite A','Satellites B');
xlabel('km');ylabel('km');
% Plotting Earth
load('topo.mat','topo','topomap1');
colormap(topomap1);
% Create the surface.
radius_earth=6378;
[x,y,z] = sphere(50);
x =radius_earth*x;
y =radius_earth*y;
z =radius_earth*z;
props.AmbientStrength = 0.1;
props.DiffuseStrength = 1;
props.SpecularColorReflectance = .5;
props.SpecularExponent = 20;
props.SpecularStrength = 1;
props.FaceColor= 'texture';
props.EdgeColor = 'none';
props.FaceLighting = 'phong';
props.Cdata = topo;
surface(x,y,z,props);
hold off;

figure(2);
plot3(R_BA(:,1),R_BA(:,2),R_BA(:,3),'k');
title('The trajectory of spacecraft B relative to spacecraft A');
xlabel('km');ylabel('km');zlabel('km');

figure(3);
plot(time/3600,r_BA);
title('Distance between two satellites');
xlabel('hour');ylabel('km')

min_r = min(r_BA);
max_r = max(r_BA);
fprintf('Max distance between two satellites %6.4f km \n',max_r);
fprintf('Min distance between two satellites %6.4f km \n',min_r);
Max distance between two satellites 13850.3054 km 
Min distance between two satellites 262.0271 km

Relative_motion_sim_01 Relative_motion_sim_02 Relative_motion_sim_03

Advertisements

Two Body Problem Numerical Solution,Satellite – Earth, R & V after dT

The initial position and velocity of an earth orbiting satellite in earth centered inertial frame is known. In this example we will numerically solve fundamental equation of relative two-body motion to find the distance of the satellite from the center of the earth and its speed after 24 hours.

clc;
clear all;
R0 = [6750 0 0];          %[km]
V0 = [0 10.5 0];          %[km/s]
mu = 398600;            % Earth’s gravitational parameter [km^3/s^2]
t = 0;                  % initial time
dt = 10;                % time step [s]
dT = 24*3600;           % Time interval [s]
% Using fourth-order Runge–Kutta method to solve fundamental equation
% of relative two-body motion
F_r = @(R) -mu/(norm(R)^3)*R;
Rd = V0; R  = R0;
i = 1;
while (t <= dT)
    Rv(i,:) = R;
    tv(i) =t;
    k_1 = dt*F_r(R);
    k_2 = dt*F_r(R+0.5*k_1);
    k_3 = dt*F_r(R+0.5*k_2);
    k_4 = dt*F_r(R+k_3);
    Rd  = Rd + (1/6)*(k_1+2*k_2+2*k_3+k_4);
    Vv(i,:) = Rd;
    R   = R + Rd*dt;
    t   = t+dt;
    i = i+1;
end
rn = (Rv(:,1).^2+Rv(:,2).^2+Rv(:,3).^2).^0.5;   % Radius Vector
vn = (Vv(:,1).^2+Vv(:,2).^2+Vv(:,3).^2).^0.5;   % Speed Vector
fprintf('Distance from Earth center = %4.2f [km] \n',norm(R));
fprintf('Satellite speed= %4.4f [km/s] \n',norm(Rd));
Distance from Earth center = 77061.78 [km] 
Satellite speed= 1.5791 [km/s]

Plots

figure(1);
hold on;grid on;
plot(tv/3600,rn);
ylabel('Altitude [km]');
xlabel('Time [hour]');
title('Distance variation of the satellite from the center of the Earth');
figure(2);
hold on;grid on;
plot(tv/3600,vn);
ylabel('Speed [km/s]');
xlabel('Time [hour]');
title('Satellite speed variation');

Satellite_two_body_problem_dT_01 Satellite_two_body_problem_dT_02

Two Body Problem Numerical Solution,Satellite – Earth

Given initial position and velocity of an earth orbiting satellite at a given instant.In this example we will numerically solve fundamental equation of relative two-body motion to find the maximum altitude(apogee altitude) reached by the satellite.

clc;
clear all;
R0 = [3102 5369 2625];          %[km]
V0 = [-6.426 0.7735 5.943];     %[km/s]
r = norm(R0);           % Initial radius  [km]
v = norm(V0);           % Initial speed   [km/s]
mu = 398600;            % Earth’s gravitational parameter [km^3/s^2]
a = mu/(2*mu/r - v^2);  % Semimajor Axis [km]
T = 2*pi*a^1.5/mu^0.5;  % Orbital period [s]
R_earth = 6378;         % Earth radius [km]
dt = T/1000;            % time step [s]
t = 0;                  % initial time
% Using fourth-order Runge–Kutta method to solve fundamental equation
% of relative two-body motion
F_r = @(R) -mu/(norm(R)^3)*R;
Rd = V0; R  = R0;
i = 1;
while (t <= T)
    Rv(i,:) = R;
    tv(i) =t;
    k_1 = dt*F_r(R);
    k_2 = dt*F_r(R+0.5*k_1);
    k_3 = dt*F_r(R+0.5*k_2);
    k_4 = dt*F_r(R+k_3);
    Rd  = Rd + (1/6)*(k_1+2*k_2+2*k_3+k_4);
    R   = R + Rd*dt;
    t   = t+dt;
    i = i+1;
end
rn = (Rv(:,1).^2+Rv(:,2).^2+Rv(:,3).^2).^0.5;
Alt = rn - R_earth; % Altitude Vector
Alt_max  = max(Alt); % [km]
fprintf('Maximum Altitude = %4.2f [km] \n',Alt_max);
Maximum Altitude = 6249.10 [km]

Ploting the altitude variation during one orbit

figure(1);
hold on;grid on;
plot(tv/3600,Alt);
ylabel('Altitude [km]');
xlabel('Time [hour]');
title('Satellite altitude variation during one orbit');

Earth_Satellite_two_body_problem_01

 

%d bloggers like this: