Counting orbits taken by a satellite

20 次查看（过去 30 天）
Aditya
Aditya 2024-10-22
回答： Sahas 2024-10-23，6:01
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);
  1 个评论
William Rose
William Rose 2024-10-23，1:16
@Aditya, Please provide the full plot as a figure. The plot you posted has no axis labels and no legend. Thank you.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

William Rose
William Rose 2024-10-23，2:47
编辑：William Rose 2024-10-23，5:10
Ran in:
@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.
  2 个评论
Aditya
Aditya 2024-10-23，4:21
Really sorry about the way I framed the problem. Thank you for helping me out regardless of it.

请先登录，再进行评论。

Sahas
Sahas 2024-10-23，6:01
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!

类别

在 Help Center 和 File Exchange 中查找有关 View and Analyze Simulation Results 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by