I need help with my MATLAB code. This code produces an angular velocity plot. I want the x-axis of the plot to go from 0 to 12 hours instead of 0 to 43200 seconds. But I still want the input to be t = 0:43200. I just want to change what it says when the code produces the plot. Also, on the y-axis, I need it to go from -0.05 to 0.1. Can you help me with that? % Initial conditions mu = 398600; % km^3/s^2 R = 6800; % km I = [400; 600; 800] * 10^-6; % kg*km^2 % Initial PRP and velocity vectors lambda = [1/sqrt(3); 1/sqrt(3); 1/sqrt(3)]; theta = 3; % deg w = [0; 0; sqrt(mu/R^3)]; t = 0:43200; % sec % Finding MRP MRP = PRP2MRP(lambda, theta) % Problem 1 (e.2) % Integrate the Euler equations using ode45 options = odeset('RelTol',1e-10,'AbsTol',1e-10); [t, y] = ode45(@KDE_MRP, t, [MRP; w; I], options); % Extract the Euler parameters and angular velocities MRP2 = y(:, 1:3); w2 = y(:, 4:6); plot(t,w2, '-') xlabel('time (s)') ylabel('angular velocity (rad/s)') legend('w1', 'w2', 'w3') function MRP = PRP2MRP(lambda, theta) % Finding EP from PRP EP1 = lambda(1)*sind(theta/2); EP2 = lambda(2)*sind(theta/2); EP3 = lambda(3)*sind(theta/2); EP4 = cosd(theta/2); EP = [EP1; EP2; EP3; EP4]; % Finding MRP from EP sigma1 = EP(1)/(1+EP(4)); sigma2 = EP(2)/(1+EP(4)); sigma3 = EP(3)/(1+EP(4)); MRP = [sigma1; sigma2; sigma3]; end function dMRPwdt = KDE_MRP(~,MRPwI) MRP = MRPwI(1:3); w = MRPwI(4:6); I = MRPwI(7:9); 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; K1 = (I(2) - I(3)) / I(1); K2 = (I(3) - I(1)) / I(2); K3 = (I(1) - I(2)) / I(3); R_mag = 6800; R1 = (R_mag/(1+dot(MRP, MRP))^2) * (4*(MRP(1)^2 - MRP(2)^2 - MRP(3)^2) + (1 - dot(MRP, MRP))); R2 = (R_mag/(1+dot(MRP, MRP))^2) * (8*MRP(2)*MRP(1) - 4*MRP(3)*(1 - dot(MRP, MRP))); R3 = (R_mag/(1+dot(MRP, MRP))^2) * (8*MRP(3)*MRP(1) - 4*MRP(2)*(1 - dot(MRP, MRP))); dwdt(1) = K1*w(2)*w(3) - ((3*w(3)^2*K1*R3*R2)/R_mag^2); dwdt(2) = K2*w(1)*w(3) - ((3*w(3)^2*K2*R1*R3)/R_mag^2); dwdt(3) = K3*w(1)*w(2) - ((3*w(3)^2*K3*R2*R1)/R_mag^2); % Combine the time derivatives into a single vector dMRPwdt = [dMRPdt; dwdt; 0;0;0]; end
I need help with my MATLAB code. This code produces an angular velocity plot. I want the x-axis of the plot to go from 0 to 12 hours instead of 0 to 43200 seconds. But I still want the input to be t = 0:43200. I just want to change what it says when the code produces the plot. Also, on the y-axis, I need it to go from -0.05 to 0.1. Can you help me with that?
% Initial conditions
mu = 398600; % km^3/s^2
R = 6800; % km
I = [400; 600; 800] * 10^-6; % kg*km^2
% Initial PRP and velocity
lambda = [1/sqrt(3); 1/sqrt(3); 1/sqrt(3)];
theta = 3; % deg
w = [0; 0; sqrt(mu/R^3)];
t = 0:43200; % sec
% Finding MRP
MRP = PRP2MRP(lambda, theta)
% Problem 1 (e.2)
% Integrate the Euler equations using ode45
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t, y] = ode45(@KDE_MRP, t, [MRP; w; I], options);
% Extract the Euler parameters and angular velocities
MRP2 = y(:, 1:3);
w2 = y(:, 4:6);
plot(t,w2, '-')
xlabel('time (s)')
ylabel('angular velocity (rad/s)')
legend('w1', 'w2', 'w3')
function MRP = PRP2MRP(lambda, theta)
% Finding EP from PRP
EP1 = lambda(1)*sind(theta/2);
EP2 = lambda(2)*sind(theta/2);
EP3 = lambda(3)*sind(theta/2);
EP4 = cosd(theta/2);
EP = [EP1; EP2; EP3; EP4];
% Finding MRP from EP
sigma1 = EP(1)/(1+EP(4));
sigma2 = EP(2)/(1+EP(4));
sigma3 = EP(3)/(1+EP(4));
MRP = [sigma1; sigma2; sigma3];
end
function dMRPwdt = KDE_MRP(~,MRPwI)
MRP = MRPwI(1:3);
w = MRPwI(4:6);
I = MRPwI(7:9);
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;
K1 = (I(2) - I(3)) / I(1);
K2 = (I(3) - I(1)) / I(2);
K3 = (I(1) - I(2)) / I(3);
R_mag = 6800;
R1 = (R_mag/(1+dot(MRP, MRP))^2) * (4*(MRP(1)^2 - MRP(2)^2 - MRP(3)^2) + (1 - dot(MRP, MRP)));
R2 = (R_mag/(1+dot(MRP, MRP))^2) * (8*MRP(2)*MRP(1) - 4*MRP(3)*(1 - dot(MRP, MRP)));
R3 = (R_mag/(1+dot(MRP, MRP))^2) * (8*MRP(3)*MRP(1) - 4*MRP(2)*(1 - dot(MRP, MRP)));
dwdt(1) = K1*w(2)*w(3) - ((3*w(3)^2*K1*R3*R2)/R_mag^2);
dwdt(2) = K2*w(1)*w(3) - ((3*w(3)^2*K2*R1*R3)/R_mag^2);
dwdt(3) = K3*w(1)*w(2) - ((3*w(3)^2*K3*R2*R1)/R_mag^2);
% Combine the time derivatives into a single vector
dMRPwdt = [dMRPdt; dwdt; 0;0;0];
end
Step by step
Solved in 4 steps with 5 images
