Home » 2013 » February

# Monthly Archives: February 2013

## Rodrigues Parameters from Rotation Matrix

In this example we derive Rodrigues parameters/Gibbs Vector from the rotation matrix elements. The Rodrigues parameters have a singularity at 180 deg and use is limited for principal rotations which are less than 180 deg.

```clc; clear all;
% Rotation matrix
Rm = [0.36   0.48  -0.8;
-0.8    0.6    0;
0.48   0.64   0.6 ];

g1 = (Rm(2,3)- Rm(3,2))/(1 + trace(Rm));
g2 = (Rm(3,1)- Rm(1,3))/(1 + trace(Rm));
g3 = (Rm(1,2)- Rm(2,1))/(1 + trace(Rm));
G = [g1 g2 g3]'```
```G =

-0.2500
0.5000
0.5000```

## Non-impulsive Orbital Maneuvers, Constant Thrust

In non-impulsive orbital maneuvers the thrust acts over a significant period of time. In this case we can’t ignore the change in position vector. Given the initial state of the spacecraft we calculate the state after the burn

```clc;
clear all;
global g0 mu T Isp M0;
mu = 398600;             %[km^3/s^2] Earth’s gravitational parameter
M0 = 4000;               %[kg] Gross mass before ignition
R = [ 7200,   0,  0];    %[km] Position vector:
V = [ 0,     7.6,  0];   %[km/s] Velocity vector
g0 = 9.81;               %[km/s^2]
T   = 20;                %[kN] Thrust, const
t_b =  200;              %[s] Burn time
m =  1000;               %[kg] Mass after burn
dm = (M0 - m)/t_b;
Isp = T/(dm*g0);
tspan = [0,t_b];
% Constructing initial state vector
y0 = [R(1),R(2),R(3),V(1),V(2),V(3),M0];
% Using ODE45 to solve  ordinary differential equations of motion
% for the given initial values
[time,f_y] = ode45(@f_yt,tspan,y0);
fs = f_y(end,:);
Rf = [fs(1),fs(2),fs(3)];
Vf = [fs(4),fs(5),fs(6)];
Mf = fs(7);
fprintf('Position vector after burn R = [%4.2f %4.2f %4.2f] km\n',Rf);
fprintf('Velocity vector after burn V = [%4.4f %4.2f %4.2f] km/s\n',Vf);
fprintf('Mass after burn M = %4.2f kg\n',Mf);```
```Position vector after burn R = [7035.78 1651.64 0.00] km
Velocity vector after burn V = [-1.7367 9.26 0.00] km/s
Mass after burn M = 1000.00 kg```
```function  dfdt = f_yt(t,y)
global mu g0 T Isp;
r = (y(1)^2+y(2)^2+y(3)^2)^0.5;
v = (y(4)^2+y(5)^2+y(6)^2)^0.5;
m = y(7);
dfdt = [y(4),y(5),y(6),-mu*y(1)/r^3 + T/m*y(4)/v,...
-mu*y(2)/r^3 + T/m*y(5)/v, -mu*y(3)/r^3 + T/m*y(6)/v,...
-T/(Isp*g0)]';
return```

## Hohmann vs Bi-elliptic transfer

Contents

```% In this example we compare be-elliptic and Hohmann transfer.
% The total speed change that required for spacecraft transfer from
% geocentric circular orbit with radius Ri to a higher altitude Rf
clc;
clear all;
Rf  =  125000;   % [km] Final circular orbit
Ri  =  7200;     % [km] Initial circular orbit
Rb  =  190000;   % [km] Apogee of the transfer ellipse
mu  =  398600;   % [km^3/s^2] Earth’s gravitational parameter

% For initial circular orbit
Vc = (mu/Ri)^0.5;
a   = Rf/Ri;
b   = Rb/Ri;```

## Hohmann transfer

For Hohmann transfer total speed change

```dV_H =Vc*(1/(a)^0.5 -(2/(a*(a+1)))^0.5*(1-a) - 1);
% Semimajor axes of the Hohmann transfer ellipse
a_h = (Rf + Ri)/2;
% Time required for Hohmann transfer
t_H = pi/(mu)^0.5*(a_h^1.5);     %[s]
fprintf('Total speed change = %4.4f [km/s]\n',dV_H);
fprintf('Time required for transfer = %4.2f [hours]\n\n',t_H/3600);```
```Total speed change = 3.9878 [km/s]
Time required for transfer = 23.49 [hours]```

## Bi-elliptic transfer

For Bi-elliptic transfer total speed change

```dV_BE = Vc*((2*(a+b)/(a*b))^0.5 - (1+1/a^0.5) - ((2/(b +b^2))^0.5*(1-b)));
% Semimajor axes of the first transfer ellipse
a1 = (Ri + Rb)/2;
% Semimajor axes of the second transfer ellipse
a2 = (Rf + Rb)/2;
t_BE = pi/(mu)^0.5*(a1^1.5+a2^1.5);     %[s]
fprintf('Total speed change = %4.4f [km/s]\n',dV_BE);
fprintf('Time required for transfer = %4.2f [hours]\n',t_BE/3600);```
```Total speed change = 3.9626 [km/s]
Time required for transfer = 129.19 [hours]```

## Bi-Elliptic Hohmann Transfer

In this example we calculate the total change in speed required for a bi-elliptic Hohmann transfer from a geocentric circular orbit of 7200 km radius to circular orbit of 125000 km radius. The apogee of the first transfer ellipse is 190000 km.

``` clc;
clear all;

R_i  = 7200;     % [km]
R1_a = 190000;   % [km]
R_f  = 125000;   % [km]

mu   = 398600;   % [km^3/s^2] Earth’s gravitational parameter
% For initial circular orbit
V_i = (mu/R_i)^0.5;
% Speed at apogee and perigee for the first transfer ellipse
V1_a = (2*mu*R_i/(R1_a*(R1_a+R_i)))^0.5;
V1_p = (2*mu*R1_a/(R_i*(R1_a+R_i)))^0.5;
% Semimajor axes of the first transfer ellipse
a1 = (R_i + R1_a)/2;
% Speed at apogee and perigee for the second transfer ellipse
V2_a = (2*mu*R_f/(R1_a*(R1_a+R_f)))^0.5;
V2_p = (2*mu*R1_a/(R_f*(R1_a + R_f)))^0.5;
% Semimajor axes of the second transfer ellipse
a2 = (R_f + R1_a)/2;
% For target circular orbit
V_f = (mu/R_f)^0.5;

% For bi-elliptic maneuver the total speed change required
dV =  abs(V_i - V1_p)+ abs(V1_a - V2_a) + abs(V_f - V2_p); %[km/s]

% Time required for transfer
t_bi = pi/(mu)^0.5*(a1^1.5+a2^1.5);     %[s]

fprintf('Total speed change = %4.4f [km/s]\n',dV);
fprintf('Time required for transfer = %4.2f [hours]\n',t_bi/3600);```
```Total speed change = 3.9626 [km/s]
Time required for transfer = 129.19 [hours]```