% Large-Scale Space Debris/Satellite Dataset Analysis
% Processes and visualizes large TLE dataset from CelesTrak
% Dataset: ~14,000+ active satellites and debris
% Data source: https://celestrak.org

clear all; close all; clc;

fprintf('=== Large-Scale Space Debris Dataset Analysis ===\n');
fprintf('Data source: CelesTrak Active Satellites\n');
fprintf('URL: https://celestrak.org/NORAD/elements/\n\n');

%% Physical Constants
G = 6.67430e-11;
M_earth = 5.972e24;
R_earth = 6.371e6;
mu = G * M_earth;

%% Configuration
tle_file = 'active_satellites_full.txt';
n_visualize = 50;      % Number of satellites to animate (50 for reasonable speed)
n_save_data = 200;     % Number to save full trajectory data (200 for operator learning)
T_sim = 2 * 3600;      % 2 hours simulation
dt = 15;               % 15 seconds time step

if ~exist(tle_file, 'file')
    error('TLE file not found: %s\nPlease ensure the file is in the current directory.', tle_file);
end

%% Load and Parse ALL TLE Data
fprintf('Loading TLE data from: %s\n', tle_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);

n_total = floor(length(tle_lines) / 3);
fprintf('Total lines read: %d\n', length(tle_lines));
fprintf('Total satellites in dataset: %d\n\n', n_total);

%% Parse All Satellites - Extract Orbital Parameters
fprintf('Parsing orbital elements for all %d satellites...\n', n_total);

satellite_data = struct();
satellite_data.names = cell(n_total, 1);
satellite_data.altitude = zeros(n_total, 1);
satellite_data.inclination = zeros(n_total, 1);
satellite_data.eccentricity = zeros(n_total, 1);
satellite_data.period = zeros(n_total, 1);
satellite_data.semi_major_axis = zeros(n_total, 1);

parse_errors = 0;

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

    try
        % Extract name and TLE lines
        satellite_data.names{i} = strtrim(tle_lines{idx});
        tle2 = tle_lines{idx + 2};

        % Parse orbital elements from TLE line 2
        inc_deg = str2double(tle2(9:16));
        ecc_str = ['0.' tle2(27:33)];
        ecc = str2double(ecc_str);
        n_revs_per_day = str2double(tle2(53:63));

        % Calculate orbital parameters
        n = n_revs_per_day * 2 * pi / 86400;
        a = (mu / n^2)^(1/3);
        altitude = a - R_earth;
        period = 2 * pi / n / 60;  % minutes

        % Store
        satellite_data.altitude(i) = altitude;
        satellite_data.inclination(i) = inc_deg;
        satellite_data.eccentricity(i) = ecc;
        satellite_data.period(i) = period;
        satellite_data.semi_major_axis(i) = a;

    catch
        parse_errors = parse_errors + 1;
        satellite_data.altitude(i) = NaN;
        satellite_data.inclination(i) = NaN;
        satellite_data.eccentricity(i) = NaN;
        satellite_data.period(i) = NaN;
        satellite_data.semi_major_axis(i) = NaN;
    end

    if mod(i, 1000) == 0
        fprintf('  Parsed %d/%d satellites...\n', i, n_total);
    end
end

fprintf('Parsing complete! Parse errors: %d\n', parse_errors);

%% Dataset Statistics
fprintf('\n=== Full Dataset Statistics ===\n');

valid_mask = ~isnan(satellite_data.altitude);
alt_km = satellite_data.altitude(valid_mask) / 1000;
inc_deg = satellite_data.inclination(valid_mask);
ecc = satellite_data.eccentricity(valid_mask);

fprintf('Valid satellites: %d / %d\n', sum(valid_mask), n_total);
fprintf('\nAltitude Statistics:\n');
fprintf('  Min:    %.1f km\n', min(alt_km));
fprintf('  Max:    %.1f km\n', max(alt_km));
fprintf('  Mean:   %.1f km\n', mean(alt_km));
fprintf('  Median: %.1f km\n', median(alt_km));

fprintf('\nInclination Statistics:\n');
fprintf('  Min:    %.1f°\n', min(inc_deg));
fprintf('  Max:    %.1f°\n', max(inc_deg));
fprintf('  Mean:   %.1f°\n', mean(inc_deg));

fprintf('\nEccentricity Statistics:\n');
fprintf('  Min:    %.6f\n', min(ecc));
fprintf('  Max:    %.6f\n', max(ecc));
fprintf('  Mean:   %.6f\n', mean(ecc));

% Orbit type classification
LEO_mask = alt_km < 2000;
MEO_mask = alt_km >= 2000 & alt_km < 35786;
GEO_mask = alt_km >= 35786 & alt_km < 36000;
HEO_mask = alt_km >= 36000;

