I need help with my MATLAB code. I want to run ode45 for each set of EPs. So, I want to run it 30 times. Is there any way I can just put the ode45 line in the for loop. I dont want to run the ode45 line 30 different times with different sets of EPs. clc; clear all; % Initial conditions mu = 398600; % km^3/s^2 R = 7000; % km I = [300; 500; 700] * 10^-6; % km^3/s^2 % Initial PRP and velocity vectors lambda = [1/sqrt(3); 1/sqrt(3); 1/sqrt(3)]; theta = 1:30; % deg w = [0; 0; sqrt(mu/R^3)]; t = 0:43200; % sec % Finding EP for i = 1:length(theta) % Finding EP from PRP EP1(i) = lambda(1)*sind(theta(i)/2); EP2(i) = lambda(2)*sind(theta(i)/2); EP3(i) = lambda(3)*sind(theta(i)/2); EP4(i) = cosd(theta(i)/2); EP(:,i) = [EP1(i); EP2(i); EP3(i); EP4(i)]; end % Using ode45 to integrate KDE and EOM options = odeset('RelTol',1e-10,'AbsTol',1e-10); [t, y] = ode45(@dwdt_KDE_EP, t, [EP; w; I], options); function dqwdt = dwdt_KDE_EP(t,EPwI) EP = EPwI(1:4); w = EPwI(5:7); I = EPwI(8:10); dqdt = zeros(4,1); dwdt = zeros(3,1); dqdt(1) = 0.5*(EP(4)*w(1) - EP(3)*w(2) + EP(2)*w(3)); dqdt(2) = 0.5*(EP(3)*w(1) + EP(4)*w(2) - EP(1)*w(3)); dqdt(3) = 0.5*(-EP(2)*w(1) + EP(1)*w(2) + EP(4)*w(3)); dqdt(4) = -0.5*(EP(1)*w(1) + EP(2)*w(2) + EP(3)*w(3)); 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 - 2*EP(2)^2 - 2*EP(3)^2); R2 = R_mag * (2*EP(1)*EP(2) - EP(3)*EP(4)); R3 = R_mag * (2*EP(1)*EP(3) + EP(2)*EP(4)); 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 dqwdt = [dqdt; dwdt; 0;0;0]; end
I need help with my MATLAB code. I want to run ode45 for each set of EPs. So, I want to run it 30 times. Is there any way I can just put the ode45 line in the for loop. I dont want to run the ode45 line 30 different times with different sets of EPs.
clc;
clear all;
% Initial conditions
mu = 398600; % km^3/s^2
R = 7000; % km
I = [300; 500; 700] * 10^-6; % km^3/s^2
% Initial PRP and velocity
lambda = [1/sqrt(3); 1/sqrt(3); 1/sqrt(3)];
theta = 1:30; % deg
w = [0; 0; sqrt(mu/R^3)];
t = 0:43200; % sec
% Finding EP
for i = 1:length(theta)
% Finding EP from PRP
EP1(i) = lambda(1)*sind(theta(i)/2);
EP2(i) = lambda(2)*sind(theta(i)/2);
EP3(i) = lambda(3)*sind(theta(i)/2);
EP4(i) = cosd(theta(i)/2);
EP(:,i) = [EP1(i); EP2(i); EP3(i); EP4(i)];
end
% Using ode45 to integrate KDE and EOM
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t, y] = ode45(@dwdt_KDE_EP, t, [EP; w; I], options);
function dqwdt = dwdt_KDE_EP(t,EPwI)
EP = EPwI(1:4);
w = EPwI(5:7);
I = EPwI(8:10);
dqdt = zeros(4,1);
dwdt = zeros(3,1);
dqdt(1) = 0.5*(EP(4)*w(1) - EP(3)*w(2) + EP(2)*w(3));
dqdt(2) = 0.5*(EP(3)*w(1) + EP(4)*w(2) - EP(1)*w(3));
dqdt(3) = 0.5*(-EP(2)*w(1) + EP(1)*w(2) + EP(4)*w(3));
dqdt(4) = -0.5*(EP(1)*w(1) + EP(2)*w(2) + EP(3)*w(3));
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 - 2*EP(2)^2 - 2*EP(3)^2);
R2 = R_mag * (2*EP(1)*EP(2) - EP(3)*EP(4));
R3 = R_mag * (2*EP(1)*EP(3) + EP(2)*EP(4));
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
dqwdt = [dqdt; dwdt; 0;0;0];
end
Trending now
This is a popular solution!
Step by step
Solved in 4 steps with 1 images
