Small Satellites

Home » 2013 » January

Monthly Archives: January 2013

Rocket Dynamics, VEGA Rocket sub-orbital ascent profile(1st stage P80)

Vega is an expendable launch system in use by Arianespace. Its jointly developed by the Italian Space Agency and the European Space Agency. First time it was launched from Guiana Space Center on 13 February 2012. In this example we show how to simulate rocket launch trajectory. VEGA rocket data from VEGA Users Manual/ The full sub-orbital ascent profile simulation coming soon !!

clc;
clear all;
global m0 g0 T A Cd rh0 H0 Re hgr_turn tf md
% Launch Site: Guiana Space Center
Alt = 1;              %[m] Alt above sea level

% VEGA Rocket
m_stage_gross = [95796, 25751,10948];% 1st, 2nd,3d
% First stage(Solid Fuel)
m_prop = 88365;       % [kg] Propellant mass
Isp    = 280 ;        % [s]  Specific impulse
d      = 3;           % [m]  Diameter
g0     = 9.81;        % [m/s^2] Constant at its sea-level value
m0  = 137000;         % [kg] Initial mass
A   = pi*d^2/4;       % [m^2]Frontal area
Cd  = 0.5 ;             % Drag coefficient,assumed to have the constant value
rh0 = 1.225;          % [kg/m^3]
H0 = 7500;            % [m] Density scale height
Re = 6378e3;          % [m] Earth's radius
hgr_turn = 200;       % [m] Rocket starts the gravity turn when h = hgr_turn
tburn = 106.8;        % [s] Fuell burn time, first stage
md = (m_prop)/tburn;  % [kg/s]Propellant mass flow rate
T   = md*(Isp*g0);    % [N] Thrust (mean)
mf = m0 - m_prop;     % [kg] Final mass of the rocket(first stage is empty)
t0 = 0;               % Rocket launch time
tf = t0 + tburn;      % The time when propellant is completely burned
%and the thrust goes to zero
t_range     = [t0,tf];  % Integration interval

% Launch initial conditions:
gamma0 = 89.5/180*pi;       % Initial flight path angle
v0 = 0;   % Velocity (m/s)  % Earth's Rotation considered in eq of motion.
x0 = 0;   % Downrange distance [km]
h0 = Alt; % Launch site altitude [km]
vD0 = 0;  % Loss due to drag (Velocity)[m/s]
vG0 = 0;  % Loss due to gravity (Velocity)[m/s]
state0   = [v0, gamma0, x0, h0, vD0, vG0];
% Solve initial value problem for ordinary differential equations
[t,state] = ode45(@RocketDynEq,t_range,state0) ;
v     = state(:,1)/1000;      % Velocity [km/s]
gamma = state(:,2)*180/pi;    % Flight path angle  [deg]
x     = state(:,3)/1000;      % Downrange distance [km]
h     = state(:,4)/1000;      % Altitude[km]
vD    = -state(:,5)/1000;     % Loss due to drag (Velocity)[m/s]
vG    = -state(:,6)/1000;     % Loss due to gravity (Velocity)[m/s]
plot(t,h,'b');
hold on;
grid on;
plot(t,h,'.b');
title('Rocket Dynamics,VEGA Rocket sub-orbital ascent profile(1st stage P80)');
xlabel('time[s]');
ylabel('Altitude[km]');
 text(80,5,'smallsats.org','Color',[0 0 1], 'VerticalAlignment','middle',...
	'HorizontalAlignment','left','FontSize',14 );

% VEGA Rocket: First Stage P80
fprintf('\n VEGA Rocket: First Stage P80\n')
fprintf('\n Propellant mass           = %4.2f [kg]',m_prop)
fprintf('\n Gross mass                = %4.2f [kg]',m_stage_gross(1))
fprintf('\n Isp                       = %4.2f [s]',Isp)
fprintf('\n Thrust(mean)              = %4.2f [kN]',T/1000)
fprintf('\n Initial flight path angle = %4.2f [deg]',gamma0*180/pi)
fprintf('\n Final speed               = %4.2f [km/s]',v(end))
fprintf('\n Final flight path angle   = %4.2f [deg]',gamma(end))
fprintf('\n Altitude                  = %4.2f [km]',h(end))
fprintf('\n Downrange distance        = %4.2f [km]',x(end))
fprintf('\n Drag loss                 = %4.2f [km/s]',vD(end))
fprintf('\n Gravity loss              = %4.2f [km/s]',vG(end))
fprintf('\n');
 VEGA Rocket: First Stage P80

 Propellant mass           = 88365.00 [kg]
 Gross mass                = 95796.00 [kg]
 Isp                       = 280.00 [s]
 Thrust(mean)              = 2272.67 [kN]
 Initial flight path angle = 89.50 [deg]
 Final speed               = 1.76 [km/s]
 Final flight path angle   = 80.11 [deg]
 Altitude                  = 71.51 [km]
 Downrange distance        = 54.55 [km]
 Drag loss                 = 0.05 [km/s]
 Gravity loss              = 1.04 [km/s]
