% This script calculates the trajectory of a ballistic projectile.
% Assumptions:
%   Projectile starts at location (0,0) in the x-y plane
%   Initial velocity is 125 m/s at a 60 degree angle to the horizontal
%   No air resistance
%   Flat surface
 
clc

% Ballistic projectile equations:
%   x(t) = 1/2*ax*t^2 + vx*t + x0
%   y(t) = 1/2*ay*t^2 + vy*t + y0
 
% Initialize variables with physical values
g = 9.81;    % m/s^2 (gravitational constant)
theta = 60;  % degrees (projectile launch angle)
v0 = 125;    % m/s (projectile launch speed)
 
ax = 0;      % m/s^2 (assume no acceleration in x direction)
ay = -g;     % m/s^2 (assume only gravity affects projectile)
 
vx = v0 * cosd(theta);  % m/s (component of velocity in x direction)
vy = v0 * sind(theta);  % m/s (component of velocity in y direction)
 
x0 = 0;      % m (initial projectile x coordinate)
y0 = 0;      % m (initial projectile y coordinate)
 
% ***** Determine maximum projectile range *****
 
% Find t at which y(t) is zero.
%   0 = 1/2*ay*t^2 + vy*t + y0

Troots = roots([ay/2  vy  y0]);
t = max(Troots(2));
 
Rmax = 1/2*ax*t^2 + vx*t + x0; % Maximum range
fprintf('Maximum range = %0.2f m\n\n', Rmax)
 
% ***** Determine maximum projectile height *****
 
% Find t at which vertical velocity is zero.
%   0 = ay*t + vy
 
t = -vy / ay;
 
Hmax = 1/2*ay*t^2 + vy*t + y0; % Maximum height
fprintf('Maximum height = %0.2f m\n', Hmax)
 
 
% ***** Plot the trajectory H versus R *****
 
R = [0:50:1400];
H = (ay/2)/(vx^2)*R.^2 + (vy/vx)*R + y0;

clf
plot(R,H)
 
grid on
xlabel('Range R (m)')
ylabel('Height H (m)')
title([{'Projectile Trajectory','v_{0} = 125 m/s @ 60 deg angle'}])