I need help with my MATLAB code. I wanted to numerically integrate MRPs. When there was a singularity in MRP, my code would stop it using an event function. Now, I want to switch to the MRP shadow set and run the ode45 again with the updated initial conditions. How do I do that? clc; clear all; w0 = [0; 0; 0.3]; MRP0 = [0; 0; -0.198912367379658]; tspan = [0, 100]; % Integrate the Euler equations using ode45 with an event function options = odeset('RelTol', 1e-10, 'AbsTol', 1e-10, 'Events', @singularityEvent); [t, y, te, ye, ie] = ode45(@KDE_MRP, tspan, [MRP0; w0], options); % Extract the Euler parameters and angular velocities MRP = y(:, 1:3); w = y(:, 4:6); % plotting MRP vs time plot(t,MRP, '-') xlabel('time (s)') ylabel('MRP') legend('MRP1', 'MRP2', 'MRP3') function [value, isterminal, direction] = singularityEvent(t, MRPw) % Event function to detect when the norm of MRP approaches 2 value = norm(MRPw(1:3)) - 1; isterminal = 1; % Stop the integration when the event is detected direction = 0; % Detect events in both directions end function dMRPwdt = KDE_MRP(t, MRPw) I = [0.3; 0.2; 0.4]; L = [0; 0; 0]; MRP = MRPw(1:3); w = MRPw(4:6); dMRPdt = zeros(3, 1); dwdt = zeros(3, 1); C11 = 1 - dot(MRP, MRP) + 2 * MRP(1)^2; C12 = 2 * (MRP(1) * MRP(2) - MRP(3)); C13 = 2 * (MRP(1) * MRP(3) + MRP(2)); C21 = 2 * (MRP(2) * MRP(1) + MRP(3)); C22 = 1 - dot(MRP, MRP) + 2 * MRP(2)^2; C23 = 2 * (MRP(2) * MRP(3) - MRP(1)); C31 = 2 * (MRP(3) * MRP(1) - MRP(2)); C32 = 2 * (MRP(3) * MRP(2) + MRP(1)); C33 = 1 - dot(MRP, MRP) + 2 * MRP(3)^2; C = [C11 C12 C13; C21 C22 C23; C31 C32 C33]; dMRPdt = 0.25 * C * w; dwdt(1) = (-(I(3) - I(2)) * w(2) * w(3) + L(1)) / I(1); dwdt(2) = (-(I(1) - I(3)) * w(3) * w(1) + L(2)) / I(2); dwdt(3) = (-(I(2) - I(1)) * w(1) * w(2) + L(3)) / I(3); % Combine the time derivatives into a single vector dMRPwdt = [dMRPdt; dwdt]; end
I need help with my MATLAB code. I wanted to numerically integrate MRPs. When there was a singularity in MRP, my code would stop it using an event function. Now, I want to switch to the MRP shadow set and run the ode45 again with the updated initial conditions. How do I do that?
clc;
clear all;
w0 = [0; 0; 0.3];
MRP0 = [0; 0; -0.198912367379658];
tspan = [0, 100];
% Integrate the Euler equations using ode45 with an event function
options = odeset('RelTol', 1e-10, 'AbsTol', 1e-10, 'Events', @singularityEvent);
[t, y, te, ye, ie] = ode45(@KDE_MRP, tspan, [MRP0; w0], options);
% Extract the Euler parameters and angular velocities
MRP = y(:, 1:3);
w = y(:, 4:6);
% plotting MRP vs time
plot(t,MRP, '-')
xlabel('time (s)')
ylabel('MRP')
legend('MRP1', 'MRP2', 'MRP3')
function [value, isterminal, direction] = singularityEvent(t, MRPw)
% Event function to detect when the norm of MRP approaches 2
value = norm(MRPw(1:3)) - 1;
isterminal = 1; % Stop the integration when the event is detected
direction = 0; % Detect events in both directions
end
function dMRPwdt = KDE_MRP(t, MRPw)
I = [0.3; 0.2; 0.4];
L = [0; 0; 0];
MRP = MRPw(1:3);
w = MRPw(4:6);
dMRPdt = zeros(3, 1);
dwdt = zeros(3, 1);
C11 = 1 - dot(MRP, MRP) + 2 * MRP(1)^2;
C12 = 2 * (MRP(1) * MRP(2) - MRP(3));
C13 = 2 * (MRP(1) * MRP(3) + MRP(2));
C21 = 2 * (MRP(2) * MRP(1) + MRP(3));
C22 = 1 - dot(MRP, MRP) + 2 * MRP(2)^2;
C23 = 2 * (MRP(2) * MRP(3) - MRP(1));
C31 = 2 * (MRP(3) * MRP(1) - MRP(2));
C32 = 2 * (MRP(3) * MRP(2) + MRP(1));
C33 = 1 - dot(MRP, MRP) + 2 * MRP(3)^2;
C = [C11 C12 C13;
C21 C22 C23;
C31 C32 C33];
dMRPdt = 0.25 * C * w;
dwdt(1) = (-(I(3) - I(2)) * w(2) * w(3) + L(1)) / I(1);
dwdt(2) = (-(I(1) - I(3)) * w(3) * w(1) + L(2)) / I(2);
dwdt(3) = (-(I(2) - I(1)) * w(1) * w(2) + L(3)) / I(3);
% Combine the time derivatives into a single
dMRPwdt = [dMRPdt; dwdt];
end
Trending now
This is a popular solution!
Step by step
Solved in 4 steps
