## Felicia_Brimigion_M​odule_6_HW4

Version 1.0.0 (2.8 KB) by
Homework 4 Module 6

Updated 7 May 2022

%Felicia Brimigion
%SSA
%Spring 2022
%HOMEWORK 4-MODULE 6
close all
clear all
clc
%% PROBLEM 1
%We are using drag for this assignment since it is the most straight
%forward perturbation force and easiest to work with.
tic
%Initial conditions
alt=250; %altitude in km
r=alt+r_Earth; %altitude of object from Earth's center in km
r_bar=[0;r;0]; %position vector of object in km
G=6.674e-11; %Universal gravitational constant (N*m^2)/kg^2
m_Earth=5.972e24; %mass of Earth kg
v=sqrt((G*m_Earth)/(r*1000)); %velocity of object relatvie to Earth in m/s
v_bar=[-v;0;0]/1000; %v_bar vector, velocity vector of object km/s @ t=0
omega_E=[0;0;360]; %Earth's rotation rate, in degrees, per sidereal day
SDR=86164.1; %sidereal day in seconds
omega_E=omega_E/SDR; %degrees
mu=3.986e5; %km^3/s^2 %the graviational parameter of Earth ~=GM where G=6.674E-11 in Nm^2/kg^2 and M>>m
x_0=[r_bar;v_bar];
%Defining the object's characteristics
m_o=2.5; %Object's mass in kg
A=0.5^2; %Object's area in m^2
C_D=1; %drag coefficient
%Calculating the object's density and velocity from table provided in lecture
rho_0=2.789e-10; % reference density in kg/m^3
h_0=200; % reference altitude in km
H=37.105; % scale altitude in km
rho=rho_0*exp(-(alt-h_0)/H); %density in kg/m^3
v_bar_rel=v_bar-cross(omega_E,r_bar);
v_rel=norm(v_bar_rel); %returning vector magnitude
%Calulcating a_bar_drag for the object orbiting Earth
a_bar_drag=-(1/2)*((C_D*A)/m_o)*rho*v_rel*v_bar_rel;
%Setting Time
t=[0,1.577e8]; %5 years in seconds
ACL=@(r)-(mu*r)/norm(r)^3; %propagation of object's motion due to 2-body&drag forces (r vector containing x,y,x)
x_dot=@(t,x)[x(4:6);ACL(x(1:3))+(1/2)*((C_D*A)/m_o)*rho*(x(4:6)-cross(omega_E,r_bar))*norm(x(4:6)-cross(omega_E,r_bar))];
[t,x_t]=ode45(x_dot,t,x_0,odeset('RelTol',1.00e-8)); %Setting relative tolerance to 1e-8
%Plotting required figures as requested in prompt.
figure('Name','X/Y vs Time')
plot(t/3600,x_t(:,1:2))
title('X/Y Position vs Time')
xlabel('Time (Hours)')
ylabel('Position (km)')
figure('Name','Top View of the Object''s Orbit')
plot(x_t(:,1),x_t(:,2))
title('Top View of Objects Orbit')
xlabel('Position (km)')
ylabel('Position (km)')
figure('Name','Full Object''s Orbit around Earth')
plot3(x_t(:,1),x_t(:,2),x_t(:,3),'-r') %3D plot
title('Full Object''s Orbit around Earth')
axis equal
hold on
[X,Y,Z]=sphere(100);
surf(X*r_Earth,Y*r_Earth,Z*r_Earth,'EdgeColor','none'); %Surface plot, turning off edge color as to not inundate the plot with black
xlabel('Position (km)')
ylabel('Position (km)')
zlabel('Position (km)')
fprintf(['\tBy the plots shown, one can see that the object is actually',...
' reducing its orbit. This is illustrated by the thickness of the line becoming smaller.'])
%% PROBLEM 2
%Initial Conditions
alt=200; %km
r_Earth=6378; %km
r=alt+r_Earth; %km
r_bar=[0;r;0]; %km
G=6.674e-11; %(N*m^2)/kg^2
m_Earth=5.972e24; %kg
v=sqrt((G*m_Earth)/(r*1000)); %m/s
v_bar=[0;0;v]/1000; %km/s @ t=0
omega_E=[0;0;360]; %degrees per sidereal day
SDR=86164.1; %sidereal day in seconds
omega_E=omega_E/SDR; %degrees
mu=3.986e5; %km^3/sec^2
x_0=[r_bar;v_bar];%x_0 vector with r_bar and v_bar
%J2 Perturbation Characteristics
J_2=0.00108;
z=r;
%Setting Time
t=[0 1.577e8]; %5 years in seconds
ACL=@(r)-(mu*r)/norm(r)^3+((3*mu*J_2*r_Earth^2)/(2*norm(r)^5))*(((5*r(3)^2)/(norm(r)^2))-1);
x_dot=@(t,x)[x(4:6);ACL(x(1:3))];
[t,x_t]=ode45(x_dot,t,x_0,odeset('RelTol',1.00e-8)); %Setting relative tolerance to 1e-8
figure('Name','The X/Y vs Time')
plot(t,x_t(:,2:3))
title('X/Y Position vs Time')
xlabel('Time (Hours)')
ylabel('Position (km)')
figure('Name','Top View of the Object''s Orbit')
plot(x_t(:,2),x_t(:,3))
title('Top View of Objects Orbit')
xlabel('Position (km)')
ylabel('Position (km)')
figure('Name','Complete Orbit aroud Earth')
plot3(x_t(:,1),x_t(:,2),x_t(:,3),'-r') %3D plot
axis equal
hold on
[X,Y,Z]=sphere(100);
surf(X*r_Earth,Y*r_Earth,Z*r_Earth,'EdgeColor','none'); %Surface plot, turning off edge color as to not inundate the plot with black
title('Full Object''s Orbit around Earth')
xlabel('Position (km)')
ylabel('Position (km)')
zlabel('Position (km)')
fprintf(['\n\tAs with the results of problem 1, it can be seen by the graphs that the object''s orbit ',...
'is reducing illustrated by the thickness of the line decreasing as the orbit carries on. ',...
'One can also see between the graphs that the J_2 lines are much thicker than the lines ',...
'representing drag. It can also be seen that the J_2 lines are much more prominent at ',...
'circular polar orbit than the drag is. All of the results discussed indicate that the script ',...
'aligns with the results given in lecture.\n'])
toc

### Cite As

Felicia Brimigion (2023). Felicia_Brimigion_Module_6_HW4 (https://www.mathworks.com/matlabcentral/fileexchange/111330-felicia_brimigion_module_6_hw4), MATLAB Central File Exchange. Retrieved .

##### MATLAB Release Compatibility
Created with R2022a
Compatible with any release
##### Platform Compatibility
Windows macOS Linux