% Space Debris/Satellite Visualization with Real TLE Data
% Loads and animates real satellite trajectories from TLE (Two-Line Element) data
% Data source: CelesTrak (https://celestrak.org)

clear all; close all; clc;

fprintf('=== Real Space Debris/Satellite Visualization ===\n');
fprintf('Data source: CelesTrak - https://celestrak.org\n\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]

%% Read TLE Data from File
tle_file = 'iss_tle_data.txt';

if ~exist(tle_file, 'file')
    error('TLE data file not found: %s', tle_file);
end

fprintf('Loading TLE data from: %s\n', tle_file);

% Read the file
fid = fopen(tle_file, 'r');
tle_lines = {};
while ~feof(fid)
    line = fgetl(fid);
    if ischar(line)
        tle_lines{end+1} = line;
    end
end
fclose(fid);

fprintf('Read %d lines from TLE file\n', length(tle_lines));

%% Parse TLE Data
% TLE format: 3 lines per object
% Line 0: Name
% Line 1: TLE line 1 (orbital elements part 1)
% Line 2: TLE line 2 (orbital elements part 2)

n_objects = floor(length(tle_lines) / 3);
fprintf('Found %d satellite/debris objects\n\n', n_objects);

% Select first 10 objects for visualization
n_display = min(10, n_objects);

satellite_names = cell(n_display, 1);
positions = zeros(n_display, 3);
velocities = zeros(n_display, 3);

fprintf('Parsing TLE data for %d objects:\n', n_display);

for i = 1:n_display
    idx = (i-1) * 3 + 1;

    % Extract name
    satellite_names{i} = strtrim(tle_lines{idx});

    % Extract TLE lines
    tle1 = tle_lines{idx + 1};
    tle2 = tle_lines{idx + 2};

    % Parse TLE elements
    try
        % From TLE line 2: inclination, RAAN, eccentricity, arg of perigee, mean anomaly, mean motion
        inc_deg = str2double(tle2(9:16));        % Inclination [deg]
        RAAN_deg = str2double(tle2(18:25));      % Right Ascension of Ascending Node [deg]
        ecc_str = ['0.' tle2(27:33)];            % Eccentricity
        ecc = str2double(ecc_str);
        omega_deg = str2double(tle2(35:42));     % Argument of perigee [deg]
        M_deg = str2double(tle2(44:51));         % Mean anomaly [deg]
        n_revs_per_day = str2double(tle2(53:63)); % Mean motion [revs/day]

        % Convert to radians
        inc = deg2rad(inc_deg);
        RAAN = deg2rad(RAAN_deg);
        omega = deg2rad(omega_deg);
        M = deg2rad(M_deg);

        % Calculate orbital parameters
        n = n_revs_per_day * 2 * pi / 86400;  % Mean motion [rad/s]
        a = (mu / n^2)^(1/3);                  % Semi-major axis [m]

        % Convert mean anomaly to eccentric anomaly (Newton's method)
        E = M;  % Initial guess
        for iter = 1:10
            E = E - (E - ecc * sin(E) - M) / (1 - ecc * cos(E));
        end

        % Convert to true anomaly
        nu = 2 * atan2(sqrt(1 + ecc) * sin(E/2), sqrt(1 - ecc) * cos(E/2));

        % Compute position and velocity in orbital plane
        r_orb = a * (1 - ecc * cos(E));

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

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

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

        % Rotation matrices 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
        R = R3_RAAN * R1_inc * R3_omega;

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

        % Calculate altitude
        altitude_km = (r_orb - R_earth) / 1000;

        fprintf('  %2d. %-20s: alt=%.0f km, inc=%.1f°, ecc=%.4f\n', ...
            i, satellite_names{i}, altitude_km, inc_deg, ecc);

    catch ME
        fprintf('  %2d. %-20s: ERROR parsing TLE\n', i, satellite_names{i});
        warning('Error parsing TLE for object %d: %s', i, ME.message);
        % Set to default circular orbit
        positions(i, :) = [R_earth + 400e3, 0, 0];
        v_circ = sqrt(mu / (R_earth + 400e3));
        velocities(i, :) = [0, v_circ, 0];
    end
end

%% Numerical Integration of Orbits

T_sim = 2 * 3600;          % 2 hours simulation
dt = 10;                   % 10 seconds time step
t_span = 0:dt:T_sim;
N_steps = length(t_span);

fprintf('\nPropagating orbits...\n');
fprintf('  Duration: %.1f hours\n', T_sim/3600);
fprintf('  Time steps: %d\n', N_steps);

% Storage for trajectories
trajectories = zeros(n_display, N_steps, 3);

for i = 1:n_display
    % Initial state
    y0 = [positions(i, :), velocities(i, :)];

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

    trajectories(i, :, :) = y(:, 1:3);

    if mod(i, 3) == 0 || i == n_display
        fprintf('  Propagated %d/%d orbits\n', i, n_display);
    end
end

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

%% Save Trajectories to Data File

output_file = 'real_satellite_trajectories.dat';
fprintf('\nSaving trajectory data to: %s\n', output_file);

fid = fopen(output_file, 'w');
fprintf(fid, '%% Real Satellite/Debris Trajectory Data\n');
fprintf(fid, '%% Data source: CelesTrak - https://celestrak.org\n');
fprintf(fid, '%% Downloaded: %s\n', datestr(now));
fprintf(fid, '%% Columns: Time(s), Sat_ID, Satellite_Name, X(m), Y(m), Z(m)\n');
fprintf(fid, '%% Number of objects: %d\n', n_display);
fprintf(fid, '%% Simulation duration: %.1f hours\n\n', T_sim/3600);

for k = 1:N_steps
    for i = 1:n_display
        fprintf(fid, '%.2f\t%d\t%s\t%.6e\t%.6e\t%.6e\n', ...
            t_span(k), i, satellite_names{i}, ...
            trajectories(i, k, 1), ...
            trajectories(i, k, 2), ...
            trajectories(i, k, 3));
    end
end

fclose(fid);
fprintf('Data saved!\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]);

% 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;

% Colors
colors = jet(n_display);

% Animation parameters
frame_skip = 20;
trail_length = 100;

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 orbits
    trail_start = max(1, k - trail_length);

    for i = 1:n_display
        % 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);

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

        % Label for key satellites
        if i <= 3  % Label first 3
            text(trajectories(i, k, 1), ...
                 trajectories(i, k, 2), ...
                 trajectories(i, k, 3), ...
                 ['  ' satellite_names{i}], ...
                 'Color', colors(i, :), 'FontSize', 9, 'FontWeight', 'bold');
        end
    end

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

    % Axis limits
    max_r = max(max(max(abs(trajectories))));
    lim = max_r * 1.1;
    xlim([-lim, lim]);
    ylim([-lim, lim]);
    zlim([-lim, lim]);

    title(sprintf('Real Satellite/Debris Orbits - Time: %.1f min', t_span(k)/60), ...
        'FontSize', 13);

    view(45, 20);
    lighting gouraud;
    camlight('headlight');

    hold(ax, 'off');
    drawnow;

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

fprintf('Animation complete!\n');

%% Summary Plots

figure('Position', [150, 150, 1400, 800]);

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

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

for i = 1:n_display
    plot3(trajectories(i, :, 1), ...
          trajectories(i, :, 2), ...
          trajectories(i, :, 3), ...
        '-', 'Color', colors(i, :), 'LineWidth', 1.5, ...
        'DisplayName', satellite_names{i});
end

axis equal;
grid on;
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Z (m)');
title('Complete Orbital Trajectories (Real Data)');
view(45, 20);
lighting gouraud;
camlight('headlight');
legend('Location', 'eastoutside', 'FontSize', 8);
hold off;

% Altitude evolution
subplot(1, 2, 2);
hold on;
for i = 1:n_display
    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', satellite_names{i});
end
grid on;
xlabel('Time (minutes)');
ylabel('Altitude (km)');
title('Altitude vs Time');
legend('Location', 'eastoutside', 'FontSize', 8);
hold off;

sgtitle('Real Satellite/Debris Data from CelesTrak', 'FontSize', 14, 'FontWeight', 'bold');

%% Display Statistics

fprintf('\n=== Orbital Statistics ===\n');

for i = 1:n_display
    r = sqrt(trajectories(i, :, 1).^2 + ...
             trajectories(i, :, 2).^2 + ...
             trajectories(i, :, 3).^2);

    alt_km = (r - R_earth) / 1000;
    alt_min = min(alt_km);
    alt_max = max(alt_km);
    alt_mean = mean(alt_km);

    % Estimate period
    T_period = 2 * pi * sqrt(mean(r)^3 / mu) / 60;  % minutes

    fprintf('%-20s: Alt %.0f-%.0f km (mean %.0f km), Period %.1f min\n', ...
        satellite_names{i}, alt_min, alt_max, alt_mean, T_period);
end

fprintf('\nData files created:\n');
fprintf('  - TLE input: %s\n', tle_file);
fprintf('  - Trajectory output: %s\n', output_file);
fprintf('\nVisualization complete!\n');

%% Helper Function
function dydt = orbital_dynamics(t, y, mu)
    % Orbital dynamics under Earth's gravity
    r = y(1:3);
    v = y(4:6);
    r_mag = norm(r);
    a = -mu / r_mag^3 * r;
    dydt = [v; a];
end
