clear
reset(0)
echo on
%This program illustrates the use of the shooting
%method.  As an example, we shall look at the angle
%of elevation for a cannon necessary to hit a target.
%The shooting method is used to convert a boundary
%value problem into an initial value problem.  In
%the case of a cannon, we are following the trajectory
%of the cannon ball.  We know where it starts from
%(the location of the cannon) and its initial speed
%(determined by the amount of gunpowder used), but
%we have to figure out what angle of elevation is
%necessary to hit a target some distance away.  We
%define the vector y = [x,z,vx,vz] and write down
%evolution equations describing their change in the
%attached function cannonball.m.
pause
%We must specify the initial velocity:
v0=100;   %initial velocity (m/s)
%
%and the location of the target:
dist=200;   %distance (m)
height=10;   %elevation (m)
%and the discretization in time:
dt=.2;
%and we add some graphical bits:
figure(1)
plot(dist,height,'xg')
set(gca,'FontSize',18)
xlabel('x')
ylabel('height')
hold on
pause
%finally the initial angle in radians:
%(We'll use a manual input to start with)

theta=input('input elevation angle = ');
while theta~=0,
%We use this to calculate the velocities:
y(1)=0;
y(2)=0;
y(3)=v0*cos(theta);
y(4)=v0*sin(theta);
plot(y(1),y(2),'o')
set(gca,'FontSize',18)


%And now for the integration procedure itself.
%We shall use a 3-stage runge-kutta method:
t=0;
while (dist>y(1)&y(2)>=0),
 k1=dt*cannonball(t,y);
 k2=dt*cannonball(t+dt,y+k1);
 k3=dt*cannonball(t+dt/2,y+(k1+k2)/4);
 y=y+(k1+4*k3+k2)/6;
 t=t+dt;
 figure(1)
 plot(y(1),y(2),'o')
 set(gca,'FontSize',18)
 echo off
end;
theta=input('input elevation angle = ');
end;
echo on
hold off
%This program illustrates manual adjustment of the
%shooting parameter.  Manual adjustment is often
%quite effective if there is only one parameter to
%play with.  Certainly it allows you to get a good
%idea of the behavior of the function, and to get an
%estimate of the correct value of the shooting
%parameter.  A better way to do it, however, is to
%use a root finding technique.  We can employ the
%matlab routine fzero to find the angle which yields
%the correct height at the required distance.
pause
%To do this we must put the integration routine in
%a function which returns the difference between the
%height of the cannonball and the target at the distance
%of the target.  We will call this function deltah(theta).
%The matlab routine fzero will adjust theta until the
%value returned by this function is zero.
pause
plot(dist,height,'xg')
set(gca,'FontSize',18)
hold on
opttheta=fzero('deltah',0.1)
%and that's all there is to it!  Life gets a bit more
%complicated if you have two shooting parameters, but you
%can solve this as a 2-D minimization procedure where
%the function to be minimized is the sum of the square
%of the deviations.  If you do this, be sure to scale 
%everything carefully!
hold off
echo off


________________________


function f=cannonball(t,y)
%This function returns the derivative of x and z
%positions and velocity with time.  We have the
%parameter values:
m=5;   %Mass of cannon ball (kg)
cd=0.44;   %Drag coefficient
a=0.03;   %area of cannon ball (x-section, m^2)
rho=1.0;   %density of air (kg/m^3)
g=9.8;   %acceleration due to gravity (m/s^2)
%
%and the derivatives:
f(1)=y(3);
f(2)=y(4);
f(3)=-cd*a*rho/2/m*(y(3)^2+y(4)^2)^.5*y(3);
f(4)=-cd*a*rho/2/m*(y(3)^2+y(4)^2)^.5*y(4)-g;
%Note that the total drag is distributed between
%the two components of velocity proportional to the
%velocity in each direction.


________________________


function dh=deltah(theta)
%This function takes an input value of the
%initial angle and returns the difference in
%height between the cannonball and the target
%at the target range.  If this is zero, the
%target will be hit.  Note that after the
%integration procedure we have to do a bit of
%interpolation to get the height just right.

%We must specify the initial velocity:
v0=100;   %initial velocity (m/s)

%and the location of the target:
dist=200;   %distance (m)
height=10;   %elevation (m)
%and the discretization in time:
dt=.2;

%We have the initial positions:
y(1)=0;
y(2)=0;
y(3)=v0*cos(theta);
y(4)=v0*sin(theta);
plot(y(1),y(2),'o')

%And now for the integration procedure itself.
%We shall use a 3-stage runge-kutta method:
t=0;
while dist>y(1),
 ylast=y;
 k1=dt*cannonball(t,y);
 k2=dt*cannonball(t+dt,y+k1);
 k3=dt*cannonball(t+dt/2,y+(k1+k2)/4);
 y=y+(k1+4*k3+k2)/6;
 t=t+dt;
 plot(y(1),y(2),'o')
end;
drawnow
%Now for the interpolation bits.  The time when
%the cannonball passed the target is a little less
%than the final time t.  We can get at this
%difference:
delta=-(y(1)-dist)/(y(1)-ylast(1))*dt;
%and now we use a couple of terms in a Taylor series
%to get a more accurate measure of the height:
hfinal=y(2)+y(4)*delta+0.5*delta^2*(y(4)-ylast(4))/dt;

%Finally, we return the difference in height:
dh=hfinal-height;