% Space Debris/Satellite Orbital Simulation
% Simulates orbital dynamics of debris around Earth under gravitational force
% with 3D animation and trajectory data output

clear all; close all; clc;

fprintf('=== Space Debris Orbital Simulation ===\n');

%% Physical Constants
G = 6.67430e-11;           % Gravitational constant [m^3 kg^-1 s^-2]
M_earth = 5.972e24;        % Mass of Earth [kg]
R_earth = 6.371e6;         % Radius of Earth [m]
mu = G * M_earth;          % Standard gravitational parameter [m^3/s^2]

fprintf('Gravitational parameter: %.4e m^3/s^2\n', mu);
fprintf('Earth radius: %.1f km\n', R_earth/1000);

%% Simulation Parameters
n_debris = 10;             % Number of debris/satellites
T_sim = 2 * 3600;          % Simulation time [seconds] (2 hours)
dt = 10;                   % Time step [seconds]
t_span = 0:dt:T_sim;
N_steps = length(t_span);

fprintf('\nSimulation settings:\n');
fprintf('Number of objects: %d\n', n_debris);
fprintf('Simulation duration: %.1f hours\n', T_sim/3600);
fprintf('Time step: %d seconds\n', dt);
fprintf('Total steps: %d\n', N_steps);

%% Generate Random Initial Conditions for Debris
% We'll create orbits with random:
% - Orbital altitude (LEO to MEO range)
% - Orbital inclination
% - Orientation (longitude of ascending node, argument of periapsis)
% - Eccentricity (mostly circular with small variations)

fprintf('\nGenerating initial conditions...\n');

rng(42);  % For reproducibility

% Storage for initial conditions
positions = zeros(n_debris, 3);  % [x, y, z] in meters
velocities = zeros(n_debris, 3); % [vx, vy, vz] in m/s

% Define orbital parameter ranges
alt_min = 300e3;           % Minimum altitude: 300 km (LEO)
alt_max = 2000e3;          % Maximum altitude: 2000 km (MEO)
ecc_max = 0.3;             % Maximum eccentricity

for i = 1:n_debris
    % Random orbital altitude
    altitude = alt_min + (alt_max - alt_min) * rand();
    r_orbit = R_earth + altitude;

    % Random eccentricity (mostly circular orbits)
    ecc = ecc_max * rand();

    % Random orbital inclination [0, 180] degrees
    inc = pi * rand();

    % Random longitude of ascending node [0, 360] degrees
    RAAN = 2 * pi * rand();

    % Random argument of periapsis [0, 360] degrees
    omega = 2 * pi * rand();

    % Random true anomaly (position along orbit) [0, 360] degrees
    nu = 2 * pi * rand();

    % Compute semi-major axis from mean radius and eccentricity
    a = r_orbit / (1 + ecc * cos(nu));

    % Compute position and velocity in orbital plane
    % Using orbital elements to Cartesian conversion
    r_orb = a * (1 - ecc^2) / (1 + ecc * cos(nu));

    % Position in orbital plane
    r_peri = [r_orb * cos(nu); r_orb * sin(nu); 0];

    % Velocity in orbital plane
    h = sqrt(mu * a * (1 - ecc^2));  % Specific angular momentum
    v_r = mu / h * ecc * sin(nu);     % Radial velocity
    v_theta = h / r_orb;              % Tangential velocity

    v_peri = [v_r * cos(nu) - v_theta * sin(nu);
              v_r * sin(nu) + v_theta * cos(nu);
              0];

    % Rotation matrices to transform to inertial frame
    R3_omega = [cos(omega), -sin(omega), 0;
                sin(omega), cos(omega), 0;
                0, 0, 1];

    R1_inc = [1, 0, 0;
              0, cos(inc), -sin(inc);
              0, sin(inc), cos(inc)];

    R3_RAAN = [cos(RAAN), -sin(RAAN), 0;
               sin(RAAN), cos(RAAN), 0;
               0, 0, 1];

    % Combined rotation matrix
    R = R3_RAAN * R1_inc * R3_omega;

    % Transform to inertial frame
    positions(i, :) = (R * r_peri)';
    velocities(i, :) = (R * v_peri)';

    % Display info
    fprintf('Debris %d: alt=%.0f km, ecc=%.3f, inc=%.1f deg\n', ...
        i, altitude/1000, ecc, rad2deg(inc));
