Home » Space Flight/Orbital Mechanics (Page 2)

# Category Archives: Space Flight/Orbital Mechanics

## Orbital Inclination Change

Transfer from a LEO 350 km circular orbit with 53.4 deg inclination to a Geostationary Equatorial Orbit(GEO)

```clear all; clc;
close all;```

## Input

```R_LEO = 6378 + 350;      % km
R_GEO = 42164;           % km
mu    = 398600;          % km^3/s^2
incl  = 53.4;            % deg```

## Method 1: Hohman transfer from LEO to GEO and after inclination change

```Rp  = R_LEO;
Ra  = R_GEO;
e   = (Ra - Rp)/(Ra + Rp);  % transfer orbit eccentricity
a   = (Ra + Rp)/2;          % transfer orbit semimajor axis
V_LEO   = (mu/R_LEO)^0.5;
Vp      = (2*mu/R_LEO - mu/a)^0.5;
Va      = (2*mu/R_GEO - mu/a)^0.5;
V_GEO   = (mu/R_GEO)^0.5;
dV_LEO  = abs(Vp - V_LEO);
dV_GEO  = abs(V_GEO - Va);
dV_Hoff = dV_GEO + dV_LEO;

% Inclination change at GEO
dV_incl = 2*V_GEO*sind(incl/2);
dV_total1 = dV_incl + dV_Hoff;

fprintf('Method 1: Hohman transfer from LEO to GEO and after inclination change\n');
fprintf('dV_Hoff  = %6.2f [km/s] \n',dV_Hoff);
fprintf('dV_incl  = %6.2f [km/s] \n',dV_incl);
fprintf('dV_total = %6.2f [km/s] \n',dV_total1);```
```Method 1: Hohman transfer from LEO to GEO and after inclination change
dV_Hoff  =   3.87 [km/s]
dV_incl  =   2.76 [km/s]
dV_total =   6.64 [km/s]```

## Method 2: Inclination change in LEO and after Hohman transfer to GEO

```dV_incl   = 2*V_LEO*sind(incl/2);
dV_total2 = dV_incl + dV_Hoff;

fprintf('\nMethod 2: Inclination change in LEO and after Hohman transfer to GEO\n');
fprintf('dV_incl  = %6.2f [km/s] \n',dV_incl);
fprintf('dV_Hoff  = %6.2f [km/s] \n',dV_Hoff);
fprintf('dV_total = %6.2f [km/s] \n\n',dV_total2);
fprintf('Method 2/Method 1 = %6.2f \n',dV_total2/dV_total1);```
```Method 2: Inclination change in LEO and after Hohman transfer to GEO
dV_incl  =   6.92 [km/s]
dV_Hoff  =   3.87 [km/s]
dV_total =  10.79 [km/s]

Method 2/Method 1 =   1.63```

## Kepler Equation Solver Without Transcendental Function Evaluations

clc;clear all;

