am trying to plot the Euler parameters constraint equation vs revolutions and plot nutation, precession, and spin angles vs time. I am unable to run this code with the [q1; 0; 0; q4;0;1;0] initial conditions. I cannot run the code with initial condition [q1; 0; 0; q4;0;1;0] or [q1; 0; 0; q4;0;-1;0]. But it works fine for any higher whole number such as [q1; 0; 0; q4;0;2;0]. Why is that? How can I fix this? clc; clear all; q1 = sin((6/2)*pi/180); q4 = cos((6/2)*pi/180); t = 0:0.01:(2*pi*4); options=odeset('RelTol',1e-12,'AbsTol',1e-12); [t,y] = ode45(@euler_eqns, t, [q1; 0; 0; q4;0;1;0], options); q = y(:,1:4) w = y(:,5:7) K = q(:,1).^2 + q(:,2).^2 + q(:,3).^2 + q(:,4).^2; K0 = 1; c = K - K0; v = t / (2*pi); plot(v,c); xlabel('revs '); ylabel('K-K0'); nutation = acos(1 - 2*q(:,1).^2 - 2*q(:,3).^2); precession = asin((2*(q(:,1).*q(:,2) - q(:,3).*q(:,4)))./(sin(nutation))); spin = asin((2*(q(:,1).*q(:,2) + q(:,3).*q(:,4)))./(sin(nutation))); psi = rad2deg(precession); theta = rad2deg(nutation); phi = rad2deg(spin); % Plot the angles over time plot(v, psi, v, theta, v, phi) legend('Precession', 'Nutation', 'Spin') xlabel('revs') ylabel('Angles (deg)') function dqwdt = euler_eqns(t,qwd, q) q = qwd(1:4); w = qwd(5:7); I = 500; J = 125; K = (I-J)/I; T = 35; ohm = 1; qstar = 1/ohm; dqdt = zeros(4,1); dqdt(1) = pi*(q(4)*w(1) + q(3)*(qstar - w(2) - ohm) + q(2)*w(3)); dqdt(2) = pi*(q(3)*w(1) + q(4)*(w(2) - ohm - qstar) - q(1)*w(3)); dqdt(3) = pi*(-q(2)*w(1) - q(1)*(qstar - w(2) - ohm) + q(4)*w(3)); dqdt(4) = pi*(-q(1)*w(1) - q(2)*(w(2) - ohm - qstar) - q(3)*w(3)); dwdt = zeros(3,1); dwdt(1) = 12*pi*K*(q(2)*q(3) + q(1)*q(4))*(1 - 2*q(1)^2 - 2*q(2)^2) ... - 2*pi*K*w(2)*w(3) + 2*pi*qstar*w(3); dwdt(2) = 0; dwdt(3) = -24*pi*K*(q(3)*q(1) - q(2)*q(4))*(q(2)*q(3) + q(1)*q(4)) ... + 2*pi*K*w(1)*w(2) + 2*pi*qstar*w(1); % Combine the time derivatives into a single vector dqwdt = [dqdt; dwdt]; end
I am trying to plot the Euler parameters constraint equation vs revolutions and plot nutation, precession, and spin angles vs time. I am unable to run this code with the [q1; 0; 0; q4;0;1;0] initial conditions. I cannot run the code with initial condition [q1; 0; 0; q4;0;1;0] or [q1; 0; 0; q4;0;-1;0]. But it works fine for any higher whole number such as [q1; 0; 0; q4;0;2;0]. Why is that? How can I fix this?
clc;
clear all;
q1 = sin((6/2)*pi/180);
q4 = cos((6/2)*pi/180);
t = 0:0.01:(2*pi*4);
options=odeset('RelTol',1e-12,'AbsTol',1e-12);
[t,y] = ode45(@euler_eqns, t, [q1; 0; 0; q4;0;1;0], options);
q = y(:,1:4)
w = y(:,5:7)
K = q(:,1).^2 + q(:,2).^2 + q(:,3).^2 + q(:,4).^2;
K0 = 1;
c = K - K0;
v = t / (2*pi);
plot(v,c);
xlabel('revs ');
ylabel('K-K0');
nutation = acos(1 - 2*q(:,1).^2 - 2*q(:,3).^2);
precession = asin((2*(q(:,1).*q(:,2) - q(:,3).*q(:,4)))./(sin(nutation)));
spin = asin((2*(q(:,1).*q(:,2) + q(:,3).*q(:,4)))./(sin(nutation)));
psi = rad2deg(precession);
theta = rad2deg(nutation);
phi = rad2deg(spin);
% Plot the angles over time
plot(v, psi, v, theta, v, phi)
legend('Precession', 'Nutation', 'Spin')
xlabel('revs')
ylabel('Angles (deg)')
function dqwdt = euler_eqns(t,qwd, q)
q = qwd(1:4);
w = qwd(5:7);
I = 500;
J = 125;
K = (I-J)/I;
T = 35;
ohm = 1;
qstar = 1/ohm;
dqdt = zeros(4,1);
dqdt(1) = pi*(q(4)*w(1) + q(3)*(qstar - w(2) - ohm) + q(2)*w(3));
dqdt(2) = pi*(q(3)*w(1) + q(4)*(w(2) - ohm - qstar) - q(1)*w(3));
dqdt(3) = pi*(-q(2)*w(1) - q(1)*(qstar - w(2) - ohm) + q(4)*w(3));
dqdt(4) = pi*(-q(1)*w(1) - q(2)*(w(2) - ohm - qstar) - q(3)*w(3));
dwdt = zeros(3,1);
dwdt(1) = 12*pi*K*(q(2)*q(3) + q(1)*q(4))*(1 - 2*q(1)^2 - 2*q(2)^2) ...
- 2*pi*K*w(2)*w(3) + 2*pi*qstar*w(3);
dwdt(2) = 0;
dwdt(3) = -24*pi*K*(q(3)*q(1) - q(2)*q(4))*(q(2)*q(3) + q(1)*q(4)) ...
+ 2*pi*K*w(1)*w(2) + 2*pi*qstar*w(1);
% Combine the time derivatives into a single
dqwdt = [dqdt; dwdt];
end
Trending now
This is a popular solution!
Step by step
Solved in 3 steps
