Problems with FTLE code -- extracting highest FTLE for plot is an empty index. How can this be fixed?
4 次查看(过去 30 天)
显示 更早的评论
I am trying to execute the code below, a finite time lyapunov exponent code with plot generations.
However, at this point:
% Use the indices to get the states at these map returns.
state_lowFTLE = xReturns_running(idx1, :);
state_lowFTLE(1) = state_lowFTLE(1) + 1e-5;
state_highFTLE = xReturns_running(idx2, :);
state_highFTLE(1) = state_highFTLE(1) + 1e-5;
the code provides a value for idx2 (342) that is not possible with xReturns_runnning (a 96x6 matrix). I cannot figure out how to fix this. The error is provided as well.
Index in position 1 exceeds array bounds. Index must
not exceed 96.
Error in FTLE_Main (line 134)
state_highFTLE = xReturns_running(idx2, :);
---------------------------------------------------------------------------------------------------------------
%% THIS CODE PROPAGATES A PERIODIC ORBIT IN CR3BP DYNAMICS
% INITIALIZE WITH:
% 1. guess for the 6 state vector initial conditions (pos,vel) [ndim].
% 2. mass parameter (MU) value for the system [ndim].
% 3. guess for orbit period [ndim].
% THE CODE THEN:
% 1. propagates the reference orbit.
% 2. calculates the FTLE for each time step.
% 3. repeat 1-2 with slight perturbations in order to generate a Poincare
% Map.
%% Clear the MATLAB Workspace.
clear
clc
close all
format compact
format long
%% Constants, initial conditions, and period of orbit.
MU_Earth_Moon = 0.0121505856; % Earth-Moon system mass parameter.
n = 1.0; % Mean-motion of the orbit of the Moon
% around the Earth.
% S/c orbit initial state and time-period.
s0_sc = [ 1.17;
0;
0;
0;
-0.489780292125578;
0];
T = 3.02;
%T = 3.042534324464009;
%% Propagate orbit with corrected conditions.
% Set event function to y == 0, but not terminal so that you can save the
% piercings of the Poincare surface.
ode_options = ...
odeset('Events', @poincare_event_func, 'RelTol', 1e-13, 'AbsTol', 1e-22);
num_orbits = 15;
%num_orbits = 50;
tspan = [0 T*num_orbits];
[t, s_first, t_event, s_event, idx] = ...
ode45(@(t,x) CR3BP(t, x, MU_Earth_Moon, n), tspan, s0_sc, ode_options);
px = s_event(:, 1); % x pos
py = s_event(:, 2); % y pos
pz = s_event(:, 4); % x vel
%% Calculate the Poincare Maps.
f1 = figure(1);
set(0, 'CurrentFigure', f1), plot(px, pz, '.k', 'Markersize', 10),
axis square, xlabel('x [NDU]'), ylabel('xdot [NDU/NTU]'),
title('Poincare Map'), grid on, hold on, set(gca, 'fontsize', 22),
set(gca, 'fontweight', 'bold'), set(gca, 'gridalpha', 0.25),
f1.Color = 'white';
f2 = figure;
% Create running arrays for data saving.
px_running = [];
py_running = [];
pz_running = [];
FTLE_running = [];
xReturns_running = [];
eps = 0.01;
for i = 1 : 3 % Propagate for this many different perturbation distances.
disp(i)
s0_sc = [ s0_sc(1) + i*eps;
0;
0;
0;
s0_sc(5);
0];
num_orbits = 25;
tspan = [0 T*num_orbits];
[t, s, t_event, s_event, idx] = ...
ode45(@(t,x) CR3BP(t, x, MU_Earth_Moon, n), tspan, s0_sc, ode_options);
% Pass Poincare Map piercings to get FTLE at each event.
T_ftleprop = 10; % ndim
[px, py, pz, ftle_spec] = poincare(s_event, T_ftleprop, MU_Earth_Moon, n);
% Append to arrays.
px_running = [px_running; px];
py_running = [py_running; py];
pz_running = [pz_running; pz];
FTLE_running = [FTLE_running ftle_spec];
xReturns_running = [xReturns_running; s_event];
% Plot the Poincare Map piercings.
set(0, 'CurrentFigure', f1), plot(px, pz, '.k', 'Markersize', 10)
axis square, xlabel('x [NDU]'), ylabel('xdot [NDU/NTU]'), grid on,
hold on, set(gca, 'fontsize', 22), set(gca, 'fontweight', 'bold'),
set(gca, 'gridalpha', 0.25), f1.Color = 'white';
% Plot the Poincare Map piercings FTLE.
set(0, 'CurrentFigure', f2), scatter(px, pz, 10, 'c', 'filled'),
axis square, colormap jet, c = colorbar; c.Label.String = 'FTLE Value';
xlabel('x [ndim]'), ylabel('xdot [ndim]'), grid on, hold on
set(gca, 'fontsize', 22), set(gca, 'fontweight', 'bold'),
set(gca, 'gridalpha', 0.25), f2.Color = 'white';
end
hold off
%% Pick the highest and lowest (nonzero) FTLE and propagate the motion.
FTLE_running(FTLE_running == 0) = NaN;
idx1 = find(FTLE_running < 0.45, 1);
idx2 = find(FTLE_running > 0.8, 1);
% Use the indices to get the states at these map returns.
state_lowFTLE = xReturns_running(idx1, :);
state_lowFTLE(1) = state_lowFTLE(1) + 1e-5;
state_highFTLE = xReturns_running(idx2, :);
state_highFTLE(1) = state_highFTLE(1) + 1e-5;
% Propagate these states.
ode_options = odeset('Events', @poincare_event_func, 'RelTol',1e-13, 'AbsTol', 1e-22);
num_orbits = 25;
prop_time = T*num_orbits;
tspan = [0 prop_time];
[t_lowFTLE, s_lowFTLE] = ...
ode45(@(t, x) CR3BP(t, x, MU_Earth_Moon, n), tspan, state_lowFTLE, ode_options);
[t_highFTLE, s_highFTLE] = ...
ode45(@(t, x) CR3BP(t, x, MU_Earth_Moon, n), tspan, state_highFTLE, ode_options);
% Earth-Moon Lagrange Point Locations.
L1x = 0.8369;
L2x = 1.1557;
%% Plot the reference DRO and L1 and L2 points.
figure(3)
plot3((s_first(:,1)), s_first(:,2), s_first(:,4), 'linewidth', 2),
hold on, plot3(L1x, 0, 0, 'o', 'markerfacecolor', 'r', 'markersize', 5),
plot3(L2x, 0, 0, 'o', 'markerfacecolor', 'r', 'markersize', 5),
plot3(1 - MU_Earth_Moon, 0, 0, 'o', 'markerfacecolor', 'k', 'markersize', 10),
xlabel('x [NDU]'), ylabel('y [NDU]'), zlabel('dx/dt [NDU/NTU]'),
legend('Trajectory', 'L1', 'L2', 'Moon'), axis square, grid on,
set(gca, 'fontsize', 22), set(gca, 'fontweight', 'bold'),
set(gca, 'gridalpha', 0.25), f3.Color = 'white';
%% Plot the low and high FTLE returns.
figure(4)
plot3(s_lowFTLE(:,1), s_lowFTLE(:,2), s_lowFTLE(:,4), 'linewidth', 2),
hold on, plot3(L1x, 0, 0, 'o', 'markerfacecolor', 'r', 'markersize', 5),
plot3(L2x, 0, 0, 'o', 'markerfacecolor', 'g', 'markersize', 5),
plot3(1 - MU_Earth_Moon, 0, 0, 'o', 'markerfacecolor', 'k', 'markersize', 10),
xlabel('x [NDU]'), ylabel('y [NDU]'), zlabel('dx/dt [NDU/NTU]'),
axis square, legend('Trajectory', 'L1', 'L2', 'Moon'), grid on,
set(gca, 'fontsize', 22), set(gca, 'fontweight', 'bold'),
set(gca, 'gridalpha', 0.25), f4.Color = 'white';
figure(5)
plot3(s_highFTLE(:,1), s_highFTLE(:,2), s_highFTLE(:,4), 'linewidth', 2),
hold on, plot3(L1x, 0, 0, 'o', 'markerfacecolor', 'r', 'markersize', 5),
plot3(L2x, 0, 0, 'o', 'markerfacecolor', 'g', 'markersize', 5),
plot3(1 - MU_Earth_Moon, 0, 0, 'o', 'markerfacecolor', 'k', 'markersize', 10),
xlabel('x [NDU]'), ylabel('y [NDU]'), zlabel('dx/dt [NDU/NTU]'),
axis square, legend('Trajectory', 'L1', 'L2', 'Moon'),
grid on, set(gca, 'fontsize', 22), set(gca, 'fontweight', 'bold'),
set(gca, 'gridalpha', 0.25), f5.Color = 'white';
1 个评论
dpb
2022-8-26
FTLE_running = [FTLE_running ftle_spec];
xReturns_running = [xReturns_running; s_event];
...
%loop containing above ends and then
...
idx1 = find(FTLE_running < 0.45, 1);
idx2 = find(FTLE_running > 0.8, 1);
...
state_highFTLE = xReturns_running(idx2, :);
...
There's insufficient code provided to see what is done and insufficient comments to have a klew about what is intended to be happening here, but in the above the FTLE_running vector is built and used to find a particular set of locations. BUT, the xReturns_running array is (apparently) being built as an array by catenating whatever s_event is vertically while FTLE_running is (also apparently) a row vector. It may be possible/true that the total number of points in the two is the same, but the number of rows in the 2D array wouldn't match the linear position in a single vector.
Set breakpoint and study your code logic caretully in this area; it clearly isn't doing what you really intended; I'm guessing there's a storage orientation issue going on...
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Dates and Time 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!