Small Satellites

Home » Posts tagged 'Aerospace Engineering' (Page 3)

Tag Archives: Aerospace Engineering

Julian Day

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;

fprintf('Julian Day = %6.4f [days] \n',JD)
Julian Day = 2456394.3413 [days]

 

Sun-Synchronous Circular Orbit, Inclination vs Altitude (LEO,J2 perturbed)

This code is a MATLAB script that can be used to design and analyze Sun-synchronous orbits. A Sun-synchronous orbit is a geocentric orbit which combines altitude and inclination in such a way that an object in this orbit has an a nodal regression rate which is equals to Earth’s orbital rotation speed around the Sun. The object in this orbit constantly illuminated by the Sun.

Output: Inclination vs Altitude Plot

clc;
clear all;
mu    = 398600.440;      % Earth’s gravitational parameter [km^3/s^2]
Re = 6378;               % Earth radius [km]
J2  = 0.0010826269;      % Second zonal gravity harmonic of the Earth
we = 1.99106e-7;    % Mean motion of the Earth in its orbit around the Sun [rad/s]
% Input
Alt = 250:5:1000;     % Altitude,Low Earth orbit (LEO)
a   = Alt + Re;       % Mean semimajor axis [km]
e   = 0.0;            % Eccentricity
h = a*(1 - e^2);     % [km]
n = (mu./a.^3).^0.5; % Mean motion [s-1]
tol = 1e-10;         % Error tolerance
% Initial guess for the orbital inclination
i0 = 180/pi*acos(-2/3*(h/Re).^2*we./(n*J2));
err = 1e1;
while(err >= tol )
    % J2 perturbed mean motion
    np  = n.*(1 + 1.5*J2*(Re./h).^2.*(1 - e^2)^0.5.*(1 - 3/2*sind(i0).^2));
    i = 180/pi*acos(-2/3*(h/Re).^2*we./(np*J2));
    err = abs(i - i0);
    i0 = i;
end
plot(Alt,i,'.b');
grid on;hold on;
xlabel('Altitude,Low Earth orbit (LEO)');
ylabel('Mean orbital inclination');
title('Sun-Synchronous Circular Orbit,Inclination vs Altitude(LEO,J2 perturbed)');
hold off;

Sun_sync_vs_Alt

Sun-Synchronous (Heliosynchronous) Orbit, Mean Orbital Inclination (J2 perturbed)

This code is a MATLAB script that can be used to design and analyze Sun-synchronous orbits. A Sun-synchronous orbit is a geocentric orbit which combines altitude and inclination in such a way that an object in this orbit has an a nodal regression rate which is equals to Earth’s orbital rotation speed around the Sun. The object in this orbit constantly illuminated by the Sun. Input: a = Mean semimajor axis, e = Eccentricity

Output: i = Mean orbital inclination

clc;
clear all;
mu    = 398600.440;      % Earth’s gravitational parameter [km^3/s^2]
Re = 6378;               % Earth radius [km]
J2  = 0.0010826269;      % Second zonal gravity harmonic of the Earth
we = 1.99106e-7;    % Mean motion of the Earth in its orbit around the Sun [rad/s]
% Input
a = 7378;           % Mean semimajor axis [km]
e = 0.05;           % Eccentricity
h = a*(1 - e^2);    % [km]
n = (mu/a^3)^0.5;   % Mean motion [s-1]
tol = 1e-10;        % Error tolerance
% Initial guess for the orbital inclination
i0 = 180/pi*acos(-2/3*(h/Re)^2*we/(n*J2));
err = 1e1;
while(err >= tol )
    % J2 perturbed mean motion
    np  = n*(1 + 1.5*J2*(Re/h)^2*(1 - e^2)^0.5*(1 - 3/2*sind(i0)^2));
    i = 180/pi*acos(-2/3*(h/Re)^2*we/(np*J2));
    err = abs(i - i0);
    i0 = i;
end
fprintf('Mean orbital inclination %4.5f [deg] \n',i);
Mean orbital inclination 99.43667 [deg]

Gibbs method of preliminary orbit determination

This example will show how to use 3 earth-centered position vectors of a satellite R1,R2,R3 at successive times to determine satellite orbit.

clear all; clc;close all;
mu = 398600;               %[km^3/s^2] Earth’s gravitational parameter

R1 = [5887 -3520 -1204];   %[km]
R2 = [5572 -3457 -2376];   %[km]
R3 = [5088 -3289 -3480];   %[km]

r1 = norm(R1); r2 = norm(R2); r3 = norm(R3);
% Vector cross products
R12 = cross(R1,R2);
R23 = cross(R2,R3);
R31 = cross(R3,R1);

