% DEMO of using Matlab ODE solvers
%
% This script m-file solves for the initial velocity needed for a pojectile
% to hit a target a specified distance from the lauch point. The program
% includes air drag and height dependent gravity.
%
% There are many ways to solve this problem and we will use the Matlab ODE
% solvers with event detection (the event in this case is the z-value
% hitting zero).
%
global targdist g0 dgdh rho0 scaleht mass xarea cd
g0 = -9.806; % m/sec^2
rho0 = 1.29; % kg/m^3
%----- Get the inout that we need -----
prompt = {'Distance (km)', 'Error (m)', 'Launch Angle (deg)', 'Mass (kg)', ...
'CrossSection Area (m^2)', 'Drag Coefficient', ...
'AirScale height (m)', 'dg/ht (m/sec^2)/m'};
title = 'Ballistic Trajectory Parameters';
nlines = 1;
defaults = {'100','1','45','500','0.01','0.1','7500.e6', '3.0786e-6'};
% Get user inputs and save in the appropriate valuyes
result = inputdlg(prompt,title,nlines, defaults,'on');
%
targdist = eval(result{1})*1000.0;
terr = eval(result{2});
ltheta = eval(result{3})*pi/180;
mass = eval(result{4}) ;
xarea = eval(result{5}) ;
cd = eval(result{6}) ;
scaleht = eval(result{7}) ;
dgdh = eval(result{8}) ;
%
% Based on no-drag, no dght get initial velocity
vx = sqrt(-targdist*g0/(2*tan(ltheta)));
vz = vx*tan(ltheta);
vtot = vx/cos(ltheta);
dvxdD = sqrt(-g0/(2*tan(ltheta)*targdist))/2;
dvdD = dvxdD/cos(ltheta);
tmax = targdist/vx*2*(1+cd);

% Set the ODE options
options = odeset('AbsTol',[terr 1 1 1],'Events','Lec18_hit');

%%% Loop until target distance is hit
done = 0; prevDD = 0; iter = 0;
while done == 0
% Set the initial conditions for the ODE
iter = iter + 1;
y0 = [0.0; 0.0; vx; vz];
[t,y,te,ye,ie] = ode23(@Lec18_bacc,[0:1:tmax],y0,options);
% Check the time
if te/tmax > 0.5 , tmax = 2*te; end
% See how close we got to target
dd = targdist - ye(1);
fprintf('Result te %f ye %f dd %f \n',te, ye(1), dd);
if abs(dd) < terr
done = 1; % We are finished
else
% Adjust the initial velocity. Initially we use the analytic value
% for the directivative but this updated once we start
if prevDD == 0
vtot = vtot + dvdD*dd;
prevDD = dd;
else
dvdD = dvdD*(1+dd/prevDD);
vtot = vtot + dvdD*dd;
end
vx = vtot*cos(ltheta); vz = vtot*sin(ltheta);
end
if vtot > 3.e8
fprintf('\nVelocity faster than speed of light!\n')
done = 1;
end

end
%
% Output the final results
fprintf('TARGET DISTANCE: %8.3f km, within %6.1f m, Launch Angle %6.3f deg\n', ...
targdist/1000, terr, ltheta*180/pi)
fprintf('MASS %7.1f kg, XArea %7.4f m**2, Drag Coefficient %7.3f\n', ...
mass, xarea, cd)
fprintf('GRAVITY %8.5f m/s**2 %11.4e (m/sec**2)/m\n', ...
g0, dgdh)
fprintf('RHO AIR %8.5f kg/m**3, Scale Height %11.5e m\n', ...
rho0, scaleht)

fprintf('Initial Velocity needed %8.3f n/s\n',vtot)
fprintf('Calculation took %d iterations with final error %8.2f m\n',iter,dd)

% Now animate (make as function)
Lec18_animate(t,y);
fprintf('To show animation: animate(t,y)\n');
%