```tic;
j = 1;
for M = 0:0.005:pi
i = 1;
for e = 0:0.005:1
c3 = 5/2 + 560*e;
a = 15/c3*(1-e);
b = -M/c3;
y = (b^2/4+ a^3/27)^0.5;
x = (-b/2 + y)^(1/3) - (b/2 + y)^(1/3);
c15 = 3003/14336 + 16384*e;
c13 = 3465/13312 - 61440*e;
c11 = 945/2816 + 92160*e;
c9 = 175/384 -70400*e;
c7 = 75/112 + 28800*e;
c5 = 9/8-6048*e;
x2 = x^2;x3 = x2*x;x4 = x3*x;
x5 = x4*x;x6 = x5*x;x7 = x6*x;
x8 = x7*x;x9 = x8*x;
x10 =x9*x; x11 = x10*x; x12 = x11*x;
x13 = x12*x; x14 = x13*x; x15 = x14*x;

f = c15*x15 + c13*x13 + c11*x11 + c9*x9 + c7*x7 + c5*x5 + c3*x3 +...
15*(1 - e)*x -M;
f1 = 15* c15 *x14 + 13* c13 *x12 + 11* c11 *x10 + 9* c9 *x8 + 7* c7 *x6 + ...
5 *c5 *x4 + 3 *c3 *x2 + 15*(1 - e);
f2 = 210* c15 *x13 + 156* c13 *x11 + 110* c11 *x9 + 72* c9 *x7 + ...
42* c7 *x5 + 20* c5 *x3 + 6* c3 *x;
f3 = 2730*c15 *x12 + 1716* c13 *x10 + 990* c11 *x8 + 504* c9 *x6 +...
210* c7 *x4 + 60* c5 *x2 + 6* c3;
f4 = 32760* c15 *x11 + 17160* c13 *x9 + 7920* c11 *x7 + 3024* c9 *x5 +...
840* c7 *x3 + 120* c5 *x;
f5 = 360360* c15 *x10 + 154440* c13 *x8 + 55440* c11 *x6 + 15120* c9 *x4 ...
+ 2520* c7 *x2 + 120* c5;
f6 = 3603600* c15 *x9 + 1235520* c13 *x7 + 332640* c11 *x5 + 60480* c9 *x3...
+ 5040* c7 *x;
f7 = 32432400* c15 *x8 + 8648640* c13 *x6 + 1663200* c11 *x4 + 181440* c9 *x2 ...
+ 5040* c7;
f8 = 25945920* c15 *x7 + 51891840* c13 *x5 + 6652800* c11 *x3 + ...
362880* c9 *x;
f9 = 1.8162144e9* c15 *x6 + 259459200* c13 *x4 + 19958400* c11 *x2 + ...
362880* c9;
f10 = 1.08972864e10* c15 *x5 + 1.0378368e9*c13 *x3 + 39916800* c11 *x;
f11 = 5.4486432e10* c15 *x4 + 3.1135104e9* c13 *x2 + 39916800* c11;
f12 = 2.17945728e11* c15 *x3 + 6.2270208e9* c13 *x;
f13 = 6.53837184e11* c15 *x2 + 6.2270208e9* c13;
f14 = 1.307674368e13* c15 *x;
f15 = 1.307674368e13* c15;

g1 = 1/2; g2 = 1/6; g3 = 1/24; g4 = 1/120;
g5 = 1/720; g6 = 1/5040; g7 = 1/40320; g8 = 1/362880;
g9 = 1/3628800; g10 = 1/39916800; g11 = 1/479001600;
g12 =1/6.2270208e9 ; g13 = 1/8.71782912e10; g14 = 1/1.307674368e12;

u1 = -f/f1;
h2 = f1 + g1* u1* f2;
u2 = -f/h2;
h3 = f1 + g1* u2* f2 + g2* u2^2*f3;
u3 = -f/h3;
h4 = f1 + g1* u3* f2 + g2* u3^2*f3 + g3* u3^3*f4;
u4 = -f/h4;
h5 = f1 + g1* u4* f2 + g2* u4^2*f3 + g3* u4^3*f4 + g4* u4^4*f5;
u5 = -f/h5;
h6 = f1 + g1* u5* f2 + g2* u5^2*f3 + g3* u5^3*f4 + g4* u5^4*f5 +...
g5* u5^5*f6;
u6 = -f/h6;
h7 = f1 + g1* u6* f2 + g2* u6^2*f3 + g3* u6^3*f4 + g4* u6^4*f5 +...
g5* u6^5*f6 + g6* u6^6*f7;
u7 = -f/h7;
h8 = f1 + g1* u7* f2 + g2* u7^2*f3 + g3* u7^3*f4 + g4* u7^4*f5 +...
g5* u7^5*f6 + g6* u7^6*f7 + g7* u7^7*f8;
u8 = -f/h8;
h9 = f1 + g1* u8* f2 + g2* u8^2*f3 + g3* u8^3*f4 + g4* u8^4*f5 +...
g5* u8^5*f6 + g6* u8^6*f7+ g7* u8^7*f8 + g8* u8^8*f9;
u9 = - f/h9;
h10 = f1 + g1* u9* f2 + g2* u9^2*f3 + g3* u9^3*f4 + g4* u9^4*f5 +...
g5* u9^5*f6 + g6* u9^6*f7+ g7* u9^7*f8 + g8* u9^8*f9 + ...
g9* u9^9*f10;
u10 = -f/h10;
h11 = f1 + g1* u10* f2 + g2* u10^2* f3 + g3* u10^3* f4 + g4* u10^4*f5...
+ g5* u10^5* f6 + g6* u10^6*f7+ g7* u10^7* f8 + g8* u10^8* f9...
+ g9* u10^9* f10 + g10* u10^10* f11;
u11 = -f/h11;
h12 = f1 + g1* u11* f2 + g2* u11^2* f3 + g3* u11^3* f4 + ...
g4* u11^4* f5 + g5* u11^5* f6 + g6* u11^6* f7...
+ g7* u11^7* f8 + g8* u11^8* f9 + g9* u11^9* f10 + g10* u11^10* f11...
+ g11* u11^11* f12;
u12 = -f/h12;
h13 = f1 + g1* u12* f2 + g2* u12^2* f3 + g3* u12^3* f4 +...
g4* u12^4* f5 + g5* u12^5* f6 + g6* u12^6* f7...
+ g7* u12^7* f8 + g8* u12^8* f9 + g9* u12^9* f10 +...
g10* u12^10* f11 + g11* u12^11* f12 + g12* u12^12* f13;
u13 = -f/h13;
h14 = f1 + g1* u13* f2 + g2* u13^2* f3 + g3* u13^3* f4 + ...
g4* u13^4* f5 + g5* u13^5* f6 + g6* u13^6* f7...
+ g7* u13^7* f8 + g8* u13^8* f9 + g9* u13^9* f10 + g10* u13^10*f11...
+ g11* u13^11* f12+ g12* u13^12* f13 + g13* u13^13* f14;
u14 = -f/h14;
h15 = f1 + g1* u14* f2 + g2* u14^2* f3 + g3* u14^3* f4 + ...
g4* u14^4* f5 + g5* u14^5* f6 + g6* u14^6* f7...
+ g7* u14^7* f8 + g8* u14^8* f9 + g9* u14^9* f10 +...
g10* u14^10* f11 + g11*u14^11* f12+ g12* u14^12* f13 + ...
g13* u14^13* f14 + g14* u14^14* f15;
u15 = -f/h15;

x =  x + u15;
w = x - 0.01171875*x^17/(1 + e);
E = M + e*(-16384*w^15 + 61440*w^13 - 92160*w^11+ 70400*w^9 -...
28800*w^7 + 6048*w^5 - 560*w^3 + 15*w);
E_nit(i,j) = E;
i = i + 1;
end
j = j + 1;
end
t_nit = toc;
fprintf('Calculation time = %4.6f [s] \n',t_nit);```
```Calculation time = 2.591335 [s]
```

