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]
Small Satellites 