Home » Posts tagged 'Orbital Mechanics' (Page 2)

# Tag Archives: Orbital Mechanics

## 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);
xlabel('M - Mean anomaly [rad]');
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]```

## Greenwich Sidereal Time

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

```clc;
clear all;
% Input Date: April 11,2013.  UT time 20:11:30
year  = 2013;  month = 4; day   = 11;
hour  = 20; min   = 11; sec   = 30;```
```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);
JC = (J0 - 2451545.0)/36525;```
`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]

```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]`

## Reduced,Modified,Truncated,Dublin Julian Days, Mars Solar Date

Julian Day is defined as the number of days since noon UT on January 1, 4713 BC.

```clc;
clear all;
% Input Date: April 11,2013.  UT time 20:11:30
year  = 2013;  month = 4; day   = 11;
hour  = 20; min   = 11; sec   = 30;```
```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```
```RJD = JD - 2400000;           % Reduced JD
MJD = JD - 2400000.5;         % Modified JD, Introduced by SAO in 1957
TJD = JD - 2440000.5;         % Truncated JD, Introduced by NASA in 1979
DJD = JD - 2415020;           % Dublin JD, Introduced by the IAU in 1955
MSD = (JD - 2405522)/1.02749; % Mars Solar Date```
```fprintf('Julian Day   = %6.4f [days] \n',JD)
fprintf('Reduced JD   = %6.4f [days] \n',RJD)
fprintf('Modified JD  = %6.4f [days] \n',MJD)
fprintf('Truncated JD = %6.4f [days] \n',TJD)
fprintf('Dublin JD    = %6.4f [days] \n',DJD)
fprintf('Mars Solar Date = %6.4f \n',MSD)```
```Julian Day   = 2456394.3413 [days]
Reduced JD   = 56394.3413 [days]
Modified JD  = 56393.8413 [days]
Truncated JD = 16393.8413 [days]
Dublin JD    = 41374.3413 [days]
Mars Solar Date = 49511.2763```