I need help with my MATLAB code. In the below code, I am trying to run the ode45 for different sets of EP's for initial theta = 0:30. Then, I am converting EP's to lambda's and theta's. I want to plot the max theta for each different set of EP's vs the corresponding initial theta. So, I want to know the max theta for initial theta = 0:30 and plot against each other. Can you help me with that % 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 % Preallocate space for EP matrix EP = zeros(4, length(theta)); for i = 1:length(theta) % Finding EP from PRP EP(1,i) = lambda(1)*sind(theta(i)/2); EP(2,i) = lambda(2)*sind(theta(i)/2); EP(3,i) = lambda(3)*sind(theta(i)/2); EP(4,i) = cosd(theta(i)/2); end % Using ode45 to integrate KDE and EOM % Preallocate space for results results = cell(1, length(theta)); options = odeset('RelTol',1e-10,'AbsTol',1e-10); for i = 1:length(theta) % Set initial conditions for each iteration initial_conditions = [EP(:, i); w; I]; % Solve ODE for each set of EPs [t, y] = ode45(@dwdt_KDE_EP, t, initial_conditions, options); % Store the results results{i}.t = t; results{i}.y = y; end % % Access and plot the results for the first set of EPs % result_idx = 1; % Change this index to access results for a different set of EPs % % % Extract time and solution vectors % t_result = results{result_idx}.t; % y_result = results{result_idx}.y; % EP1 = y_result(:,1); % EP2 = y_result(:,2); % EP3 = y_result(:,3); % EP4 = y_result(:,4); for i = 1:length(ttheta) result_idx(i) = i; % Change this index to access results for a different set of EPs % Extract time and solution vectors t_result(i) = results{result_idx(i)}.t; y_result(i) = results{result_idx(i)}.y; EP1(i) = y_result(:,1); EP2(i) = y_result(:,2); EP3(i) = y_result(:,3); EP4(i) = y_result(:,4); lambda1(i) = EP1(i) / (sqrt(EP1(i)^2 + EP2(i)^2 + EP3(i)^2)); lambda2(i) = EP2(i) / (sqrt(EP1(i)^2 + EP2(i)^2 + EP3(i)^2)); lambda3(i) = EP3(i) / (sqrt(EP1(i)^2 + EP2(i)^2 + EP3(i)^2)); theta(i) = 2*acosd(EP4(i)); max_theta(i) = max(theta(i)); end 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. In the below code, I am trying to run the ode45 for different sets of EP's for initial theta = 0:30. Then, I am converting EP's to lambda's and theta's. I want to plot the max theta for each different set of EP's vs the corresponding initial theta. So, I want to know the max theta for initial theta = 0:30 and plot against each other. Can you help me with that
% 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
% Preallocate space for EP matrix
EP = zeros(4, length(theta));
for i = 1:length(theta)
% Finding EP from PRP
EP(1,i) = lambda(1)*sind(theta(i)/2);
EP(2,i) = lambda(2)*sind(theta(i)/2);
EP(3,i) = lambda(3)*sind(theta(i)/2);
EP(4,i) = cosd(theta(i)/2);
end
% Using ode45 to integrate KDE and EOM
% Preallocate space for results
results = cell(1, length(theta));
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
for i = 1:length(theta)
% Set initial conditions for each iteration
initial_conditions = [EP(:, i); w; I];
% Solve ODE for each set of EPs
[t, y] = ode45(@dwdt_KDE_EP, t, initial_conditions, options);
% Store the results
results{i}.t = t;
results{i}.y = y;
end
% % Access and plot the results for the first set of EPs
% result_idx = 1; % Change this index to access results for a different set of EPs
%
% % Extract time and solution vectors
% t_result = results{result_idx}.t;
% y_result = results{result_idx}.y;
% EP1 = y_result(:,1);
% EP2 = y_result(:,2);
% EP3 = y_result(:,3);
% EP4 = y_result(:,4);
for i = 1:length(ttheta)
result_idx(i) = i; % Change this index to access results for a different set of EPs
% Extract time and solution vectors
t_result(i) = results{result_idx(i)}.t;
y_result(i) = results{result_idx(i)}.y;
EP1(i) = y_result(:,1);
EP2(i) = y_result(:,2);
EP3(i) = y_result(:,3);
EP4(i) = y_result(:,4);
lambda1(i) = EP1(i) / (sqrt(EP1(i)^2 + EP2(i)^2 + EP3(i)^2));
lambda2(i) = EP2(i) / (sqrt(EP1(i)^2 + EP2(i)^2 + EP3(i)^2));
lambda3(i) = EP3(i) / (sqrt(EP1(i)^2 + EP2(i)^2 + EP3(i)^2));
theta(i) = 2*acosd(EP4(i));
max_theta(i) = max(theta(i));
end
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
Step by step
Solved in 3 steps with 4 images