end

%% Numerical Integration of Equations of Motion
% Using ode45 for each debris

fprintf('\nIntegrating orbits...\n');

% Storage for all trajectories
trajectories = zeros(n_debris, N_steps, 3);  % [debris, time, xyz]

for i = 1:n_debris
    % Initial state vector [x, y, z, vx, vy, vz]
    y0 = [positions(i, :), velocities(i, :)];

    % Integrate equations of motion
    options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
    [~, y] = ode45(@(t, y) orbital_dynamics(t, y, mu), t_span, y0, options);

    % Store trajectory
    trajectories(i, :, :) = y(:, 1:3);

    if mod(i, 2) == 0
        fprintf('  Computed orbit for debris %d/%d\n', i, n_debris);
    end
end

fprintf('Orbital integration complete!\n');

%% Save Data to File
filename = 'space_debris_trajectories.dat';
fprintf('\nSaving trajectory data to %s...\n', filename);

% Open file for writing
fid = fopen(filename, 'w');

% Write header
fprintf(fid, '%% Space Debris Trajectory Data\n');
fprintf(fid, '%% Columns: Time(s), Debris_ID, X(m), Y(m), Z(m)\n');
fprintf(fid, '%% Number of debris: %d\n', n_debris);
fprintf(fid, '%% Simulation duration: %.1f hours\n', T_sim/3600);
fprintf(fid, '%% Time step: %d seconds\n\n', dt);

% Write data
for k = 1:N_steps
    for i = 1:n_debris
        fprintf(fid, '%.2f\t%d\t%.6e\t%.6e\t%.6e\n', ...
            t_span(k), i, ...
            trajectories(i, k, 1), ...
            trajectories(i, k, 2), ...
            trajectories(i, k, 3));
    end
end

fclose(fid);
fprintf('Data saved successfully!\n');

%% Create Animation

fprintf('\nCreating animation...\n');

fig = figure('Position', [100, 100, 1200, 900]);
ax = axes('Position', [0.05, 0.05, 0.9, 0.9]);

% Create Earth sphere
[x_sphere, y_sphere, z_sphere] = sphere(50);
x_sphere = x_sphere * R_earth;
y_sphere = y_sphere * R_earth;
z_sphere = z_sphere * R_earth;

% Color for debris trails
colors = jet(n_debris);

% Animation parameters
frame_skip = 20;  % Display every 20th frame for speed
trail_length = 100;  % Number of points in trail

for k = 1:frame_skip:N_steps
    cla(ax);
    hold(ax, 'on');

    % Draw Earth
    surf(ax, x_sphere, y_sphere, z_sphere, ...
        'FaceColor', [0.2, 0.4, 0.8], ...
        'EdgeColor', 'none', ...
        'FaceAlpha', 0.9);

    % Draw trajectories (trails)
    trail_start = max(1, k - trail_length);

    for i = 1:n_debris
        % Trail
        trail_x = squeeze(trajectories(i, trail_start:k, 1));
        trail_y = squeeze(trajectories(i, trail_start:k, 2));
        trail_z = squeeze(trajectories(i, trail_start:k, 3));

        plot3(ax, trail_x, trail_y, trail_z, '-', ...
            'Color', colors(i, :), 'LineWidth', 1.5, 'HandleVisibility', 'off');

        % Current position (marker)
        plot3(ax, trajectories(i, k, 1), ...
                  trajectories(i, k, 2), ...
                  trajectories(i, k, 3), ...
            'o', 'MarkerSize', 8, 'MarkerFaceColor', colors(i, :), ...
            'MarkerEdgeColor', 'k', 'LineWidth', 1, ...
            'DisplayName', sprintf('Debris %d', i));
    end

    % Settings
    axis equal;
    grid on;
    xlabel('X (m)', 'FontSize', 11);
    ylabel('Y (m)', 'FontSize', 11);
    zlabel('Z (m)', 'FontSize', 11);

    % Set axis limits (slightly larger than max orbit)
    max_r = max(max(max(abs(trajectories))));
    lim = max_r * 1.1;
    xlim([-lim, lim]);
    ylim([-lim, lim]);
    zlim([-lim, lim]);

    title(sprintf('Space Debris Orbital Simulation - Time: %.1f min (%.1f%%)', ...
        t_span(k)/60, 100*k/N_steps), 'FontSize', 13);

    % View angle
    view(45, 20);

    % Lighting
    lighting gouraud;
    camlight('headlight');

    hold(ax, 'off');

    drawnow;

    % Display progress
    if mod(k, 100) == 0 || k == N_steps
        fprintf('  Animation progress: %.1f%%\n', 100*k/N_steps);
    end
