basketball3

Trajectory of ball w/ quadratic drag -- to be called by compareIntegration.m
by R.Sonnenfeld, New Mexico Tech Physics -- Version 1.7  10/30/2008.
USAGE basketball3(1,100,3);
newplot=1 means -- Make a new plot (0 means just plot on top of existing plot)
npts=100 means -- Divide time interval into 100 points
If solvertype==1, Use Euler method
If solvertype==2, use RK2 method for ODEs
If solvertype==3, use RK4 method
If solvertype==4, use Matlab built-in function ODE45

The Runge-Kutta notation here is the same as used by Weisstein in
Weisstein, Eric W. "Runge-Kutta Method." 
From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Runge-KuttaMethod.html
See also http://en.wikipedia.org/wiki/Runge-Kutta#Examples

================
SET INITIAL CONDITIONS

function basketball3(newplot,npts,solvertype)
%Trajectory of ball w/ quadratic drag -- to be called by compareIntegration.m
%by R.Sonnenfeld, New Mexico Tech Physics -- Version 1.7  10/30/2008.
%USAGE basketball3(1,100,3);
%newplot=1 means -- Make a new plot (0 means just plot on top of existing plot)
%npts=100 means -- Divide time interval into 100 points
%If solvertype==1, Use Euler method
%If solvertype==2, use RK2 method for ODEs
%If solvertype==3, use RK4 method
%If solvertype==4, use Matlab built-in function ODE45
%
%The Runge-Kutta notation here is the same as used by Weisstein in
%Weisstein, Eric W. "Runge-Kutta Method."
%From MathWorld--A Wolfram Web Resource.
%http://mathworld.wolfram.com/Runge-KuttaMethod.html
%See also http://en.wikipedia.org/wiki/Runge-Kutta#Examples
%
%================
%SET INITIAL CONDITIONS
theta=40;   %degrees
v0=25;      %m/s
vx0=v0*cosd(theta);
vy0=v0*sind(theta);
x0=0;  y0=0; 
z0=[x0 y0 vx0 vy0];
%================
%SET PARAMETERS
% Fd=(1/16)*pi*D^2*rho*v^2=c*v^2
global m g c
D=0.24; % m
rho=1.3; % kg/m^3
m=0.624;  % kg
g=9.8;  % m/s*s
c=(1/16)*pi*D^2*rho;  %

%================
%SET TIME RANGE
t0=0;
tf=3.3;
%================
%ACTUALLY SOLVE ODE
t=linspace(0,tf,npts);

switch solvertype
    case 1
    [t,z]=Euler(t,z0); 
    case 2
    [t,z]=RK2(t,z0);
    case 3
    [t,z]=RK4(t,z0);
    case 4
    [t,z]=ode45(@quaddrag,t,z0);
    otherwise
    Warning('Invalid solver type.  What are you, nuts?')
end
%Note: The initial conditions go in z0. z0 is an array, not just a single #
% For this problem (a second order ODE) z0 contains initial x, y positions and vx, vy
% velocities.
%================
%UNPACK THE SOLUTIONS INTO CONVENIENTLY NAMED COLUMN VECTORS
x=z(:,1); y=z(:,2);
vx=z(:,3); vy=z(:,4);
%================
c=0; %No drag
%ACTUALLY SOLVE ODE
switch solvertype
    case 1
    [tvac,zvac]=Euler(t,z0);
    case 2
    [tvac,zvac]=RK2(t,z0);
    case 3
    [tvac,zvac]=RK4(t,z0);
    case 4
    [vact,zvac]=ode45(@quaddrag,t,z0);
    otherwise
    Warning('Invalid solver type.  What are you, nuts?')
end

%================
%UNPACK THE SOLUTIONS INTO CONVENIENTLY NAMED COLUMN VECTORS
xvac=zvac(:,1); yvac=zvac(:,2);
%================
%PLOT IT
if newplot==0
    hold on
    plot(x,y,'rx',xvac,yvac,'r+','MarkerSize',14);
    hold off
end

