%% Numerical Integration of Coupled Differential Equations for Quadratic
%% Drag Force
% ODE45 is a function for the numerical solution of ordinary
% differential equations.  It employs variable step size Runge-Kutta
% integration methods.  ODE45 uses a 4th and 5th order pair of formulas
% for high accuracy.  This demo shows its use.
%
% Dale Gary 
% $Revision: 1.0 $  $Date: 2008/09/01 13:49:00 $

%%
% Consider the pair of first order ordinary differential equations for
% motion of a projectile (baseball) under the quadratic drag force
%
%     vdot_x = -(c/m)*sqrt(v_x^2 + v_y^2)*v_x
%
%     vdot_y = -g - (c/m)*sqrt(v_x^2 + v_y^2)*v_y
%
% The functions v_x and v_y are the horizontal and vertical velocities of the
% baseball, respectively.  The cross term accounts for the interations between
% the velocities.

%%
% To simulate a system, create a function M-file that returns a column vector of
% velocity derivatives, given velocity and time values.  For this example, we've
% created a file called QUAD_DRAG.M.

type quad_drag

%%
% To simulate the differential equation defined in QUAD_DRAG over the interval 0 < t
% < 8, invoke ODE45.  Use the default relative accuracy of 1e-3 (0.1 percent).

% Define initial conditions.
t0 = 0;          % Initial time (0 s)
tfinal = 8;      % Final time (8 s)
theta = 50*pi/180;         % Initial angle (converted to radians)
v_x = 30*cos(theta);       % Initial horizontal velocity (for a start angle of 50 degrees)
v_y = 30*sin(theta);       % Initial vertical velocity (for a start angle of 50 degrees)
y0 = [v_x; v_y];           % Create initial velocity value array  
% y0 = [0; 0]              % Uncomment for case of dropping a baseball off a building
% Simulate the differential equation.
tfinal = tfinal+eps;
[T,V] = ode45('quad_drag',[t0 tfinal],y0);

%%
% Plot the result of the velocity simulation two different ways.

subplot(1,2,1)
plot(T,V)
title('Time history')
xlabel('Time [s]')
ylabel('Velocity [m/s]')

subplot(1,2,2)
plot(V(:,1),V(:,2))
title('Velocity Vy vs. Vx plot')
xlabel('Vx [m/s]')
ylabel('Vy [m/s]')

%%
% Now integrate the velocities as a function of time to determine the
% position.  This is a crude approximation assuming constant velocities in
% each time bin.

type int_yp

pos = int_yp(T,V);
figure
subplot(1,1,1)
plot(v_x*T,v_y*T-4.9*T.^2);  % Plot vacuum case
xlabel('X position [m]');
ylabel('Y position [m]');
title('Trajectory plot')
hold on
plot(pos(:,1),pos(:,2),'color','red');  % Overplot quadratic drag case
hold off


%displayEndOfDemoMessage(mfilename)