Counting orbits taken by a satellite

An orbit can be classified as good or an irregular orbit. I think it is better explained with the plot. How to count the number of good orbits and irregular orbits while also keeping a stamp of each orbit that is taken by their starting and ending values. The logic I am using right now is at the end of the code with the two for loops. Note that the indicies of rising and fallilng edges are out of order.

y_foo = readtable("aug33_orbit236_31aug2024_UVA.csv");
y_values = y_foo.Counts;
k = 1;  % Sensitivity parameter
differences = diff(y_values);
% Mean and std deviation
mu = mean(differences);
sigma = std(differences);
% Thresholds
T_high = mu + k * sigma;
T_low = mu - k * sigma;
% Categorizing based on threshold values
states = zeros(1, length(differences));  % Initialize states with zeros (default to steady)
states(differences > T_high) = 1;  % Rising
states(differences < T_low) = -1;   % Falling
% y_values(states == 0 & y_values < 100) = NaN;
% finding the index of the points on rising edge
counts_uva = y_values;
rising_uva_id = [find((diff(counts_uva))>10); length(counts_uva)];
rising_edge_id = rising_uva_id(diff(rising_uva_id)>3400);
% selecting for indices where the difference is more than 100
edge_id = rising_uva_id(diff(rising_uva_id)>100);
falling_uva_id = [find((diff(counts_uva))<-10); length(counts_uva)];
falling_edge_id = falling_uva_id(diff(falling_uva_id)>3400);
% num_normal_orbit = length(rising_edge_id);
% num_irr_orbit = length(falling_edge_id) - length(rising_edge_id);
%trying to remove kink from the rising side
num_rising_edges = length(rising_edge_id);
for i = 1:num_rising_edges
     % ish start points
     current_index = rising_edge_id(i);
     end_index = min(current_index + 800, length(states));
     states(current_index:end_index) = 1;
     prev_index = min(current_index - 175, length(states));
     states(prev_index:current_index) = 1;  % Set to rising state
end
%trying to remove kink from the falling side
num_falling_edges = length(falling_edge_id);
for i = 1:num_falling_edges
     current_index = falling_edge_id(i);
     end_index = min(current_index + 50, length(states));
     states(current_index:end_index) = -1;  % Set to falling state
     prev_index = max(current_index - 800, 1);
     states(prev_index:current_index) = -1;  % Set to falling state
end

% Initialize three vectors of NaNs
N = numel(y_values);
y_rising = NaN(1,N);
y_falling = NaN(1,N);
y_on = NaN(1,N);
% Populate the rising data line
idx = find(states == 1);
y_rising(idx) = y_values(idx);
y_rising(idx+1) = y_values(idx+1);
% Populate the falling data line
idx = find(states == -1);
y_falling(idx) = y_values(idx);
y_falling(idx+1) = y_values(idx+1);
% Populate the steady (on-duty) data line
idx = find(states == 0);
y_on(idx) = y_values(idx);
y_on(idx+1) = y_values(idx+1);

figure
hold on
plot(y_rising,'g','LineWidth',2)  % Rising state in green
plot(y_falling,'r','LineWidth',2)  % Falling state in red
plot(y_on,'b','LineWidth',2)  % Steady state in blue
% Plot original values in dashed lines
plot(y_values, 'k--', 'LineWidth', 1);  % Black dashed line for original values
xlabel('Sample Number');
ylabel('y(t)');
title('Classification of Changes in Square Wave with Kink Handling');





% Count the number of valid orbits
num_correct = 0;

get_rising = [];
% debug = [];
for i = 1:num_rising_edges %indicies of rising edges
     current_inedge = rising_edge_id(i);
     for j = 1:num_falling_edges
         if (abs(current_inedge - falling_edge_id(j))<3400) %&& (current_inedge - falling_edge_id(j)<0)
%             num_weird = num_weird + 1;
%             debug = [debug falling_edge_id(j) current_inedge];
             continue
         end
         if (abs(current_inedge - falling_edge_id(j))>3400)
             num_correct = num_correct + 1;
             get_rising = [get_rising current_inedge];
             break
         end
     end
end

disp(num_correct);
num_weird = num_falling_edges - num_correct;
disp(num_weird)
disp(get_rising);

[William Rose]

@Aditya, Please provide the full plot as a figure. The plot you posted has no axis labels and no legend. Thank you.

Answer 1

[William Rose]

@Aditya,

[Edit: fix typos. No changes to code.]

If you explain your data, you are morelikely to receive useful assistance. You said "it is better explained with the plot". However, it is hard to understand a plot with three colors, no legend, and no axis labels.

