% impact.m
% Finds the point of impact of a parabolic trajectory.
% J.M.Soler. Nov.2012

% Set terrain data
clear
h = 200;                % hill height (m)
a = 800;                % hill width (m)

% Set trajectory data
x0 = 0; y0 = h;         % initial position (top of hill)
v0 = 100;               % initial velocity (m/s)
theta = 30*pi/180;      % shot angle (rad)
g = 9.8;                % gravitational acceleration (m/s^2)

% Find terrain height
x = 0:10:2500;                   % horizontal mesh for plot (m)
yground = h*exp(-0.5*(x/a).^2);  % ground heights

% Find trajectory
vx = v0*cos(theta);     % horizontal velocity
vy = v0*sin(theta);     % initial vertical velocity
t = (x-x0)/vx;          % times at points x
y = y0+vy*t-0.5*g*t.^2; % vertical position at times t

% Plot trajectory and terrain
plot( x,yground, x,y )
legend('terrain','trajectory')
% axis equal            % uncomment for true aspect ratio
axis([0,2500,0,400])
xlabel('x (m)')
ylabel('y (m)')
grid on

% Set data for root search
x1 = 1;                 % min. distance to impact (m)
x2 = 5000;              % max. distance to impact (m)
dxMax = 1.e-3;          % required precision for impact point (m)

% Find point of impact (root of height function)
x_impact = bisectionRoot( @height, x1, x2, dxMax )   % simpler
x_impact = bisectionRoot( @(x)generalHeight(x,h,a,v0,theta), ...
                          x1, x2, dxMax )            % more general
y_impact = h*exp(-0.5*(x_impact/a).^2)
t_impact = (x_impact-x0)/vx