if newplot==1
    plot(x,y,'k-',xvac,yvac,'k--');
    xlabel('X-position (m)','Fontsize',14);
    ylabel('Y-position (m)','Fontsize',14);
    h=legend('Range in Air','Range in a vacuum');
    set(h,'Fontsize',12)
    axis([0 max(xvac) 0 20])
    xtx=8; ytx=16;  %xtx and ytx are variables to hold the location of the
    %text annotations on the plot.
    text(xtx,ytx-2,'\rho_{air} = 1.3kg/m^3', 'Fontsize', [16])
    text(xtx,ytx,'D=0.24m m=0.624kg ', 'Fontsize', [16])

end
end

function [t,z]=Euler(t,z0)
global m g c 
b=c/m;
len=length(t); 
z=zeros(len,4); %Set up a matrix for x, y, vx, vy
x=zeros(len,1);y=zeros(len,1);vx=zeros(len,1);vy=zeros(len,1);
%Plug initial conditions into convenient variables
x(1)=z0(1);y(1)=z0(2);
vx(1)=z0(3);vy(1)=z0(4);
for n=1:len-1
    dt=t(n+1)-t(n);
    %k1
    magv1=sqrt(vx(n)^2+vy(n)^2);
    k1vx=dt*(-b*vx(n)*magv1);
    k1vy=dt*(-g-b*vy(n)*magv1);
    k1x=dt*vx(n);
    k1y=dt*vy(n);
   
    vx(n+1)=vx(n)+k1vx;
    x(n+1)=x(n)+k1x;
    vy(n+1)=vy(n)+k1vy;      
    y(n+1)=y(n)+k1y;    
end
%keyboard
z(:,1)=x;  %1st column
z(:,2)=y;  %2nd column
z(:,3)=vx; 
z(:,4)=vy;
end

function [t,z]=RK2(t,z0)
global m g c
 
b=c/m;
len=length(t); 
z=zeros(len,4); %Set up a matrix for x, y, vx, vy
x=zeros(len,1);y=zeros(len,1);vx=zeros(len,1);vy=zeros(len,1);
%Plug initial conditions into convenient variables
x(1)=z0(1);y(1)=z0(2);
vx(1)=z0(3);vy(1)=z0(4);
for n=1:len-1
    dt=t(n+1)-t(n);
    %k1
    magv1=sqrt(vx(n)^2+vy(n)^2);
    k1vx=dt*(-b*vx(n)*magv1);
    k1vy=dt*(-g-b*vy(n)*magv1);
    k1x=dt*vx(n);
    k1y=dt*vy(n);
    %k2
    magv2=sqrt((vx(n)+k1vx/2)^2+(vy(n)+k1vy/2)^2);
    k2vx=dt*(-b*(vx(n)+k1vx/2)*magv2);
    k2vy=dt*(-g-b*(vy(n)+k1vy/2)*magv2);
    k2x=dt*(vx(n)+k1vx/2);
    k2y=dt*(vy(n)+k1vy/2);
 
    vx(n+1)=vx(n)+k2vx;
    x(n+1)=x(n)+k2x;
    vy(n+1)=vy(n)+k2vy;      
    y(n+1)=y(n)+k2y;    
end
%keyboard
z(:,1)=x;  %1st column
z(:,2)=y;  %2nd column
z(:,3)=vx; 
z(:,4)=vy;
end