The spreadsheet column labels are date/time in column 1, "counts" and "irradiance" in other columns, and three columns that are all zeros. "Counts" and "irradiance" are exactly proportional to one another, so they do not provide independent information.

A standard way to count Earth orbits is to count crossings of the equator, from S to N or from N to S. Do you have such data? I guess not.

The time stamps are at 1 second intervals. However, the last sixty rows have time stamps that are earlier than the preceding 79,000 rows, and there is a gap between the times of the last sixty rows, and all the preceding rows. Therefore I recommend ignoring the last sixty rows.

Are the unexpected dips in counts caused by an orbital disturbance, or by something else, such as occultation by another object, or movement of panels into an unusual position, etc.?

Define start times of "orbits" as times when "counts" make an upward crossing of a threshold. I suspect this an unreliable definition of orbit, so I will use orbit in quotes when defined this way. (I suspect the unexpected dips in counts are not due to a suddent change in orbit, but rather to some other factor.)

To find the number of "orbits", and the time stamps of the start of each orbit, do the following:

% read data
data = readtable("aug33orbit236.csv");
counts = data.Counts(1:end-60);  % skip last 60 rows since time stamps don't fit
N=length(counts);
dt=1;                   % sampling interval (s)
t=(0:N-1)*dt;           % time vector (s)
thresh=12000;           % threshold
k=find(diff(counts>thresh)>0);  % indices of upward threshold crossings
Nk=length(k);           % number of upward threshold crossings
tcross=t(k);            % times of upward crossings

Find the duration of each "orbit":

period=diff(tcross);

Plot counts versus elapsed time, and show the upward threshold crossings:

ycross=repmat(thresh,[1,Nk]);  % y values for plotting threshold crossings
figure
subplot(211), plot(t,counts,'-k',tcross,ycross,'rx'); 
grid on; xlabel('Elasped Time (s)'); ylabel('Counts'); title('+ Threshold Crossings')

Plot the time between threshold crossings, versus elapsed time.

subplot(212), plot(tcross(2:end),period,'rx'); 
grid on; xlabel('Elapsed Time (s)'); ylabel('Period (s)'); title('Time between Crossings')

The second plot above shows that most "orbits" have period~=5500 s, but some "orbits" are shorter. Using your terminology, these correspond to "regular" and "irregular" orbits, respectively. I think that terminology is misleading, but it is your choice, so I will use it.

Count and display the number of regular and irregular orbits:

regOrbitDur=5000;               % minimum duration of "regular orbit" (s)
k=find(period>=regOrbitDur);
Nreg=length(k);
k=find(period<regOrbitDur);
Nirr=length(k);
fprintf('There are %d regular orbits and %d irregular orbits.\n',Nreg,Nirr)

There are 11 regular orbits and 7 irregular orbits.

@Aditya, good luck with your analysis.

[Aditya]

Really sorry about the way I framed the problem. Thank you for helping me out regardless of it.

[William Rose]

@Aditya, no problem.

Answer 2

[Sahas]

Hi @Aditya,

As per my understanding, you are trying to get good and irregular orbits through the data in the attached CSV file and want to count the number of good orbits and irregular orbits while also keeping a stamp of each orbit that is taken by their starting and ending values.

Refer to the sample code below to achieve the required functionality. I have changed the logic of the code where we identify good and irregular objects and assumed the threshold to be "3400" as provided in the code in the question.

% Initialize counters
num_good_orbits = 0;
num_irregular_orbits = 0;

good_orbit_starts = [];
good_orbit_ends = [];

for i = 1:num_rising_edges
    current_rising = rising_edge_id(i);

    % Find the first falling edge that occurs after the current rising edge
    for j = 1:num_falling_edges
        current_falling = falling_edge_id(j);

        % Check if the falling edge is sufficiently far from the rising edge
        if current_falling > current_rising && (current_falling - current_rising) > 3400
            num_good_orbits = num_good_orbits + 1;
            good_orbit_starts = [good_orbit_starts, current_rising];
            good_orbit_ends = [good_orbit_ends, current_falling];
            break;
        end
    end

end

num_irregular_orbits = num_falling_edges - num_good_orbits;

disp('Number of Good Orbits:');
disp(num_good_orbits);
disp('Number of Bad Orbits:');
disp(num_irregular_orbits);
disp('Good Orbit Start Indices:');
disp(good_orbit_starts);
disp('Good Orbit End Indices:');
disp(good_orbit_ends);

I hope this is beneficial!