## Newton Iterative Solver

```tic;
j = 1;
eps = 1E-15;        % Tolerance
for M = 0:0.005:pi
i = 1;
for e = 0:0.005:1
En  = M;
Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En));
while ( abs(Ens-En) > eps )
En = Ens;
Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En));
end;
E_it(i,j) = Ens;
i = i + 1;
end
j = j + 1;
end
t_it = toc;
fprintf('Calculation time = %4.6f [s] \n',t_it);```
`Calculation time = 3.520123 [s]`
```Error = log10(abs(E_it - E_nit));
e = 0:0.005:1;
M = (0:0.005:pi);
contourf (M,e,Error);
ylabel('e - Eccentricity');
title('Error in E calculation on Log scale [rad] ');
colorbar;

REFERENCES
Adonis Pimienta-Pe˜nalver and John L. Crassidis†,ACCURATE KEPLER EQUATION
SOLVER WITHOUT TRANSCENDENTAL FUNCTION EVALUATIONS,
University at Buffalo, State University of New York, Amherst, NY, 14260-4400```

## Kepler’s equation, Iterative and Non-Iterative Solver Comparison

Kepler’s equation M = E-e*sin(E) M – Mean anomaly [0..pi] rad e – Eccentricity [0..1] E – Eccentric anomaly [rad]

```clear all;
clc;```

## Non-Iterative Solver

```tic;
j = 1;
for M = 0:0.01:pi
i = 1;
for e = 0:0.01:1
a   =  (1 - e)*3/(4*e+0.5);
b   =  -M/(4*e+0.5);
y   =  (b^2/4 + a^3/27)^0.5;
x   =  (-b/2 + y)^(1/3) - (b/2 + y)^(1/3);
w   =  x - 0.078*x^5/(1 + e);
E   =  M + e*(3*w - 4*w^3);
% Two times applying 5th-order Newton correction
for k = 1:2
f   =  (E - e*sin(E)- M);
fd  =  1 - e*cos(E);
f2d =  e*sin(E);
f3d = -e*cos(E);
f4d =  e*sin(E);
E = E -f/fd*(1 + f*f2d/(2*fd^2) + f^2*(3*f2d^2 - fd*f3d)/(6*fd^4)+...
(10*fd*f2d*f3d - 15*f2d^3-fd^2*f4d)*f^3/(24*fd^6));
end
E_nit(i,j) = E;
i = i + 1;
end
j = j + 1;
end
t_nit = toc;
fprintf('Calculation time = %4.6f [s] \n',t_nit);
Calculation time = 0.149860 [s]```

## Newton Iterative Solver