fprintf('\nOrbit Classification:\n');
fprintf('  LEO (< 2000 km):          %d satellites (%.1f%%)\n', sum(LEO_mask), 100*sum(LEO_mask)/sum(valid_mask));
fprintf('  MEO (2000-35786 km):      %d satellites (%.1f%%)\n', sum(MEO_mask), 100*sum(MEO_mask)/sum(valid_mask));
fprintf('  GEO (~35786 km):          %d satellites (%.1f%%)\n', sum(GEO_mask), 100*sum(GEO_mask)/sum(valid_mask));
fprintf('  HEO (> 36000 km):         %d satellites (%.1f%%)\n', sum(HEO_mask), 100*sum(HEO_mask)/sum(valid_mask));

%% Select Satellites for Detailed Analysis

% Strategy: Select diverse sample covering different orbit types
fprintf('\n=== Selecting Satellites for Simulation ===\n');

% Select from different altitude ranges
n_leo = min(30, sum(LEO_mask));
n_meo = min(10, sum(MEO_mask));
n_geo = min(5, sum(GEO_mask));
n_heo = min(5, sum(HEO_mask));

selected_indices = [];

% Select LEO
leo_indices = find(LEO_mask);
if ~isempty(leo_indices)
    selected_indices = [selected_indices; leo_indices(randperm(length(leo_indices), n_leo))];
end

% Select MEO
meo_indices = find(MEO_mask);
if ~isempty(meo_indices)
    selected_indices = [selected_indices; meo_indices(randperm(length(meo_indices), n_meo))];
end

% Select GEO
geo_indices = find(GEO_mask);
if ~isempty(geo_indices)
    selected_indices = [selected_indices; geo_indices(randperm(length(geo_indices), n_geo))];
end

% Select HEO
heo_indices = find(HEO_mask);
if ~isempty(heo_indices)
    selected_indices = [selected_indices; heo_indices(randperm(length(heo_indices), n_heo))];
end

n_selected = length(selected_indices);
fprintf('Selected %d satellites for simulation:\n', n_selected);
fprintf('  LEO: %d, MEO: %d, GEO: %d, HEO: %d\n', n_leo, n_meo, n_geo, n_heo);

%% Convert TLE to Cartesian Coordinates for Selected Satellites

fprintf('\nConverting TLE to Cartesian coordinates...\n');

positions = zeros(n_selected, 3);
velocities = zeros(n_selected, 3);
selected_names = cell(n_selected, 1);

for i = 1:n_selected
    sat_idx = selected_indices(i);
    idx = (sat_idx-1) * 3 + 1;

    selected_names{i} = strtrim(tle_lines{idx});
    tle2 = tle_lines{idx + 2};

    try
        % Parse TLE
        inc_deg = str2double(tle2(9:16));
        RAAN_deg = str2double(tle2(18:25));
        ecc_str = ['0.' tle2(27:33)];
        ecc = str2double(ecc_str);
        omega_deg = str2double(tle2(35:42));
        M_deg = str2double(tle2(44:51));
        n_revs_per_day = str2double(tle2(53:63));

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

        % Calculate semi-major axis
        n = n_revs_per_day * 2 * pi / 86400;
        a = (mu / n^2)^(1/3);

        % Mean to eccentric anomaly
        E = M;
        for iter = 1:10
            E = E - (E - ecc * sin(E) - M) / (1 - ecc * cos(E));
        end

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

        % Position in orbital frame
        r_orb = a * (1 - ecc * cos(E));
        r_peri = [r_orb * cos(nu); r_orb * sin(nu); 0];

        % Velocity in orbital frame
        h = sqrt(mu * a * (1 - ecc^2));
        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
        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];
        R = R3_RAAN * R1_inc * R3_omega;

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

    catch ME
        warning('Failed to convert satellite %d: %s', i, ME.message);
        positions(i, :) = [R_earth + 400e3, 0, 0];
        velocities(i, :) = [0, sqrt(mu/(R_earth+400e3)), 0];
    end
end

fprintf('Conversion complete!\n');

%% Propagate Orbits

t_span = 0:dt:T_sim;
N_steps = length(t_span);

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

trajectories = zeros(n_selected, N_steps, 3);

for i = 1:n_selected
    y0 = [positions(i, :), velocities(i, :)];
    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, 10) == 0
        fprintf('  Propagated %d/%d\n', i, n_selected);
    end
end

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

%% Save Large Dataset

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

