Untitled

%ISS NORAD ID (25544)
%2019/08/13
%2023/08/12
clc;
clear;
close all;
%% real time
% List of TLE file names
%tleFiles = ["satellite1.tle", "satellite2.tle", "satellite3.tle"];
for z = 1:7758
    fileName = sprintf('longlong%d.tle',z);  % Generate the file name
    tleFiles{z} = fileName;                 % Read the file contents
end
% Define the time range for propagsation
startTime = datetime(2019,8,13);
stopTime = datetime(2023,8,12);
sampleTime =86400;
timeRange = startTime:seconds(16261):stopTime;
% Create the satellite scenario
sc = satelliteScenario(startTime, stopTime, sampleTime);

% Loop through each TLE file    
for i = 1:7758
    tleFile = tleFiles{i};
    tleFilePath = tleFile;
    [epochYearArray, month, day] = extractTLEData(tleFilePath);
    
    try
        % Create the satellite object using the TLE file
        sat = satellite(sc, tleFile);
        
        % Calculate the year for the specified date
        year = 2000 + epochYearArray;
        
        % Specify the time for which you want to get position and velocity
        time = datetime(year, month, day, 12, 0, 0);

        % Compute position and velocity in ECEF coordinate frame
        [position, velocity] = states(sat, time, "CoordinateFrame", "ecef");

        % Display the results
        disp("TLE File: " + tleFile);
        position_realn(i) = (norm(position') / 1000) - 6371;
        velocity_realn(i) = norm(velocity') / 1000;
        disp(time);
        disp("Position (ECEF) [km]: " + position_realn(i));
        disp("Velocity (ECEF) [km/s]: " + velocity_realn(i));
        disp("-------------------------------");

    catch exception
        fprintf('Inconsistent data in iteration %d, skipped.\n', i);
        disp("Error Message: " + exception.message);
    end
end
%play(sc)
%% SGP4

startTime = datetime(2019,8,13);
stopTime = datetime(2023,8,12);
sampleTime =86400;
sc1 = satelliteScenario(startTime,stopTime,sampleTime);

tleFile = "long1.tle";


satSGP4 = satellite(sc1,tleFile, ...
    "Name","satSGP4", ...
    "OrbitPropagator","sgp4");

satSDP4 = satellite(sc1,tleFile, ...
    "Name","satSDP4", ...
    "OrbitPropagator","sdp4");

[positionSGP4,velocitySGP4] = states(satSGP4,"CoordinateFrame", "ecef");
[positionSDP4,velocitySDP4] = states(satSDP4,"CoordinateFrame", "ecef");
for x=1:7758
    positionSGP4_n(x) = (norm(positionSGP4(:,x)')/1000)-6371;
    velocitySGP4_n(x) = norm(velocitySGP4(:,x)')/1000;
        
    
    positionSDP4_n(x) = (norm(positionSDP4(:,x)')/1000)-6371;
    velocitySDP4_n(x) = norm(velocitySDP4(:,x)')/1000;
end 
%v = satelliteScenarioViewer(sc1);
%play(sc1)
%% RelativePosition
re_pos_sgp4=abs(position_realn-positionSGP4_n);
re_pos_sgp4_mean = mean(re_pos_sgp4);
    
re_pos_sdp4=abs(position_realn-positionSDP4_n);
re_pos_sdp4_mean = mean(re_pos_sdp4);

test_r_p=(positionSDP4_n-positionSGP4_n);

figure;

subplot(2,1,1)
plot(timeRange,re_pos_sgp4,'r-',timeRange,re_pos_sgp4_mean,'b->')
xlim("auto");
ylim("auto");
ylabel('km');
title('RelativePosition:SGP4');

subplot(2,1,2)
plot(timeRange,re_pos_sdp4,'r-',timeRange,re_pos_sdp4_mean,'b->')
xlim("auto")
ylim("auto");
ylabel('km');
title('RelativePosition:SDP4');