function [t,z]=RK4(t,z0)
global m g c
% Note: The initial conditions go in z0.  z0 is an array, not just a single #
% For this problem (a second order ODE)
% z0 contains initial x, y positions and vx, vy velocities.
% The array t is not assumed to be evenly spaced (but it can be);
% The first returned data point corresponds to t(1) and the last
% to t(length(t)).
b=c/m;
len=length(t); 
z=zeros(len,4); %Set up a matrix for x, y, vx, vy
x=zeros(len,1);y=zeros(len,1);vx=zeros(len,1);vy=zeros(len,1);
%Plug initial conditions into convenient variables
x(1)=z0(1);y(1)=z0(2);
vx(1)=z0(3);vy(1)=z0(4);
for n=1:len-1
    dt=t(n+1)-t(n);
    %k1
    magv1=sqrt(vx(n)^2+vy(n)^2);
    k1vx=dt*(-b*vx(n)*magv1);
    k1vy=dt*(-g-b*vy(n)*magv1);
    k1x=dt*vx(n);
    k1y=dt*vy(n);
    %k2
    magv2=sqrt((vx(n)+k1vx/2)^2+(vy(n)+k1vy/2)^2);
    k2vx=dt*(-b*(vx(n)+k1vx/2)*magv2);
    k2vy=dt*(-g-b*(vy(n)+k1vy/2)*magv2);
    k2x=dt*(vx(n)+k1vx/2);
    k2y=dt*(vy(n)+k1vy/2);
    %k3
    magv3=sqrt((vx(n)+k2vx/2)^2+(vy(n)+k2vy/2)^2);
    k3vx=dt*(-b*(vx(n)+k2vx/2)*magv3);
    k3vy=dt*(-g-b*(vy(n)+k2vy/2)*magv3);
    k3x=dt*(vx(n)+k2vx/2);
    k3y=dt*(vy(n)+k2vy/2);
    %k4
    magv4=sqrt((vx(n)+k3vx)^2+(vy(n)+k3vy)^2);
    k4vx=dt*(-b*(vx(n)+k3vx)*magv4);
    k4vy=dt*(-g-b*(vy(n)+k3vy)*magv4);
    k4x=dt*(vx(n)+k3vx);
    k4y=dt*(vy(n)+k3vy);
    %%
    vx(n+1)=vx(n)+k1vx/6+k2vx/3+k3vx/3+k4vx/6;
    x(n+1)=x(n)+k1x/6+k2x/3+k3x/3+k4x/6;
    vy(n+1)=vy(n)+k1vy/6+k2vy/3+k3vy/3+k4vy/6;      
    y(n+1)=y(n)+k1y/6+k2y/3+k3y/3+k4y/6;    
end
z(:,1)=x;  %1st column
z(:,2)=y;  %2nd column
z(:,3)=vx; z(:,4)=vy;
end

function dz=quaddrag(t,z)
%THIS IS THE FUNCTION THAT THE ODE solver INTEGRATES.
%It is almost all that needs to be changed to do a different physics problem.
%x=z(:,1); y=z(:,2);
%vx_dot=z(:,3); vy_dot=z(:,4);

global m g c
%General equation. Standard vector form
% -->  mv_dot=-mg(y-hat)-c(v_x^2+v_y^2)v-hat
% Step A: Put in component form:
%   x_dot=vx
%   y_dot=vy
%   magv=(vx^2+vy^2)^0.5
%   m*vx_dot=-c*magv*vx;
%   m*vy_dot=-mg-c*magv*vy;
%
% Step B: Keep only derivatives on LHS:
% Divide through by m
%   x_dot=vx
%   y_dot=vy
%   magv=(vx^2+vy^2)^0.5
%   vx_dot=-c/m*magv*vx;
%   vy_dot=-g-c/m*magv*vy;
%
% Step C: Replace x, y, vx and vy by z1,z2,z3, and z4
%   x_dot=z3
%   y_dot=z4
%   magv=(x_dot^2+y_dot^2)^0.5
%   vx_dot=-c/m*magv*x_dot;
%   vy_dot=-g-c/m*magv*y_dot;
%
% Step D: Use these new subscripts to put in form of Matlab vectors
% In Matlab Vector Form as follows
% Fd=(1/16)*pi*D^2*rho*v^2=c*v^2

x_dot=z(3); %First derivative in x
y_dot=z(4); %First derivative in y
magv=(x_dot^2+y_dot^2)^0.5;
vx_dot=-(c/m)*x_dot*magv;  %Second deriv in X
vy_dot=-g-(c/m)*y_dot*magv; %Second deriv in Y

%Put them in a column matrix
dz=[x_dot;y_dot;vx_dot;vy_dot];
end