VEGA Rocket sub-orbital ascent profile(1st stage P8

VEGA Rocket sub-orbital ascent profile(1st stage P8

RocketDynEq.m
function dfdt = RocketDynEq(t,y)
global m0 g0 T A Cd rh0 H0 Re hgr_turn md
v  =  y(1);     % Velocity
gm =  y(2);     % Flight path angle
x  =  y(3);     % Downrange distance
h  =  y(4);     % Altitude
vD =  y(5);     % Velocity loss due to drag
vG =  y(6);     % Velocity loss due to gravity
% Equations of motion of a gravity turn trajectory
      m = m0 - md*t;  % Vehicle mass
% else
%     m = mf;          % Burnout mass
%     T = 0;           % No more thrust is generated
% end
g  = g0/(1 + h/Re)^2;          % Gravitational variation with altitude
rh = rh0*exp(-h/H0);            % Atmospheric density exponential model
D = 1/2 * rh*v^2 * A * Cd;      % Drag force

% Rocket starts the gravity turn when h = hgr_turn
if h <= hgr_turn % Vertical flight
    dv_dt  = T/m - D/m - g;
    dgm_dt = 0;
    dx_dt  = 0;
    dh_dt  = v;
    dvG_dt = -g;
else                                  % Gravity turn
    dv_dt  = T/m - D/m - g*sin(gm);
    dgm_dt = -1/v*(g - v^2/(Re + h))*cos(gm);
    dx_dt  = Re/(Re + h)*v*cos(gm) + 463*sin(gm);
    dh_dt  = v*sin(gm) + 463*cos(gm); % Adding earth's rotation speed
    dvG_dt = -g*sin(gm);              % Gravity loss rate [m/s^2]
end
    dvD_dt = -D/m;           % Drag loss rate [m/s^2]
dfdt = [ dv_dt,dgm_dt, dx_dt,dh_dt, dvD_dt, dvG_dt]';
return

Advertisements

Gibbs method of preliminary orbit determination

This example will show how to use 3 earth-centered position vectors of a satellite R1,R2,R3 at successive times to determine satellite orbit.

clear all; clc;close all;
mu = 398600;               %[km^3/s^2] Earth’s gravitational parameter

R1 = [5887 -3520 -1204];   %[km]
R2 = [5572 -3457 -2376];   %[km]
R3 = [5088 -3289 -3480];   %[km]

r1 = norm(R1); r2 = norm(R2); r3 = norm(R3);
% Vector cross products
R12 = cross(R1,R2);
R23 = cross(R2,R3);
R31 = cross(R3,R1);

Nv  = r1*R23 + r2*R31 + r3*R12;
Dv  = R23 + R31 + R12;
Sv  = (r2-r3)*R1 + (r3-r1)*R2 + (r1-r2)*R3;

N   = norm(Nv);
D   = norm(Dv);
% Velocity vector
V1 = (mu/(N*D))^0.5*(cross(Dv,R1)/r1 + Sv);
[ h i RAAN e omega theta ] = State2Orbital(R1,V1);
OE = [h i RAAN e omega theta];

fprintf('h [km^2/s]    i [deg]     RAAN [deg]  e    omega[deg]   theta [deg] \n');
fprintf('%4.2f     %4.2f       %4.2f     %4.4f    %4.2f    %4.2f  \n',OE);
% Check, Only true anomaly
V2 = (mu/(N*D))^0.5*(cross(Dv,R2)/r2 + Sv);
[ h i RAAN e omega theta ] = State2Orbital(R2,V2);
OE = [h i RAAN e omega theta];

fprintf('h [km^2/s]    i [deg]     RAAN [deg]  e    omega[deg]   theta [deg] \n');
fprintf('%4.2f     %4.2f       %4.2f     %4.4f    %4.2f    %4.2f  \n',OE);
h [km^2/s]    i [deg]     RAAN [deg]  e    omega[deg]   theta [deg] 
52948.88     95.03       150.01     0.0127    151.69    38.30  
h [km^2/s]    i [deg]     RAAN [deg]  e    omega[deg]   theta [deg] 
52948.88     95.01       150.00     0.0127    151.69    48.31
EarthTopographicMap(1,820,420);
[Rv, alfa,delta ] = J2PropagR( h, i, RAAN, e,omega,theta) ;
scatter(alfa,delta,'.r');
text(270,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',...
	'HorizontalAlignment','left','FontSize',14 );
title('Satellite ground track');
% Earth 3D Plot
Earth3DPlot(2);
scatter3(Rv(:,1),Rv(:,2),Rv(:,3),'.r');
% Plot position vectors
% Direction Cosines
l = R1/r1; m = R2/r2; n = R3/r3;
lv =(0:1:r1)'*l; mv =(0:1:r2)'*m;nv =(0:1:r3)'*n;
figure(3);
scatter3(Rv(:,1),Rv(:,2),Rv(:,3),'.k');
hold on;
scatter3(lv(:,1),lv(:,2),lv(:,3),'.r');
scatter3(mv(:,1),mv(:,2),mv(:,3),'.g');
scatter3(nv(:,1),nv(:,2),nv(:,3),'.b');
legend('Orbit','R1','R2','R3');
title('Satellite Orbit and R1,R2,R3 successive position vectors');
zoom(2)
% Equations from the book
% Orbital Mechanics for Engineering Students, 2nd Edition, Aerospace EngineeringGibbsMethod_01

 

Universal anomaly, Lagrange f and g coefficients

Contents

Universal anomaly, Lagrange f and g coefficients

Given the initial position R0 and velocity V0 vector of a satelite at time t0 in earth centered earth fixed reference frame. In this example we show how to calculate position and velocity vectors of the satellite dt minutes later using Lagrange f and g coefficients in terms of the universal anomaly

clc; clear all;
R0 = [7200, -13200 0];     %[km]   Initial position
V0 = [3.500, 2.500 1.2];   %[km/s] Initial velocity
dt = 36000;                %[s] Time interval

mu = 398600;               %[km^3/s^2] Earth’s gravitational parameter
r0  = norm(R0);
v0  = norm(V0);
vr0 = R0*V0'/r0;
alfa   = 2/r0 - v0^2/mu;

if     (alfa > 0 )
    fprintf('The trajectory is an ellipse\n');
elseif (alfa < 0)
    fprintf('The trajectory is an hyperbola\n');
elseif (alfa == 0)
    fprintf('The trajectory is an parabola\n');
end

X0 = mu^0.5*abs(alfa)*dt;  %[km^0.5]Initial estimate of X0
Xi = X0;
tol = 1E-10;              % Tolerance
while(1)
    zi = alfa*Xi^2;
    [ Cz,Sz] = Stumpff( zi );
    fX  = r0*vr0/(mu)^0.5*Xi^2*Cz + (1 - alfa*r0)*Xi^3*Sz + r0*Xi -(mu)^0.5*dt;
    fdX = r0*vr0/(mu)^0.5*Xi*(1 - alfa*Xi^2*Sz) + (1 - alfa*r0)*Xi^2*Cz + r0;
    eps = fX/fdX;
    Xi = Xi - eps;
    if(abs(eps) < tol )
        break
    end
end
fprintf('Universal anomaly X = %4.3f [km^0.5] \n\n',Xi)
% Lagrange f and g coefficients in terms of the universal anomaly

f  =  1 - Xi^2/r0*Cz;
g  = dt - 1/mu^0.5*Xi^3*Sz;
R  = f*R0 + g*V0;
r = norm(R);
df = mu^0.5/(r*r0)*(alfa*Xi^3*Sz - Xi);
dg = 1 - Xi^2/r*Cz;
V = df*R0 + dg*V0;
fprintf('Position after %4.2f [min] R = %4.2f*i + %4.2f*j + %4.2f*k[km] \n',dt/60,R);
fprintf('Velocity after %4.2f [min] V = %4.3f*i + %4.3f*j + %4.2f*k[km/s] \n',dt/60,V);
The trajectory is an ellipse
Universal anomaly X = 1922.210 [km^0.5] 

Position after 600.00 [min] R = -6781.27*i + -11870.72*j + -3270.69*k[km] 
Velocity after 600.00 [min] V = 3.488*i + -3.362*j + 0.41*k[km/s]

Ploting trajectory

t      = 0;
step   = 100;     %[s] Time step
ind    = 1;
while (t < dt)
    [R V] = UniversalLagrange(R0, V0, t);
    Ri(ind,:) = R;
    [Long(ind,:), Lat(ind,:)] = R2RA_Dec(R);
    ind = ind + 1;
    t = t + step;
end
% Earth 3D Plot
Earth3DPlot(1);
scatter3(Ri(:,1),Ri(:,2),Ri(:,3),'.r');
hold off;
% Earth topographic map
EarthTopographicMap(2,820,420);
scatter(Long ,Lat,'.r');
text(280,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
 title('Satellite ground track');
 hold off;
%

Lagrange Universal Anomaly_02

Lagrange Universal Anomaly_01

Hyperbolic trajectory

clear Ri Long Lat
R0 = [7200, -6200 0];      %[km]   Initial position
V0 = [5.500, 7.500 1.2];   %[km/s] Initial velocity
dt = 12000;                %[s] Time interval

% Ploting trajectory
t      = 0;
step   = 100;     %[s] Time step
ind    = 1;
while (t < dt)
    [R V] = UniversalLagrange(R0, V0, t);
    Ri(ind,:) = R;
    [Long(ind,:), Lat(ind,:)] = R2RA_Dec(R);
    ind = ind + 1;
    t = t + step;
end
% Earth 3D Plot
Earth3DPlot(3);
scatter3(Ri(:,1),Ri(:,2),Ri(:,3),'.r');
% Earth topographic map
EarthTopographicMap(4,820,420);
scatter(Long ,Lat,'.r');
text(280,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
 title('Satellite ground track');

Lagrange Universal Anomaly_03Lagrange Universal Anomaly_04

Equations from the book Orbital Mechanics for Engineering Students, Second Edition Aerospace_Engineering

 

Universal Kepler’s equation

In this example we solve the universal Kepler’s equation to determine universal anomaly after dt time with a given initial radius r0, velocity v0 and semimajor axis of a spacecraft.

clc;
clear all;
mu       = 398600;        % Earth’s gravitational parameter [km^3/s^2]
% Initial conditions
r0       = 12000;  %[km] Initial radius
vr0      = 3;      %[km/s] Initial radial speed
dt       = 7200;   %[s] Time interval
a        = -20000; %[km] Semimajor axis
alfa = 1/a;        %[km^-1] Reciprocal of the semimajor axis;

X0 = mu^0.5*abs(alfa)*dt;  %[km^0.5]Initial estimate of X0
Xi = X0;
tol = 1E-10;              % Tolerance
while(1)
    zi = alfa*Xi^2;
    [ Cz,Sz] = Stumpff( zi );
    fX  = r0*vr0/(mu)^0.5*Xi^2*Cz + (1 - alfa*r0)*Xi^3*Sz + r0*Xi -(mu)^0.5*dt;
    fdX = r0*vr0/(mu)^0.5*Xi*(1 - alfa*Xi^2*Sz) + (1 - alfa*r0)*Xi^2*Cz + r0;
    eps = fX/fdX;
    Xi = Xi - eps;
    if(abs(eps) < tol )
        break
    end
end
fprintf('Universal anomaly X = %4.3f [km^0.5] \n',Xi)

% Equations from the book
% Orbital Mechanics for Engineering Students,2nd Edition,Aerospace Engineering
Universal anomaly X = 173.324 [km^0.5]

Stumpff functions

The Stumpff functions C(z),S(z) developed by Karl Stumpff, are used for analyzing orbits using the universal variable formulation.

 clc;
 clear all;
 z = -25:0.1:400;
 n = size(z);
 for i = 1:n(2)
     if     (z(i) > 0)
         S(i) = (z(i)^0.5 -sin(z(i)^0.5))/z(i)^1.5;
         C(i) = (1 - cos(z(i)^0.5))/z(i);
     elseif (z(i) < 0)
         S(i) = (sinh((-z(i))^0.5) - (-z(i))^0.5)/((-z(i))^1.5);
         C(i) = (cosh((-z(i))^0.5) - 1)/(-z(i));
     else
         S(i) = 1/6;
         C(i) = 0.5;
     end
 end
% Plot
figure(1);
scatter(z(1:500),C(1:500),'.b');
hold on;grid on;
scatter(z(1:500),S(1:500),'.r');
xlabel('z');
ylabel('C(z),S(z)');
legend('C(z)','S(z)');
title('Stumpff functions C(z),S(z)');
text(10,0.4,'smallsats.org','Color',[0 0 0], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
hold off;

figure(2);
scatter(z(501:n(2)),C(501:n(2)),'.b');
hold on;grid on;
scatter(z(501:n(2)),S(501:n(2)),'.r');
xlabel('z');
ylabel('C(z),S(z)');
legend('C(z)','S(z)');
title('Stumpff functions C(z),S(z)');
text(250,0.0010,'smallsats.org','Color',[0 0 0], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
hold off;

Stumpff functionsExStumpff functionsEx

Equations from the book Orbital Mechanics for Engineering Students, Second Edition, Aerospace Engineering

%d bloggers like this: