Stopping plot when y is less than zero
조회 수: 7 (최근 30일)
이전 댓글 표시
I have this code for projectile motion, but I need to stop the plot when y<0. How could I do this?
global m cd A rho g theta v0 x0 y0
m= [0.145 0.42 0.43]; % mass of projectile [kg]
cd= [0.3 0.5 0.25]; % drag coefficients
A= [0.00426 0.02318 0.038]; % cross sectional area of projectile [m^2]
V= [0.0002194 0.00462115 0.00573547]; % volume of projectile
rho= [m(1)/V(1) m(2)/V(2) m(3)/V(3)]; % density
g= 9.8; % gravity [m/s^2]
theta= 30*pi/180; % initial launch angle [rad]
x0= [0 0 0]; % initial horizontal position
y0= [1.8288 1.8288 0]; % initial vertical position
v0= [18 18 22]; % initial velocity (m/s)
%set projectile, 1=baseball, 2=football, 3=soccer ball
proj=3;
m= m(proj);
cd= cd(proj);
A= A(proj);
rho= rho(proj);
x0= x0(proj);
y0= y0(proj);
v0= v0(proj);
% call function
[t y]= ode45(@my_fun,[0 2],[x0 y0 v0 theta]);
% Verification
max_y= max(y(:,2)); % maximum vertical height of projectile [m]
% plotting results
figure(1)
plot(t,y(:,1))
grid on
xlabel('time [s]')
ylabel('horizontal position [m]')
title('plot of horizontal location vs time')
figure(2)
plot(t,y(:,2))
grid on
xlabel('time [s]')
ylabel('vertical position [m]')
title('plot of vertical location vs time')
figure(3)
plot(y(:,1),y(:,2))
grid on
xlabel('horizontal position [m]')
ylabel('vertical position [m]')
title('plot of vertical location vs horizontal position')
end
%% Function
function dy=my_fun(t,y) % y is a matrix containing [x y v theta]
global m cd A rho g theta x0 y0 v0
dy=zeros(4,1);
dy(1)= y(2); %dx/dt
dy(2)= y(3); %dy/dt
dy(3)= -(cd*A*rho)/(2*m)*sqrt((dy(1))^2+(dy(2))^2)*dy(1); % d2x/dt2
dy(4)= -g-(cd*A*rho)/(2*m)*sqrt((dy(1))^2+(dy(2))^2)*dy(2); % d2y/dt2
end
댓글 수: 0
채택된 답변
Ameer Hamza
2020년 4월 29일
You can control it using events in ODE: https://www.mathworks.com/help/releases/R2020a/matlab/ref/odeset.html#d120e857760
Define the following function in your code
function [event,isterminal,direction] = odeEvent(t,y)
event = y(2); % observe y(2)=0
isterminal = 1; % stop on observing y(2)=0
direction = -1; % observe that y(2) is decreasing before the event y(2)=0
end
and then call ode45 like this
odeOpt = odeset('Events', @odeEvent);
[t y]= ode45(@my_fun,[0 2],[x0 y0 v0 theta], odeOpt);
댓글 수: 0
추가 답변 (1개)
Image Analyst
2020년 4월 29일
DON'T use cd as the name of your variable since it's the name of a built-in function.
Try masking by using a logical index that's 1 (true) where y is >0 and 0 (false) where it's below 0:
global m cd A rho g theta v0 x0 y0
m = [0.145 0.42 0.43]; % mass of projectile [kg]
cd = [0.3 0.5 0.25]; % drag coefficients
A = [0.00426 0.02318 0.038]; % cross sectional area of projectile [m^2]
V = [0.0002194 0.00462115 0.00573547]; % volume of projectile
rho = [m(1)/V(1) m(2)/V(2) m(3)/V(3)]; % density
g = 9.8; % gravity [m/s^2]
theta = 30*pi/180; % initial launch angle [rad]
x0 = [0 0 0]; % initial horizontal position
y0 = [1.8288 1.8288 0]; % initial vertical position
v0 = [18 18 22]; % initial velocity (m/s)
%set projectile, 1=baseball, 2=football, 3=soccer ball
proj=3;
m= m(proj);
cd= cd(proj);
A= A(proj);
rho= rho(proj);
x0= x0(proj);
y0= y0(proj);
v0= v0(proj);
% call function
[t y]= ode45(@my_fun,[0 2],[x0 y0 v0 theta]);
% Verification
max_y= max(y(:,2)); % maximum vertical height of projectile [m]
% plotting results
hFig = figure;
subplot(2, 2, 1);
hFig.WindowState = 'maximized';
positiveIndexes = y(:, 1) > 0;
plot(t(positiveIndexes), y(positiveIndexes,1), 'b.-', 'LineWidth', 2, 'MarkerSize', 18)
grid on
xlabel('time [s]')
ylabel('horizontal position [m]')
title('plot of horizontal location vs time')
subplot(2, 2, 2);
positiveIndexes = y(:, 2) > 0;
plot(t(positiveIndexes), y(positiveIndexes,2), 'b.-', 'LineWidth', 2, 'MarkerSize', 18)
grid on
xlabel('time [s]')
ylabel('vertical position [m]')
title('plot of vertical location vs time')
subplot(2, 2, 3);
positiveIndexes = (y(:, 1) > 0) & (y(:, 2) > 0);
plot(y(positiveIndexes,1), y(positiveIndexes,2), 'b.-', 'LineWidth', 2, 'MarkerSize', 18)
grid on
xlabel('horizontal position [m]')
ylabel('vertical position [m]')
title('plot of vertical location vs horizontal position')
fprintf('Done running %s.m.\n', mfilename);
%% Function
function dy=my_fun(t,y) % y is a matrix containing [x y v theta]
global m cd A rho g theta x0 y0 v0
dy=zeros(4,1);
dy(1)= y(2); %dx/dt
dy(2)= y(3); %dy/dt
dy(3)= -(cd*A*rho)/(2*m)*sqrt((dy(1))^2+(dy(2))^2)*dy(1); % d2x/dt2
dy(4)= -g-(cd*A*rho)/(2*m)*sqrt((dy(1))^2+(dy(2))^2)*dy(2); % d2y/dt2
end
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/288685/image.png)
Also, for fun, I'm attaching my projectile demo that computes just about everything anyone could possibly want to know about a projectile's trajectory.
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!