balle.m

% balle - Program to compute the trajectory of a baseball
% using the Euler method.
clear;  help balle;  % Clear memory and print header
r = input('Enter initial position r=[x,y] - ');  % (meters)
v = input('Enter initial velocity v=[vx,vy] - '); % (m/sec)
tau = input('Enter timestep, tau - ');  % (sec)
Cd = 0.35;      % Drag coefficient (dimensionless)
rho = 1.2;      % Density of air (kg/m^3)
area = 4.3e-3;  % Cross-sectional area of projectile (m^2)
grav = 9.81;    % Gravitational acceleration (m/s^2)
mass = 0.145;   % Mass of projectile (kg)
air_const = -0.5*Cd*rho*area/mass;  % Air resistance constant
maxstep = 1000; % Maximum number of steps
%%%%% MAIN LOOP %%%%%
for istep=1:maxstep
  xplot(istep) = r(1);   % Record trajectory for plot
  yplot(istep) = r(2);
  accel = air_const*norm(v)*v;   % Air resistance
  accel(2) = accel(2)-grav;      % Gravity
  r = r + tau*v;                 % Euler step
  v = v + tau*accel;     
  if( r(2) < 0 )  % Break out of loop when ball hits ground
    break;        % i.e. when y-position is < 0
  end 
end
xplot(istep+1) = r(1);  yplot(istep+1) = r(2);
fprintf('Maximum range is %g meters\n',r(1))
fprintf('Time of flight is %g seconds\n',istep*tau)
% Mark the location of the ground by a straight line
xground = [0 xplot(istep+1)]; yground = [0 0];
% Graph the trajectory of the baseball
plot(xplot,yplot,'+',xground,yground,'-');
xlabel('Range (m)')
ylabel('Height (m)')
title('Projectile motion')