end

fprintf('\nAnimation complete!\n');

%% Summary Plot: All Orbits
figure('Position', [150, 150, 1200, 800]);

% 3D view
subplot(1, 2, 1);
hold on;

% Draw Earth
surf(x_sphere, y_sphere, z_sphere, ...
    'FaceColor', [0.2, 0.4, 0.8], ...
    'EdgeColor', 'none', ...
    'FaceAlpha', 0.7);

% Draw all complete orbits
for i = 1:n_debris
    plot3(trajectories(i, :, 1), ...
          trajectories(i, :, 2), ...
          trajectories(i, :, 3), ...
        '-', 'Color', colors(i, :), 'LineWidth', 1.5);
end

axis equal;
grid on;
xlabel('X (m)', 'FontSize', 11);
ylabel('Y (m)', 'FontSize', 11);
zlabel('Z (m)', 'FontSize', 11);
title('Complete Orbital Trajectories', 'FontSize', 12);
view(45, 20);
lighting gouraud;
camlight('headlight');
hold off;

% Altitude vs Time
subplot(1, 2, 2);
hold on;
for i = 1:n_debris
    % Compute altitude over time
    r = sqrt(trajectories(i, :, 1).^2 + ...
             trajectories(i, :, 2).^2 + ...
             trajectories(i, :, 3).^2);
    altitude_km = (r - R_earth) / 1000;

    plot(t_span / 60, altitude_km, '-', 'Color', colors(i, :), ...
        'LineWidth', 1.5, 'DisplayName', sprintf('Debris %d', i));
end
grid on;
xlabel('Time (minutes)', 'FontSize', 11);
ylabel('Altitude (km)', 'FontSize', 11);
title('Altitude Evolution', 'FontSize', 12);
legend('Location', 'eastoutside', 'FontSize', 8);
hold off;

sgtitle('Space Debris Orbital Dynamics Summary', 'FontSize', 14, 'FontWeight', 'bold');

%% Final Statistics
fprintf('\n=== Simulation Summary ===\n');

for i = 1:n_debris
    % Compute orbital statistics
    r = sqrt(trajectories(i, :, 1).^2 + ...
             trajectories(i, :, 2).^2 + ...
             trajectories(i, :, 3).^2);

    r_min = min(r);
    r_max = max(r);
    r_mean = mean(r);

    alt_min_km = (r_min - R_earth) / 1000;
    alt_max_km = (r_max - R_earth) / 1000;
    alt_mean_km = (r_mean - R_earth) / 1000;

    % Estimate orbital period (time for one complete orbit)
    a_est = r_mean;
    T_orbit = 2 * pi * sqrt(a_est^3 / mu) / 60;  % minutes

    fprintf('Debris %d: Alt = %.0f-%.0f km (mean %.0f km), Period ≈ %.1f min\n', ...
        i, alt_min_km, alt_max_km, alt_mean_km, T_orbit);
end

fprintf('\nSimulation complete!\n');
fprintf('Data file: %s\n', filename);

%% Helper Function: Orbital Dynamics
function dydt = orbital_dynamics(t, y, mu)
    % Equations of motion for orbital dynamics
    % y = [x, y, z, vx, vy, vz]'
    % dydt = [vx, vy, vz, ax, ay, az]'

    % Extract position and velocity
    r = y(1:3);
    v = y(4:6);

    % Compute distance from Earth center
    r_mag = norm(r);

    % Gravitational acceleration: a = -mu * r / r^3
    a = -mu / r_mag^3 * r;

    % Assemble derivative
    dydt = [v; a];
end
