Untitled

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])
shading interp

% 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