fid = fopen(output_file, 'w');
fprintf(fid, '%% Large-Scale Satellite Trajectory Dataset\n');
fprintf(fid, '%% Source: CelesTrak Active Satellites (https://celestrak.org)\n');
fprintf(fid, '%% Total satellites in catalog: %d\n', n_total);
fprintf(fid, '%% Satellites with trajectories: %d\n', n_selected);
fprintf(fid, '%% Simulation duration: %.1f hours\n', T_sim/3600);
fprintf(fid, '%% Time step: %d seconds\n', dt);
fprintf(fid, '%% Columns: Time(s) | Sat_ID | X(m) | Y(m) | Z(m) | Altitude(km) | Name\n\n');

for k = 1:N_steps
    for i = 1:n_selected
        r = norm(squeeze(trajectories(i, k, :)));
        alt_km = (r - R_earth) / 1000;
        fprintf(fid, '%.2f\t%d\t%.6e\t%.6e\t%.6e\t%.2f\t%s\n', ...
            t_span(k), i, ...
            trajectories(i, k, 1), trajectories(i, k, 2), trajectories(i, k, 3), ...
            alt_km, selected_names{i});
    end
end
fclose(fid);

fprintf('Dataset saved! File size: %.1f MB\n', dir(output_file).bytes / 1e6);

%% Create Visualizations

fprintf('\nGenerating visualizations...\n');

% Figure 1: Dataset Overview
figure('Position', [50, 50, 1400, 900]);

% Altitude distribution
subplot(2, 3, 1);
histogram(alt_km, 50, 'FaceColor', [0.3, 0.5, 0.8]);
xlabel('Altitude (km)');
ylabel('Number of Satellites');
title(sprintf('Altitude Distribution (n=%d)', sum(valid_mask)));
grid on;
set(gca, 'YScale', 'log');

% Inclination distribution
subplot(2, 3, 2);
histogram(inc_deg, 50, 'FaceColor', [0.8, 0.4, 0.3]);
xlabel('Inclination (degrees)');
ylabel('Number of Satellites');
title('Inclination Distribution');
grid on;

% Eccentricity distribution
subplot(2, 3, 3);
histogram(ecc, 50, 'FaceColor', [0.3, 0.8, 0.5]);
xlabel('Eccentricity');
ylabel('Number of Satellites');
title('Eccentricity Distribution');
grid on;

% Altitude vs Inclination scatter
subplot(2, 3, 4);
scatter(inc_deg, alt_km, 10, alt_km, 'filled', 'MarkerFaceAlpha', 0.5);
xlabel('Inclination (degrees)');
ylabel('Altitude (km)');
title('Altitude vs Inclination');
colorbar;
colormap(jet);
grid on;
set(gca, 'YScale', 'log');

% Period distribution
subplot(2, 3, 5);
periods = satellite_data.period(valid_mask);
histogram(periods, 50, 'FaceColor', [0.7, 0.3, 0.7]);
xlabel('Orbital Period (minutes)');
ylabel('Number of Satellites');
title('Orbital Period Distribution');
grid on;

% Orbit types pie chart
subplot(2, 3, 6);
labels = {'LEO', 'MEO', 'GEO', 'HEO'};
values = [sum(LEO_mask), sum(MEO_mask), sum(GEO_mask), sum(HEO_mask)];
pie(values, labels);
title('Orbit Type Distribution');

sgtitle(sprintf('CelesTrak Dataset: %d Active Satellites/Debris', n_total), ...
    'FontSize', 14, 'FontWeight', 'bold');

%% Animation

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

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

[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 = jet(n_selected);
frame_skip = 25;
trail_length = 80;

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

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

    % Trajectories
    trail_start = max(1, k - trail_length);
    for i = 1:n_selected
        trail = squeeze(trajectories(i, trail_start:k, :));

        % Handle single point case
        if size(trail, 1) == 3 && size(trail, 2) == 1
            trail = trail';  % Transpose to make it 1x3
        end

        % Only plot trail if we have more than one point
        if size(trail, 1) > 1
            plot3(ax, trail(:,1), trail(:,2), trail(:,3), '-', ...
                'Color', colors(i, :), 'LineWidth', 1);
        end

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

    axis equal; grid on;
    xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');

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

    title(sprintf('Large Dataset: %d/%d Satellites - Time: %.1f min', ...
        n_selected, n_total, 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

fprintf('\n=== Processing Complete ===\n');
fprintf('Total satellites processed: %d\n', n_total);
fprintf('Trajectories computed: %d\n', n_selected);
fprintf('Output file: %s\n', output_file);
fprintf('File size: %.1f MB\n', dir(output_file).bytes / 1e6);
fprintf('\nThis dataset is ready for operator learning applications!\n');

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