Nv  = r1*R23 + r2*R31 + r3*R12;
Dv  = R23 + R31 + R12;
Sv  = (r2-r3)*R1 + (r3-r1)*R2 + (r1-r2)*R3;

N   = norm(Nv);
D   = norm(Dv);
% Velocity vector
V1 = (mu/(N*D))^0.5*(cross(Dv,R1)/r1 + Sv);
[ h i RAAN e omega theta ] = State2Orbital(R1,V1);
OE = [h i RAAN e omega theta];

fprintf('h [km^2/s]    i [deg]     RAAN [deg]  e    omega[deg]   theta [deg] \n');
fprintf('%4.2f     %4.2f       %4.2f     %4.4f    %4.2f    %4.2f  \n',OE);
% Check, Only true anomaly
V2 = (mu/(N*D))^0.5*(cross(Dv,R2)/r2 + Sv);
[ h i RAAN e omega theta ] = State2Orbital(R2,V2);
OE = [h i RAAN e omega theta];

fprintf('h [km^2/s]    i [deg]     RAAN [deg]  e    omega[deg]   theta [deg] \n');
fprintf('%4.2f     %4.2f       %4.2f     %4.4f    %4.2f    %4.2f  \n',OE);
h [km^2/s]    i [deg]     RAAN [deg]  e    omega[deg]   theta [deg] 
52948.88     95.03       150.01     0.0127    151.69    38.30  
h [km^2/s]    i [deg]     RAAN [deg]  e    omega[deg]   theta [deg] 
52948.88     95.01       150.00     0.0127    151.69    48.31
EarthTopographicMap(1,820,420);
[Rv, alfa,delta ] = J2PropagR( h, i, RAAN, e,omega,theta) ;
scatter(alfa,delta,'.r');
text(270,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',...
	'HorizontalAlignment','left','FontSize',14 );
title('Satellite ground track');
% Earth 3D Plot
Earth3DPlot(2);
scatter3(Rv(:,1),Rv(:,2),Rv(:,3),'.r');
% Plot position vectors
% Direction Cosines
l = R1/r1; m = R2/r2; n = R3/r3;
lv =(0:1:r1)'*l; mv =(0:1:r2)'*m;nv =(0:1:r3)'*n;
figure(3);
scatter3(Rv(:,1),Rv(:,2),Rv(:,3),'.k');
hold on;
scatter3(lv(:,1),lv(:,2),lv(:,3),'.r');
scatter3(mv(:,1),mv(:,2),mv(:,3),'.g');
scatter3(nv(:,1),nv(:,2),nv(:,3),'.b');
legend('Orbit','R1','R2','R3');
title('Satellite Orbit and R1,R2,R3 successive position vectors');
zoom(2)
% Equations from the book
% Orbital Mechanics for Engineering Students, 2nd Edition, Aerospace EngineeringGibbsMethod_01

 

Universal Kepler’s equation

In this example we solve the universal Kepler’s equation to determine universal anomaly after dt time with a given initial radius r0, velocity v0 and semimajor axis of a spacecraft.

clc;
clear all;
mu       = 398600;        % Earth’s gravitational parameter [km^3/s^2]
% Initial conditions
r0       = 12000;  %[km] Initial radius
vr0      = 3;      %[km/s] Initial radial speed
dt       = 7200;   %[s] Time interval
a        = -20000; %[km] Semimajor axis
alfa = 1/a;        %[km^-1] Reciprocal of the semimajor axis;

X0 = mu^0.5*abs(alfa)*dt;  %[km^0.5]Initial estimate of X0
Xi = X0;
tol = 1E-10;              % Tolerance
while(1)
    zi = alfa*Xi^2;
    [ Cz,Sz] = Stumpff( zi );
    fX  = r0*vr0/(mu)^0.5*Xi^2*Cz + (1 - alfa*r0)*Xi^3*Sz + r0*Xi -(mu)^0.5*dt;
    fdX = r0*vr0/(mu)^0.5*Xi*(1 - alfa*Xi^2*Sz) + (1 - alfa*r0)*Xi^2*Cz + r0;
    eps = fX/fdX;
    Xi = Xi - eps;
    if(abs(eps) < tol )
        break
    end
end
fprintf('Universal anomaly X = %4.3f [km^0.5] \n',Xi)

% Equations from the book
% Orbital Mechanics for Engineering Students,2nd Edition,Aerospace Engineering
Universal anomaly X = 173.324 [km^0.5]
%d bloggers like this: