# Untitled

unknown
matlab
7 months ago
3.5 kB
11
Indexable
Never
```clc;
clear;
close all;
G = 6.67430e-20; % km^3/(kg * s^2)
m1 = 5.97219e24; % kg
m2 = 62; %Mass of the satellite kg
mu = G * m1; % km^3/s^2
RE = 6378.12; % km

Cd = 2.2; % Example value for drag coefficient
A = 2.1840e-06; % Example value for cross-sectional area in square meters
r0 = [-303.7171769308391 -2451.236755737989 -6680.159742109596]; % km
v0 = [7.028651283303245 2.34824252032872 -1.194904065231528 ]; % km/s
Y0 = [r0 v0];

t0 = 0; % seconds
tf = 86400; % seconds, period of one orbit
time_step = 60;  % specify the time step in seconds
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-9, 'MaxStep', time_step);
[t, Y] = ode23s(@(t, Y) relative_motion_drag_J2(t, Y, mu, Cd, A, m2), t0:time_step:tf, Y0, options);

rvec = Y(:, 1:3);
vvec = Y(:, 4:6);

rmag = vecnorm(rvec, 2, 2);
altitude = rmag - RE;

% Plot altitude as a function of time
figure;
plot(t/86400, altitude);
xlabel('Time (s)');
ylabel('Altitude (km)');
title('Altitude vs. Time');

% Create a new figure for the 3D orbit
figure;
[xx, yy, zz] = sphere(200);
surf(RE * xx, RE * yy, RE * zz)
colormap(light_blue)
clim([-RE / 100 RE / 100])

% Draw and label the X, Y, and Z axes
line([0 2 * RE], [0 0], [0 0]);
text(2 * RE + RE / 10, 0, 0, 'x')
line([0 0], [0 2 * RE], [0 0]);
text(0, 2 * RE + RE / 10, 0, 'y')
line([0 0], [0 0], [0 2 * RE]);
text(0, 0, 2 * RE + RE / 10, 'z')

hold on
plot3(rvec(:, 1), rvec(:, 2), rvec(:, 3), 'k')
plot3(rvec(1, 1), rvec(1, 2), rvec(1, 3), 'go'); % Starting point
plot3(rvec(end, 1), rvec(end, 2), rvec(end, 3), 'ro'); % Ending point

% Specify some properties of the graph
grid on
axis equal
xlabel('km')
ylabel('km')
zlabel('km')
title('3D Orbit');

function map = light_blue
r = 0.8;
g = r;
b = 1.0;
map = [r g b; 0 0 0; r g b];
end

function Ydot = relative_motion_drag_J2(t, Y, mu, Cd, A, m)
% Earth's J2 coefficient and equatorial radius
J2 = 1.08263e-3;
RE = 6378.12; % km
Omega = 7.2921*10^-5;
Omega = [0.0, 0.0, Omega]'; %
% Extract position and velocity vectors
rvector = Y(1:3);
vvector = Y(4:6);

% Standard gravitational acceleration
r = norm(rvector);
a_gravity = -mu * rvector / r^3;

% J2 perturbation
z2 = rvector(3)^2;
r2 = r^2;
factor = 1.5 * J2 * mu * RE^2 / r2^2;
ax_J2 = rvector(1) / r * (5 * z2 / r2 - 1) * factor;
ay_J2 = rvector(2) / r * (5 * z2 / r2 - 1) * factor;
az_J2 = rvector(3) / r * (5 * z2 / r2 - 3) * factor;
a_J2 = [ax_J2; ay_J2; az_J2];

% Calculate the relative velocity
v_icrf = Y(4:6);
r_icrf = Y(1:3);
if isrow(r_icrf)
r_icrf = r_icrf';
end
if isrow(v_icrf)
v_icrf = v_icrf';
end

% Calculate the relative velocity
v_rel = v_icrf - cross(Omega, r_icrf);

% Calculate the drag acceleration
v_rel_mag = norm(v_rel);
rho = atmospheric_density(norm(r_icrf) - RE); % call to your atmospheric density function
a_drag = -0.5 * Cd * A / m *  rho * v_rel_mag^2 * (v_rel / v_rel_mag)

% Total acceleration
a_total = a_gravity + a_J2 + a_drag;

% Output derivative of state vector
Ydot = [vvector; a_total];
end

function rho = atmospheric_density(altitude)
% Simple exponential atmosphere model (or use a more complex model)
rho0 = 3.614*10^-14; % kg/m^3 at sea level
H = 8500; % scale height in meters
rho = rho0 * exp(-((altitude-700) / 88.667));
end

```