```tic;
j = 1;
eps = 1E-15;        % Tolerance
for M = 0:0.01:pi
i = 1;
for e = 0:0.01:1
En  = M;
Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En));
while ( abs(Ens-En) > eps )
En = Ens;
Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En));
end;
E_it(i,j) = Ens;
i = i + 1;
end
j = j + 1;
end
t_it = toc;
fprintf('Calculation time = %4.6f [s] \n',t_it);```
`Calculation time = 0.893673 [s]`
```Error = log10(abs(E_it - E_nit));
e = 0:0.01:1;
M = (0:0.01:pi);
contourf (M,e,Error);
ylabel('e - Eccentricity');
title('Non-Iterative Solver,Error in E calculation on Log scale [rad] ');
colorbar;

[6] S. Mikkola, “A Cubic Approximation for Kepler’s Equation,” Celestial
Mechanics and Dynamical Astronomy,Vol. 40, 1987, pp. 329–334.```

## Local Sidereal Time

This code is a MATLAB script that can be used to calculate Greenwich Sidereal Time, Local Sidereal Time, Julian Day

```clc;
clear all;```

## Input:

Date, April 11,2013, UT time 20:11:30, Longitude [deg]

## Output:

Julian day(JD), Greenwich sidereal time(GST), Local sidereal time(LST)

```year  = 2013;  month = 4;  day   = 11;
hour  = 20;    min   = 11; sec   = 30;
long  = -73.99;```

Julian day

```UT = hour + min/60 + sec/3600;
J0 = 367*year - floor(7/4*(year + floor((month+9)/12))) ...
+ floor(275*month/9) + day + 1721013.5;
JD = J0 + UT/24;              % Julian Day
fprintf('Julian day = %6.4f [days] \n',JD);```
`Julian day = 2456394.3413 [days]`

JC is Julian centuries between the Julian day J0 and J2000(2,451,545.0) Greenwich sidereal time at 0 hr UT can be calculated by this equation [Seidelmann,1992]

```JC = (J0 - 2451545.0)/36525;
GST0 = 100.4606184 + 36000.77004*JC + 0.000387933*JC^2 - 2.583e-8*JC^3; %[deg]
GST0 = mod(GST0, 360);  % GST0 range [0..360]
fprintf('Greenwich sidereal time at 0 hr UT %6.4f [deg]\n',GST0);```
`Greenwich sidereal time at 0 hr UT 199.3719 [deg]`

Greenwich sidereal time at any other UT time

```GST = GST0 + 360.98564724*UT/24;
GST = mod(GST, 360);  % GST0 range [0..360]
fprintf('Greenwich sidereal time at UT[hours] %6.4f [deg]\n',GST);```
`Greenwich sidereal time at UT[hours] 143.0761 [deg]`

Local sidereal time (LST)

```LST = GST + long;
LST = mod(LST, 360);  % LST range [0..360]
fprintf('Local sidereal time,LST %6.4f [deg]\n',LST);```
`Local sidereal time,LST 69.0861 [deg]`

## Azimuth and Elevation

This code is a MATLAB script that can be used to calculate spacecraft Azimuth and Elevation angles relative to observer position.

## Input:

R_sc – Spacecraft geocentric equatorial position vector [km],

H – Observer elevation from Sea level [km],

lat – Observer latitude [deg],

lst – Local sidereal time [deg]

## Output:

Elev – Elevation angle [deg],

Az – Azimuth angle [deg]

```clc; clear all;
R_sc    = [-2000;4500;-4500];
H       = 0.42;
lat     = 40.5;
lst     = 90.5;

Re      = 6378.137;     % Equatorial Earh's radius [km]
Rp      = 6356.7523;    % Polar Earh's radius [km]
f       = (Re - Rp)/Re; % Oblateness or flattening```
```C1   = (Re/(1 - (2*f - f^2)*sind(lat)^2)^0.5 + H)*cosd(lat);
C2   = (Re*(1 - f)^2/(1 - (2*f - f^2)*sind(lat)^2)^0.5 + H)*sind(lat);
% Position vector of the observer,GEF
R_ob = [C1*cosd(lst); C1*sind(lst);C2];
% Position vector of the spacecraft relative to the observer
R_rel = R_sc - R_ob;

% GE_TH is direction cosine matrix to transform position vector components
% from geocentric equatorial frame into the topocentric horizon fream

GE_TH = [-sind(lst)          cosd(lst)              0;
-sind(lat)*cosd(lst) -sind(lat)*sind(lst)  cosd(lat);
cosd(lat)*cosd(lst)  cosd(lat)*sind(lst)   sind(lat)
];
R_rel_TH = GE_TH*R_rel;
rv = R_rel_TH/norm(R_rel_TH);
Elev = asin(rv(3))*180/pi;      % Elevation angle
Az  =atan2(rv(1),rv(2))*180/pi; % Azimuth angle```
```fprintf('Elevation angle =  %4.2f [deg] \n',Elev);
fprintf('Azimuth angle   =  %4.2f [deg] \n',Az);```
```Elevation angle =  -41.45 [deg]
Azimuth angle   =  162.80 [deg]```