I need help with a MATLAB code. This code just keeps running and does not give me any plots. I even reduced the tolerance from 1e-9 to 1e-6. Can you help me fix this? Please make sure your solution runs. % Initial Conditions rev = 0:0.001:2; g1 = deg2rad(1); g2 = deg2rad(3); g3 = deg2rad(6); g4 = deg2rad(30); g0 = deg2rad(0); Z0 = 0; w0 = [0; Z0*cos(g0); -Z0*sin(g0)]; Z1 = 5; w1 = [0; Z1*cos(g1); -Z1*sin(g1)]; Z2 = 11; w2 = [0; Z2*cos(g2); -Z2*sin(g2)]; [v3, psi3, eta3] = Nut_angle(Z2, g2, w2); plot(v3, psi3) function dwedt = K_DDE(~, w_en) % Extracting the initial condtions to a variable % Extracting the initial condtions to a variable w = w_en(1:3); e = w_en(4:7); Z = w_en(8); I = 0.060214; J = 0.015707; x = (J/I) - 1; y = Z - 1; s = Z; % Kinematic Differential Equations dedt = zeros(4,1); dedt(1) = pi*(e(3)*(s-w(2)-1) + e(2)*w(3) + e(4)*w(1)); dedt(2) = pi*(e(4)*(w(2)-1-s) + e(3)*w(1) - e(1)*w(3)); dedt(3) = pi*(-e(1)*(s-w(2)-1) - e(2)*w(1) + e(4)*w(3)); dedt(4) = pi*(-e(2)*(w(2)-1-s) - e(1)*w(1) - e(3)*w(3)); % Dynamical Differential Equations dwdt = zeros(3,1); dwdt(1) = 2*pi*(6*(-x)*(e(2)*e(3)+e(1)*e(4))*(1-2*e(1)^2-2*e(2)^2) -... (-x)*w(2)*w(3) + s*w(3)); dwdt(2) = 0; dwdt(3) = 2*pi*(-12*(-x)*(e(3)*e(1)-e(2)*e(4))*(e(2)*e(3)+e(1)*e(4))+... (-x)*w(1)*w(2) + s*w(1)); % Combining the w and C into one output matrix dwedt = [dwdt; dedt; 0]; end function [EP, norm_EP] = DCMtoEP(c) % Finding EP^2 values to determine which is bigger EP_1 = (1/4)*(1 + 2*c(1,1) - trace(c)); EP_2 = (1/4)*(1 + 2*c(2,2) - trace(c)); EP_3 = (1/4)*(1 + 2*c(3,3) - trace(c)); EP_4 = (1/4)*(1 + trace(c)); EP_squared = [EP_1, EP_2, EP_3, EP_4]; % Max Value of EP_squared and its location in the matrix [maxVal, maxIndex] = max(EP_squared); C = c; % Using if statements to find EP based on the max value of EP^2 if maxIndex == 1 EP = [sqrt(maxVal); (C(1,2)+C(2,1))/(4*sqrt(maxVal)); (C(1,3)+C(3,1))/(4*sqrt(maxVal)); (C(2,3)-C(3,2))/(4*sqrt(maxVal))]; elseif maxIndex == 2 EP = [(C(1,2)+C(2,1))/(4*sqrt(maxVal)); sqrt(maxVal); (C(2,3)+C(3,2))/(4*sqrt(maxVal)); (C(3,1)-C(1,3))/(4*sqrt(maxVal))]; elseif maxIndex == 3 EP = [(C(3,1)+C(1,3))/(4*sqrt(maxVal)); (C(2,3)+C(3,2))/(4*sqrt(maxVal)); sqrt(maxVal); (C(1,2)-C(2,1))/(4*sqrt(maxVal))]; elseif maxIndex == 4 EP = [(C(2,3)-C(3,2))/(4*sqrt(maxVal)); (C(3,1)-C(1,3))/(4*sqrt(maxVal)); (C(1,2)-C(2,1))/(4*sqrt(maxVal)); sqrt(maxVal)]; end % Finding the norm of EP norm_EP = norm(EP); end function [v, psi, eta] = Nut_angle(n, g, w) % Initial Conditions rev = 0:0.001:2; C = [cos(g) -sin(g) 0; sin(g) cos(g) 0; 0 0 1]; [e, ~] = DCMtoEP(C); % Using ode45 to integrate the KDE and DDE options = odeset('RelTol',1e-6,'AbsTol',1e-6); result = ode45(@K_DDE, rev, [w; e; n], options); v = result.x; % Extracting information from the ode45 solver e_ode = result.y(4:7, :); % Finding C21, C22, AND C23 C21 = 2*(e_ode(1,:).*e_ode(2,:) + e_ode(3,:).*e_ode(4,:)); C22 = 1 - 2*e_ode(3,:).^2 - 2*e_ode(1,:).^2; C23 = 2*(e_ode(2,:).*e_ode(3,:) - e_ode(1,:).*e_ode(4,:)); C32 = 2*(e_ode(2,:).*e_ode(3,:) + e_ode(1,:).*e_ode(4,:)); C12 = 2*(e_ode(1,:).*e_ode(2,:) - e_ode(3,:).*e_ode(4,:)); % Initializing arrays psi = zeros(1, length(v)); eta = zeros(1, length(v)); % Finding the nutation angles for i = 1:length(v) psi(i) = acosd(C22(i)); eta(i) = atan2d(C12(i), C32(i)); end end
I need help with a MATLAB code. This code just keeps running and does not give me any plots. I even reduced the tolerance from 1e-9 to 1e-6. Can you help me fix this? Please make sure your solution runs.
% Initial Conditions
rev = 0:0.001:2;
g1 = deg2rad(1);
g2 = deg2rad(3);
g3 = deg2rad(6);
g4 = deg2rad(30);
g0 = deg2rad(0);
Z0 = 0;
w0 = [0; Z0*cos(g0); -Z0*sin(g0)];
Z1 = 5;
w1 = [0; Z1*cos(g1); -Z1*sin(g1)];
Z2 = 11;
w2 = [0; Z2*cos(g2); -Z2*sin(g2)];
[v3, psi3, eta3] = Nut_angle(Z2, g2, w2);
plot(v3, psi3)
function dwedt = K_DDE(~, w_en)
% Extracting the initial condtions to a variable
% Extracting the initial condtions to a variable
w = w_en(1:3);
e = w_en(4:7);
Z = w_en(8);
I = 0.060214;
J = 0.015707;
x = (J/I) - 1;
y = Z - 1;
s = Z;
% Kinematic Differential Equations
dedt = zeros(4,1);
dedt(1) = pi*(e(3)*(s-w(2)-1) + e(2)*w(3) + e(4)*w(1));
dedt(2) = pi*(e(4)*(w(2)-1-s) + e(3)*w(1) - e(1)*w(3));
dedt(3) = pi*(-e(1)*(s-w(2)-1) - e(2)*w(1) + e(4)*w(3));
dedt(4) = pi*(-e(2)*(w(2)-1-s) - e(1)*w(1) - e(3)*w(3));
% Dynamical Differential Equations
dwdt = zeros(3,1);
dwdt(1) = 2*pi*(6*(-x)*(e(2)*e(3)+e(1)*e(4))*(1-2*e(1)^2-2*e(2)^2) -...
(-x)*w(2)*w(3) + s*w(3));
dwdt(2) = 0;
dwdt(3) = 2*pi*(-12*(-x)*(e(3)*e(1)-e(2)*e(4))*(e(2)*e(3)+e(1)*e(4))+...
(-x)*w(1)*w(2) + s*w(1));
% Combining the w and C into one output matrix
dwedt = [dwdt; dedt; 0];
end
function [EP, norm_EP] = DCMtoEP(c)
% Finding EP^2 values to determine which is bigger
EP_1 = (1/4)*(1 + 2*c(1,1) - trace(c));
EP_2 = (1/4)*(1 + 2*c(2,2) - trace(c));
EP_3 = (1/4)*(1 + 2*c(3,3) - trace(c));
EP_4 = (1/4)*(1 + trace(c));
EP_squared = [EP_1, EP_2, EP_3, EP_4];
% Max Value of EP_squared and its location in the matrix
[maxVal, maxIndex] = max(EP_squared);
C = c;
% Using if statements to find EP based on the max value of EP^2
if maxIndex == 1
EP = [sqrt(maxVal);
(C(1,2)+C(2,1))/(4*sqrt(maxVal));
(C(1,3)+C(3,1))/(4*sqrt(maxVal));
(C(2,3)-C(3,2))/(4*sqrt(maxVal))];
elseif maxIndex == 2
EP = [(C(1,2)+C(2,1))/(4*sqrt(maxVal));
sqrt(maxVal);
(C(2,3)+C(3,2))/(4*sqrt(maxVal));
(C(3,1)-C(1,3))/(4*sqrt(maxVal))];
elseif maxIndex == 3
EP = [(C(3,1)+C(1,3))/(4*sqrt(maxVal));
(C(2,3)+C(3,2))/(4*sqrt(maxVal));
sqrt(maxVal);
(C(1,2)-C(2,1))/(4*sqrt(maxVal))];
elseif maxIndex == 4
EP = [(C(2,3)-C(3,2))/(4*sqrt(maxVal));
(C(3,1)-C(1,3))/(4*sqrt(maxVal));
(C(1,2)-C(2,1))/(4*sqrt(maxVal));
sqrt(maxVal)];
end
% Finding the norm of EP
norm_EP = norm(EP);
end
function [v, psi, eta] = Nut_angle(n, g, w)
% Initial Conditions
rev = 0:0.001:2;
C = [cos(g) -sin(g) 0;
sin(g) cos(g) 0;
0 0 1];
[e, ~] = DCMtoEP(C);
% Using ode45 to integrate the KDE and DDE
options = odeset('RelTol',1e-6,'AbsTol',1e-6);
result = ode45(@K_DDE, rev, [w; e; n], options);
v = result.x;
% Extracting information from the ode45 solver
e_ode = result.y(4:7, :);
% Finding C21, C22, AND C23
C21 = 2*(e_ode(1,:).*e_ode(2,:) + e_ode(3,:).*e_ode(4,:));
C22 = 1 - 2*e_ode(3,:).^2 - 2*e_ode(1,:).^2;
C23 = 2*(e_ode(2,:).*e_ode(3,:) - e_ode(1,:).*e_ode(4,:));
C32 = 2*(e_ode(2,:).*e_ode(3,:) + e_ode(1,:).*e_ode(4,:));
C12 = 2*(e_ode(1,:).*e_ode(2,:) - e_ode(3,:).*e_ode(4,:));
% Initializing arrays
psi = zeros(1, length(v));
eta = zeros(1, length(v));
% Finding the nutation angles
for i = 1:length(v)
psi(i) = acosd(C22(i));
eta(i) = atan2d(C12(i), C32(i));
end
end
Step by step
Solved in 2 steps
