Open source software for numerical computation

Skip to main content

Scilab @ ICATT 2016 --> Orbit Simulation Tutorial

6th_ICATT_logo_guillermo

Following the advice of our long-time partners from the French Space Agency CNES, we are presenting today at the 6th International Conference on Astrodynamics Tools and Techniques (ICATT).

To demonstrate the capacities of Scilab in this field, here is a classical problem of space mechanics solved step-by-step:

Satellite orbit around the earth

Download the script: 

Download the full tutorial:

1.   Express the physics problem 

The problem is based on the universal law of gravitation:

eq1

We write down Newton’s third law of motion in an earth-centred referential:

(1)

eq2

(2)

eq3
referential

Position of the satellite is at a distance r [x; y]

Earth mass centre is at O [0; 0]

Gravitational constant: 

G

Mass of the earth:

M

Radius of the Earth:

r_earth

2.   Translate your problem into Scilab

Scilab is a matrix-based language. Instead of expressing the system as set of 4 independent equations (along the x and y axis, for position and speed), we describe it as a single matrix equation, of dimension 4x4:

This method is a classical trick to switch from a second order scalar differential equation to a first order matrix differential equation.

eq4

with

A
u

To simplify the equation, we define the variable:

c

Open scinotes with edit myEarthRotation.sci

Define the skeleton of the function:

function udot=f(t, u)
G = 6.67D-11; //Gravitational constant
M = 5.98D24; //Mass of the Earth
c = -G * M;
r_earth = 6.378E6; //radius of the Earth
r = sqrt(u(1)^2 + u(2)^2);
// Write the relationhsip between udot and u
if r < r_earth then
udot = [0 0 0 0]';
else
A = [[0 0 1 0];
[0 0 0 1];
[c/r^3 0 0 0];
[0 c/r^3 0 0]];
udot = A*u;
end
endfunction

The condition defined by the distance r of the satellite with the centre of earth stops the simulation if it’s colliding with earth’s surface.

crash

Try out the final script with the following initial conditions in speed and altitude:

 --> geo_alt = 35784; // in kms

 --> geo_speed = 1074; // in m/s

 --> simulation_time = 24; // in hours

 --> U = earthrotation(geo_alt, geo_speed, simulation_time);

3.   Compute the results and create a visual animation

With this function, we go to the core of the problem: 

function U=earthrotation(altitude, v_init, hours)
// altitude given in km
// v_init is a vector [vx; vy] given in m/s
// hours is the number of hours for the simulation
r_earth = 6.378E6;
altitude = altitude * 1000;
U0 = [r_earth + altitude; 0; 0; v_init];
t = 0:10:(3600*hours); // simulation time, one point every 10 seconds
U = ode(U0, 0, t, f);

// Draw the earth in blue
angle = 0:0.01:2*%pi;
x_earth = 6378 * cos(angle);
y_earth = 6378 * sin(angle);
fig = scf();
a = gca();
a.isoview = "on";
plot(x_earth, y_earth, 'b--');
plot(0, 0, 'b+');

// Draw the trajectory computed
comet(U(1,:)/1000, U(2,:)/1000, "colors", 3);
endfunction

The resolution of the ordinary differential equation (ODE) is computed with the Scilab function ode.

ode solves Ordinary Different Equations defined by:

y

where y is a real vector or matrix

The simplest call of ode is: y = ode(y0,t0,t,f) where y0 is the vector of initial conditions, t0 is the initial time, t is the vector of times at which the solution y is computed and y is matrix of solution vectors y=[y(t(1)),y(t(2)),...].

Go further

To go further in numerical analysis, find out more about the solvers:

Ordinary Differential Equations with Scilab, WATS Lectures, Université de Saint-Louis, G. Sallet, 2004

You can also model and simulate it with Xcos:

img