echo off
clc
% need to use canode23.m which is set up to break out of main loop
% when y is negative
%
% consider the equation for the cannon ball from
%
% Scientific Computing and Differential Equations
% G.H. Golub and J. M. Ortega
% Academic Press, 1992
% p. 18
%
% (2.1.15) x' = v cos(theta)   y' = v sin(theta)
%
%                T - c p s v^2/2)                   mp
% (2.1.17) v' =  ---------------- - g sin(theta) - --- v
%                      m                            m
%
%                      g
% (2.1.18) theta' = - --- cos(theta)
%                      v
%
% (x,y) cartesian coordinates
%   v   velocity
% theta launch angle
%   m   mass of projectile
%   mp  m' the rate of change of mass of the projectile
%   T   thrust as a function of time
%        cpsv^2
%   D =  ------  drag force; c=coef of drag; p (rho) = air density
%          2                 s = cross-sectional area of projectile
%   g   acceleration of gravity
%
% let mp = 0
%     T = 0
% from Kahaner, Moler and Nash, problem 8-4, let
%     c = 0.2
%     p = 1.29 kg/meter
%     s = 0.25 meter^2
%     g = 9.81 meter/sec^2
%     m = 15 kg
%
% these values will be made global so that cannon.m has their values
%
c=0.2; r=1.29; s=0.25; g=9.81; D=c*r*s/2;
global g D;
%
% v(0) = 50 meter/sec
%
% z = [x y v theta]
angled = input(' enter initial launch angle (in degrees) ');
% Matlab works in radians, so convert degrees to radians
% Matlab has "special variables" like pi predefined.
angle = angled * pi / 180.0;

% don't need since canode23 has a break tfinal = input('enter final time');
tfinal = 100;

z0 = [0 0 50 angle];
t0 = 0;
tol = 1.e-5;
trace = 0;
[t,z] = canode23('cannon',t0,tfinal,z0,tol,trace);
%
mntx=[50 90 130]; mnty=[0 60 0];
plot(z(:,1),z(:,2),mntx,mnty),
stit=['Cannon shot    Launch angle=' num2str(angled) ' degrees;  Accuracy='num2str(tol)];
title(stit),xlabel('Distance'),ylabel('Height')