-
Notifications
You must be signed in to change notification settings - Fork 3
/
projectile.m
70 lines (68 loc) · 1.94 KB
/
projectile.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
% Projectile Motion:
% Program to compute the 2D projectile of a baseball using
% the Euler, Euler-Cromer, and Midpoint methods.
%
clear, clc
% Initial conditions
x0=[0 0];
theta=45*pi/180;
v0=100;
v0=[v0*cos(theta) v0*sin(theta)];
x=x0;
v=v0;
state=[x v]; % Set the initial state vector
time=0;
method=menu('Numerical method','Midpoint','Runge Kutta','Adaptive Runge Kutta');
% Other parameters
grav=[0,-9.8]; % acceleration of gravity
cd = 0.35; % drag coefficient
r = .037; % radius of baseball
area = pi*r*r; % frontal area of baseball
mass = 0.145; % mass of baseball
rho = 1.2; % air density
tau=0.1; % time step size
dat=[time x v]; % solution matrix
dragconst=0.5*rho*area/mass;
adaptErr=5e-4;
% Loop to compute trajectory
counter = 0;
while(x(2)>=0 & counter<10000) % repeat until the projectile hits the ground
counter = counter + 1;
if(method==1) % Midpoint
accel = grav - dragconst*v*norm(v);
vnew = v + tau*accel;
x = x + tau*0.5*(v+vnew);
v = vnew;
time = time + tau;
elseif(method==2) % Runge Kutta
% Call rk4
% x = ? % extract position vector from the state vector
% v = ? % extract velocity vector from the state vector
time = time + tau;
else % Adaptive Runge Kutta
% Call rka
% x = ? % extract position vector from the state vector
% v = ? % extract velocity vector from the state vector
end
dat=[dat;time x v];
end
%
if(method==1)
col='k-';
elseif(method==2)
col='bp';
else
col='rd';
end
subplot(221)
plot(dat(:,1),dat(:,4),col)
xlabel('time (sec)'), ylabel('v_x (m/s)')
hold on
subplot(222)
plot(dat(:,1),dat(:,5),col)
xlabel('time (sec)'), ylabel('v_y (m/s)')
hold on
subplot(212)
plot(dat(:,2),dat(:,3),col)
xlabel('x (m)'), ylabel('y (m)')
hold on