Felicia_Brimigion_M​odule_6_HW4

版本 1.0.0 (2.8 KB) 作者: Felicia Brimigion
Homework 4 Module 6
0.0
(0)
5.0 次下载
更新时间 2022/5/7

查看许可证

%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_Earth=6378; %Earth's radius 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

引用格式

Felicia Brimigion (2024). Felicia_Brimigion_Module_6_HW4 (https://www.mathworks.com/matlabcentral/fileexchange/111330-felicia_brimigion_module_6_hw4), MATLAB Central File Exchange. 检索时间: .

MATLAB 版本兼容性
创建方式 R2022a
兼容任何版本
平台兼容性
Windows macOS Linux
标签 添加标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

版本 已发布 发行说明
1.0.0