% Ballistics Project by Peter Howard and Jean Marie Linhart
% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
% 
% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}
\usepackage{gensymb}
\usepackage{verbatim}
\usepackage{float}
%
\renewcommand{\baselinestretch}{1.3}
\setlength{\emergencystretch}{3 em}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.
%
% Top matter is above and the beginning of the document below.
%
\begin{document}

\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote. 
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION\else STUDENT VERSION \fi\\Ballistics Modeling with a Sponge Dart}



\markboth{Ballistics}{Ballistics}


\author{Peter Howard \\
Texas A\&M University\\
College Station TX 78743 USA\\
phoward@math.tamu.edu\\
\\
Jean Marie Linhart \\
Central Washington University\\
Ellensburg WA 98926 USA\\
JeanMarie.Linhart@cwu.edu\\
}

\makeStitlePDFLaTex

\begin{abstract}
The goal of this project is for students to develop, analyze, and compare three different models for the flight of a sponge dart moving under the influences of gravity and air resistance.  The first two models are based respectively on the common simplifying assumptions of no air resistance and linear air resistance, and the third is based on the correct physical form of air resistance (which students are asked to identify).  Data is provided both for parameter estimation and for model verification, and the verification data is used to compare results of the models.  While most calculations  associated with the first two models can be carried out by hand, technology is essential for analysis of the physical model.
\end{abstract}


\section*{SCENARIO DESCRIPTION}

Ballistics is the science of projectile motion and impact, phenomena well described by Newtonian mechanics.  The number of applications of this type of analysis is staggering, ranging from such mundane issues as automobile accident simulations and optimal golfing to the critical studies of missile defense and space exploration.  Somewhat less dramatically, in this project we will use Newtonian mechanics to describe the flight of a sponge dart (a Nerf dart or an off-brand), light enough so that air resistance will play a critical role.

\section*{EXPERIMENTAL DATA}

The data for this project was collected by firing sponge darts from a toy gun (\$3.99, WalMart---it's not all that easy talking a mathematics department into financing this sort of thing). 

\begin{figure}[h!]
\begin{center}
\begin{tabular}{p{3.1 in} p{3.1in}}
\includegraphics[width=3.in]{NerfGuns.JPG} & 
\includegraphics[width=3. in]{Protractor2.jpg}\\
\end{tabular}
\parbox{9.4 true cm}{\caption{\label{fig:NerfGuns} Left: Sponge dart gun and dart. Right: A plastic protractor used for measuring angles.}}
\end{center}
\end{figure}

The angle of the gun and muzzle height were set and the gun was
fastened with duct tape in the correct position to an upright stand
for firing.  For each dart fired, an observer marked the spot where
the dart initially hit the ground, and a tape measure was used to
obtain the horizontal distance to this point from the point where the
gun was fired.  The horizontal distance was measured four times and
those four readings were averaged to make the data point.  Likewise, a
timer was used to make four measurements of the amount of time it
takes a dropped dart to reach the ground, and how long it takes a dart
fired vertically to reach the ground.  The average of these times was
recorded in the dataset included in Table \ref{table:data}.


Using this setup, first, we tabulated a set of measurements for
distance traveled (by the projectile) versus angle of inclination of
the gun, taking angles of inclination 5, 10, 15,
\ldots, 85 degrees.  The darts were fired from a height of 0.18
meters.

\begin{table}[h!]
\begin{center}
  \medskip
  \begin{tabular}{|c|c|c|}
    \hline
    Initial height (m) & Dropped/Fired & Time to reach ground (s) \\
    \hline
    4.06  & dropped & 0.95\\
    \hline
    0.39  & fired & 2.13 \\
    \hline
  \end{tabular}

  \medskip
\begin{tabular}{|l| c | c | c | c | c | c | c | c | c | c |}\hline
Angle of inclination (deg) & 5 & 10 & 15 & 20 & 25 & 30 & 35 & 40 & 45 \\
\hline
Distance traveled (m) & 4.37 & 5.23 & 6.95 & 7.84 & 8.17 & 8.69 & 8.81 & 8.99 & 8.95\\
\hline
\end{tabular}

\medskip
\begin{tabular}{|l| c | c | c | c | c | c | c | c | c | }\hline
Angle of inclination (deg) & 50 & 55 & 60 & 65 & 70 & 75 & 80 & 85 \\
\hline
Distance traveled (m) & 8.83 & 8.19 & 7.84 & 7.12 & 6.38 & 5.08 & 3.34 & 2.13 \\
\hline
\end{tabular}

\medskip
\parbox{10.7 true cm}{\caption{\label{table:data} The top table contains distance (meters) and time (seconds) data for a dropped dart and for a dart fired straight up. The second and third table contain angle of inclination in degrees and distance traveled (meters) for the sponge darts fired from a height of 0.18 meters.}}
\end{center}
\end{table}

These fairly simple measurements will suffice for the assignments in this project.

For most of the assignments, estimation of model parameters will be based on the falling dart and on the dart fired straight up.  The data in the table of angles of inclination and distances traveled will only be used to evaluate the models.  The final assignment will address incorporating the entire collection of data in parameter estimation and evaluating the resulting model.

For those wishing to replicate the original experiment, two
people are needed to fire the gun and observe and record measurements.  A
ruler or other straight edge can be taped to the gun to facilitate
angle measurements.  It is possible to use video analysis tools such
as \cite{tracker}, to analyze videos of the dart's
motion.  These tools often require having an object of a known length
such as a meter stick in the video image.




\section*{ASSIGNMENTS}

\subsection*{Analysis in the absence of air resistance}

Ignoring air resistance, Newton's equations of motion for an object
under the influence of gravity alone are quite straightforward.
Letting $y(t)$ represent the height of the object as a function of
time and $x(t)$ represent its horizontal distance traveled as a
function of time, we have\footnote{We are also ignoring the fact that
the gravitational pull on an object above the earth depends on its
height above the earth.  We take $g=9.81$m/s$^2$.}
\begin{equation}
\dfrac{d^2}{dt^2} y(t) = -g, \qquad \dfrac{d^2}{dt^2} x(t) = 0\label{eq:grav}
\end{equation}

Your tasks are:
\begin{enumerate}
\item Use the time for a dart fired straight up to hit the ground (in Table \ref{table:data}) to determine the initial velocity at which the darts are fired.
\item Solve the two second order differential equations in \eqref{eq:grav} for $x(t)$ and $y(t)$ and use a computer, e.g. Python or MATLAB, to plot the trajectory ($x$ versus $y$) obtained for angle of inclination $\theta = 35\degree$.  Observe that $x(0)$ and $y(0)$ were determined by the experimental set-up and are known to be $0$ and $0.18$ m respectively.  The remaining two initial conditions $x'(0)$ and $y'(0)$ are given by
$$x'(0) = v_0\cos \theta,\qquad \mbox{and}\qquad y'(0) = v_0\sin\theta,$$
where $v_0$ is the speed with which the gun fires its darts.  Make sure the aspect ratio on this plot is 1.
\item Determine the distance, $d(\theta)$, your object travels as a function ofthe firing angle, $\theta$, and  compute the angle that maximizes its distance.  You may use a computer to help compute the angle that maximizes the distance traveled by the dart. Your answer should be in degrees and accurate to two decimal places.  
\item Plot $d(\theta)$ using a computer, e.g. Python or MATLAB, and compare it with the experimental values $d_{\mbox{exp}}(\theta)$ given in Table \ref{table:data}.  Discuss the discrepancies.
\item Compute an error, $E$, for your model based on the sum of squared errors,
$$E = \sqrt{\sum_{k=1}^{17}\left(d_{\mbox{exp}}(\theta_i) - d(\theta_i)\right)^2}$$
where Table \ref{table:data} contains data $\theta_i$ and $d_{\mbox{exp}}(\theta_i)$.
\end{enumerate}


\subsection*{Linear Air Resistance}

A typical first attempt at introducing the effects of air resistance into a physical model is through linear air resistance, using terms of the form $-by'(t)$ and $-bx'(t)$.  Adding such terms to \eqref{eq:grav}, we obtain

\begin{equation}
\dfrac{d^2}{dt^2} y(t) = -g -b\dfrac{d}{dt} y(t), \qquad \dfrac{d^2}{dt^2} x(t) = -b \dfrac{d}{dt} x(t), \label{eq:linresist}\end{equation}
with the same initial conditions as before.  We will now work through steps similar to those in the case of no air resistance, except that as the mathematics becomes more cumbersome, we will begin to fall back more on computer tools, e.g. Python or MATLAB.

Your tasks are:
\begin{enumerate}
  \item Determine  the units on constant $b$.  Check that the signs are
  correct so that the air resistance term always works to bring the
  velocity closer to zero.
\item Develop analytic solutions to  \eqref{eq:linresist}.
\item Use your solutions to \eqref{eq:linresist} and the experimental data from Table \ref{table:data} to determine the value of $b$ for this model.
\item Find a (corrected) initial velocity $v_0$.
\item Plot the trajectory ($y$ versus $x$) for angle $\theta = 35 \degree$ along with your similar trajectory in the absence of air resistance. Make sure the aspect ratio of the plot is 1.
\item Analytically, determine an implicit equation for $d(\theta)$, the distance the dart travels given that it was shot with angle of inclination $\theta$. Use a computer to help you determine the value of $\theta$ that maximizes $d(\theta)$.
\item Plot $d(\theta)$ for this model along with the data, and $d(\theta)$ from the model with no air resistance on the same figure, and again discuss the discrepancies.
\item Compute the error $E$ for this model, and compare it to that from the model with no air resistance.
\end{enumerate}

\subsection*{Physical Air Resistance}

In general, linear air resistance is not an adequate model of projectile motion.  In this part of the project, you will work through the steps of the previous cases with a more physical form of air resistance. Unfortunately, analytic solutions for your model will be too cumbersome to derive or work with, so you will have to do most of your analysis with a computer, e.g. Python or MATLAB.

Your tasks are:

\begin{enumerate}
\item Use the method of dimensional analysis to determine a more physical form for your air resistance.  Give new differential equations with the new physical air resistance term.
  \item Check your work. Make sure the air resistance term always works to decrease the velocity.  Once again you use an unknown constant $b$ in the air resistance term; what are its units?
\item Use a computer  and the experimental data from Table \ref{table:data} to determine the value of $b$ for your model.
\item Use your new value for $b$ to determine a (further corrected) initial velocity $v_0$.
\item Use Python's {\tt SciPy.integrate.solve\_ivp()} or MATLAB's {\tt ode45()} and {\sl event location} to plot a trajectory for angle $\theta = 35\degree$, and plot it along with the trajectory without air resistance and with the trajectory with linear air resistance at this angle.
\item Plot $d(\theta)$ for this model along with the data and $d(\theta)$ with no air resistance and with linear air resistance.  Use a computer to help you determine the value of $\theta$ that maximizes $d(\theta)$ for this model.
\item Compute your error $E$ for this model and compare it with your errors
from the model without air resistance and the model with linear air resistance.
\end{enumerate}

\subsection*{Regression Approach}

If we have the same amount of confidence in each data point, we expect to obtain better values for $b$ and $v$ by using more of the data in Table \ref{table:data} to estimate these parameters.

Your tasks are:

\begin{enumerate}
\item Using the values of $b$ and $v_0$ from the physical air resistance model as starting values, carry out a nonlinear least squares regression fit for the function $D(\theta, b, v_0)$ taken to return the horizontal distance traveled by a dart launched with angle of inclination $\theta$, initial velocity $v_0$, and coefficient of air resistance $b$.  You will obtain further corrected values for $b$ and $v_0$.
\item Plot your fit for $D(\theta, b, v_0)$ along with your plots of $d(\theta)$ from the physical air resistance model.
\item Compute the error $E$ for the  model with these parameter values, and compare to the physical air resistance model.
\end{enumerate}

Wikipedia provides a good general overview of projectile motion, covering many of the ideas covered in this project \cite{Wikipedia}.


\begin{thebibliography}{98}
\bibitem{Wikipedia} Wikipedia Contributors. 2018. Projectile Motion. {\sl Wikipedia, The Free Encyclopedia\/}. Available from \href{https://en.wikipedia.org/wiki/Projectile\_motion}{https://en.wikipedia.org/wiki/Projectile\_motion}.
\end{thebibliography}

\teacher{
  \section*{NOTES FOR TEACHERS}
  This project was originally created for incorporation in a senior-level mathematics course in mathematical modeling at Texas A\&M University in College Station, Texas.
Parts of this project are taught in a quarter-long course at Central Washington University.
Instructors are encouraged to incorporate parts of this project into their own modeling and differential equations courses. The calculations for the trajectories $x(t)$ and $y(t)$ in the no air resistance and linear air resistance models, in particular, give a review of concepts and hand calculation techniques that are usually taught in introductory differential equations courses.  Learning to write the computer code for these two models is an excellent introduction to basic programming for mathematics students.

Developing the model for physical air resistance is a practical
application of concepts of multivariate calculus and vectors with
differential equations.  The programming for this model is highly
challenging and can provide a way for students to show or to develop
algorithmic thinking, while also learning how to use
computer tools to solve a complicated problem. Alternatively, an
instructor could provide results for this model once developed rather
than taking students step-by-step into the tools and methods to
develop their own code.  This is the approach taken at Central Washington University.

Students can collect their own data for all or some of this project
using simple equipment: A sponge dart gun (a Nerf gun or an off-brand),
several tape measures, a
large protractor (inclinometer), and a timer.  The original data was
collected with such equipment.  It is also possible to use a cell
phone to take video of the dart with some known length measures in the
video frame.  A cell phone video of the dart firing and its trajectory
can be analyzed frame-by-frame to find the distance the dart traveled
before hitting the ground and also the time from firing until it hit
the ground.

Our modern notions of ballistics date to the works of Galileo
\cite{Galileo}. Projectile motion, with and without air resistance,
gets only a brief treatment in many differential equations books
(\cite{BoyceDiPrima}, \cite{Edwards}, \cite{Nagle}, \cite{Noonburg}).
A more thorough treatment can be found in books on computational
physics such as \cite{Giordano} and \cite{Malthe}.  The reference
\cite{Malthe} also provides an introduction to programming with
Python.

Python 3 code for this project was created using Spyder, and is
provided in the file \linebreak[4] {\tt Ballistics.py}.  The popular
Anaconda distribution containing Python, SciPy and NumPy is freely
available \cite{anaconda}.  MATLAB code is in a variety of program
files; we discuss the relevant files and what is calculated in each at
the end of each section. The MATLAB Optimization toolbox is needed for
many of the optimization routines mentioned here, and the Symbolic
Mathematics Toolbox is needed for one program, although an alternative
implementation that avoids it is also provided. Numerical results
agree on the two platforms to three significant figures. We will
report all calculated values to 4 significant figures.  Results
mentioned in the discussion are from the Python code.

\subsection*{Analysis in the absence of air resistance}

Solutions to the differential equations in \eqref{eq:grav}
can be obtained by integrating twice with respect to $t$.  The initial
speed of the dart is $v_0$ and $\theta$ is the angle at which the dart is
fired.

\begin{equation}
  y(t) = -\frac{1}{2}gt^2 + v_0t\sin\theta  + H \qquad x(t) = v_0t\cos\theta. \label{eq:gravsoln}
\end{equation}

When the dart is fired straight up ($\theta = \frac{\pi}{2}$) from a height $H = 0.39$ meters (m) it takes $t = 2.13$ seconds (s) to hit the ground. Applying this information with the equation for $y(t)$ in \eqref{eq:gravsoln} and setting $y(2.13) = 0$, allows one to solve for $v_0 = 10.26$ m/s.  With known values for $v_0$, $\theta$, and $H$,  we can solve $y(t_f) = 0$, a quadratic, for $t_f(\theta)$ the time it takes a dart to hit the ground, in terms of $\theta$.  
\begin{equation}
t_f(\theta) = \dfrac{v_0\sin\theta + \sqrt{(v_0\sin\theta)^2 + 2gH}}{g}.\label{eq:tf}
\end{equation}
The value of $t_f(\theta)$ can be used to calculate $x(t_f(\theta)) = d(\theta)$ which gives the distance the dart travels as a function of $\theta$
\begin{equation}
 d(\theta) = v_0\cos\theta\left[\dfrac{v_0\sin\theta + \sqrt{(v_0\sin\theta)^2 + 2gH}}{g}\right]. \label{eq:dnoair}
\end{equation}
The leftmost plot in Figure \ref{fig:NoAirDistAngle} shows the trajectory when $\theta = 35^\circ$; note that the axis scaling is the same for both the $x$ and $y$ axis in this graph.  The rightmost plot compares the data for the distance the dart traveled to the model predictions.



To determine the value of $\theta$ that gives the longest trajectory, differentiate \eqref{eq:dnoair} to find $d'(\theta)$ and then find the value of $\theta$ for which $d'(\theta) = 0$.  Calculating the derivative can be done
by hand, or with the Symbolic Math Toolbox in MATLAB (or with Mathematica). 
% It may also be possible to do this with SymPy, a Python library for doing symbolic mathematics.
\begin{equation}
d'(\theta) = -v_0\sin\theta\left[ \dfrac{v_0\sin\theta + 
\sqrt{(v_0\sin\theta)^2 + 2gH}}{g}\right]+ v_0\cos\theta\left[\dfrac{v\cos\theta}{g} + 
\dfrac{v_0^2\sin\theta\cos\theta}{g\sqrt{(v_0 \sin\theta)^2 + 2gH}}\right].
\label{eq:dprime}
\end{equation}
\begin{figure}[t]
  \begin{center}
 \begin{tabular}{p{3. in} p{3. in}}
    \includegraphics[width=2.7 in]{NoAirTrajectory35.pdf} &
    \includegraphics[width=2.7 in]{NoAirDistAngle.pdf}\\
    \end{tabular}
\parbox{10.7 true cm}{ \caption{\label{fig:NoAirDistAngle} Left: foam dart trajectory for a firing angle of $\theta = 35^\circ$. Right: distance traveled vs. firing angle in the case of no air resistance.}}
\end{center}
\end{figure}
 It is difficult to solve $d'(\theta)=0$ for $\theta$ by hand, but you
 can use Python's {\tt fsolve()} or MATLAB's {\tt fzero()}
 routine to find the value of $\theta$ for which $d'(\theta) = 0$.
 The routines {\tt fsolve()} and {\tt fzero()} are nonlinear equation
 solvers; both solve an equation of the form $y(t) = 0$ for $t$.  Note
 that if, instead, you wanted to solve $y(t) = c$ you could define a
 new function $\tilde{y}(t) = y(t) - c$ and solve $\tilde{y}(t) = 0$
 to find the value of $t$ that gives $y(t) = c$. We found $\theta = 44.53^\circ$, slightly under $45^\circ$
 because the dart was not fired from ground level.  

The square root of the sum of squared errors for this model is $E = 5.099$ m. This number is truly only meaningful to compare two or more models, such as the model with no air resistance and those with air resistance.

The MATLAB code to plot the 35$^\circ$ trajectory (and the one for linear
air Resistance and physical Air Resistance) is in {\tt
  ballistics\_trajectories.m}.  A comparison of model distances and
actual distances versus firing angle (again for all 3 models) is in
{\tt ballistics\_distances.m}; this code uses the function in {\tt
  dist.m} to do some calculations for the physical air resistance
model. Code to determine the firing angle for the longest trajectory
is in {\tt noair1.m} (using MATLAB's Symbolic Toolbox) and {\tt
  noair2.m} (programming the derivative in \eqref{eq:dprime}
directly).

\subsection*{Linear Air Resistance}

We solve the differential equation for $\dfrac{d^2}{dt^2} y(t)$ in \eqref{eq:linresist} by integrating twice.  Let $w(t) = y'(t)$, then \eqref{eq:linresist} becomes

\begin{equation}
  w'(t) + b w(t) = -g. \label{eq:linresista}
\end{equation}
We solve by using an integrating factor $\mu(t)$. We multiply through by $\mu(t)$ to put the left hand side in the form $\left[\mu(t)w(t)\right]'$. We find
$\mu'(t) = b \mu(t)$ which has solution $\mu(t) = e^{bt}$.  \eqref{eq:linresista} becomes
\begin{equation}
  \left[e^{bt}w(t)\right]' = -ge^{bt}. \label{eq:linresistb}
\end{equation}
Integrating both sides, then solving for $w(t)$, we find
\begin{equation}
  y'(t) = w(t) = -\frac{g}{b} + Ce^{-bt}.
\end{equation}
Since $w(0)=y'(0) = v\sin\theta$, we find $C =v_0\sin\theta + \frac{g}{b}$, and we have
\begin{equation}
  y'(t)  = -\frac{g}{b} + \left(v_0\sin\theta + \frac{g}{b}\right)e^{-bt}.\label{eq:this}
\end{equation}
Integrating \eqref{eq:this} and using $y(0) = H$ we see that
\begin{equation}
  y(t) = H -\frac{g}{b}t + \left(\frac{v_0\sin\theta}{b} + \frac{g}{b^2}\right)(1 - e^{-bt}).\label{eq:LinAirY}
\end{equation}

Integrating \eqref{eq:linresist} with initial conditions $x'(0) = v_0\cos\theta$ and $x(0) = 0$ to find $x(t)$ is more straightforward.  First, separate variables,
\begin{equation}
\dfrac{dx'}{x'} = -b\, dt.
\end{equation}
Integrate and solve for $x'(t)$ and use $x'(0) = v_0 \cos\theta$ to find the unknown constant of integration,
\begin{equation}
  x'(t) = (v_0\cos\theta) e^{-bt}.
\end{equation}
Now integrate again, using $x(0) = 0$ to solve for the unknown constant of integration.  We find
\begin{equation}
  x(t) = \frac{v_0\cos\theta}{b} -\frac{v_0\cos\theta}{b} e^{-bt}.
\label{eq:LinAirX}
\end{equation}

To find $b$, we use \eqref{eq:LinAirY} and the data from the dropped sponge dart ($v_0 = 0$, $H = 4.06$, $t_f = 0.95$, $\theta = \frac{\pi}{2}$) with {\tt fsolve()} or {\tt fzero()} to find $b = 0.2792$. Next, we use equation
\eqref{eq:LinAirY} with the data from the sponge dart fired straight up ($H = 0.39$, $t_f = 2.13$, $\theta = \frac{\pi}{2}$) with {\tt fsolve()} or {\tt fzero()} to find the initial velocity of the gun, $v_0 = 11.23$ m/s.  Finally, we use {\tt fsolve()} or {\tt fzero()} to find the value for $t_f$ from $y(t_f) = 0$ for a given firing angle $\theta$.  Substituting this value for $t_f$ in \eqref{eq:LinAirX} yields the distance the dart travels for a firing angle $\theta$. This allows us to create the graphs in Figure \ref{fig:LinAirDistAngle}.


\begin{figure}[t]
  \begin{center}
\begin{tabular}{p{3. in} p{3. in}}
    \includegraphics[width=2.7 in]{LinTrajectory35.pdf} &
    \includegraphics[width=2.7 in]{LinDistAngle.pdf}\\
    \end{tabular}
\parbox{10.3 true cm}{ \caption{\label{fig:LinAirDistAngle} Left: foam dart trajectory for a firing angle of $\theta = 35^\circ$ in the model with no air resistance and the model with linear air resistance. Right: distance traveled vs. firing angle for the model with no air resistance and the model with linear air resistance.}}
    \end{center}
\end{figure}

To find $d(\theta)$, the distance the dart travels as a function of the firing angle, we use \eqref{eq:LinAirY} and \eqref{eq:LinAirX}, and let $t(\theta)$ be the time the dart hits the ground as a function of $\theta$.  Then $x(t(\theta)) = d(\theta)$, and
\begin{equation}
d(\theta) = \dfrac{v_0}{b} \cos\theta\left(1-e^{-bt(\theta)}\right).
\end{equation}
This implies
\begin{equation}
\left(1-e^{-bt(\theta)}\right) = \dfrac{bd(\theta)}{v_0\cos\theta}, \label{eq:LinAirD1}
\end{equation}
and
\begin{equation}
t(\theta) = -\dfrac{1}{b} \ln\left(1-\dfrac{bd(\theta)}{v\cos\theta}\right).\label{eq:LinAirD2}
\end{equation}
Setting $y(t(\theta)) = 0$, and substituting \eqref{eq:LinAirD1} and \eqref{eq:LinAirD2} into the right hand side into \eqref{eq:LinAirY} gives an implicit equation for $d(\theta)$:
\begin{equation}
0 = H + \dfrac{g}{b^2}\ln\left(1 - \dfrac{bd(\theta)}{v_0\cos\theta}\right) +
\left(\dfrac{v_0\sin\theta}{b} + \dfrac{g}{b^2}\right) \dfrac{bd(\theta)}{v\cos\theta}.\label{eq:ImpD}
\end{equation}
We can use {\tt fsolve()} or {\tt fzero()} to find $d(\theta)$ from
\eqref{eq:ImpD}, and find $-d(\theta)$.  This can be minimized using
{\tt SciPy.optimize.fmin()} or {\tt fminsearch()} with MATLAB.  Using
{\tt fmin()}, we find the angle that maximizes $d(\theta)$ is
$40.68^\circ$.  Since the horizontal component of the velocity is decreasing from air resistance, the dart must be fired at a more shallow angle to give it more time in the air to achieve the maximum range.

The square root of the sum of squared errors for this model drops to $E= 2.753$ from $E=5.099$ in the model without air resistance.

MATLAB code to find $b$ and $v$ for linear air resistance is {\tt
linair1.m}.  Code to plot the 35$^\circ$ trajectory and the one for
no air resistance and physical air resistance is in {\tt
ballistics\_trajectories.m}.  A comparison of model distances and
actual distances versus firing angle (again for all 3 models) is in
{\tt ballistics\_distances.m}, this code uses the function in {\tt
dist.m} to do some calculations for the physical air resistance model. Code to find the value of $\theta$
that maximizes $d(\theta)$ for linear air resistance is in {\tt
linair2.m}.



\subsection*{Physical Air Resistance}
Letting $L$ represent length and $T$ represent time, a very elementary
dimensional analysis tells us that the air resistance term must have
units $LT^{-2}$, and this might be enough to suggest an attempt at a
more physical formulation of the equations in terms of the velocity
squared.

More formally, we can consider the force $F$ due to air resistance, and start with the ansatz that it depends on the air density $\rho$, surface area of the object $S$, and velocity of the object $v$,
\begin{equation*}
F = \kappa \rho^a S^b v^c,
\end{equation*}
where $\kappa$ is a dimensionless constant.  This form of this ansatz is justified by Buckingham's Theorem \cite{Buckingham}.  Now we equate dimensions in this relation, $M$ represents mass, $L$ represents length, and $T$ represents time, to get
\begin{equation*}
MLT^{-2} = M^a L^{-3a} L^{2b} L^c T^{-c},
\end{equation*}
and consequently the system
\begin{align*}
&L:\quad 1 = -3a+2b +c \\
&M:\quad 1 = a \\
&T: \quad -2 = -c.
\end{align*}
We find $a=1$, $c=2$ and $b=1$ giving 
\begin{equation}
F = k\rho S v^2.
\end{equation}
Ultimately, we let 
\begin{equation*}
b = \dfrac{\kappa \rho S}{m},
\end{equation*}
and include a term that is proportional to the velocity squared.


A first attempt to incorporate a term proportional to the velocity
squared might be
\begin{equation}
  \dfrac{d^2}{dt^2}y(t) = -g -b\left(\dfrac{d}{dt}y(t)\right)^2 \qquad \dfrac{d^2x}{dt^2} = -b \left(\dfrac{d}{dt}x(t)\right)^2.\label{eq:quadresist0}
\end{equation}
There are several problems with this attempt.  First, $\left(\frac{d}{dt}y(t)\right)^2$ and $\left(\frac{d}{dt}x(t)\right)^2$ are always positive, and so the air resistance term can act in the same direction as the velocity.  Second, more importantly, this should be a vector model; position, velocity, acceleration, and force should all be understood as vectors.  In particular, the air resistance should act in the opposite direction of the velocity vector, and the above equations {\it do not} take this into account.  A better formulation is
\begin{equation}
  \dfrac{d^2}{dt^2}\vec{r}(t) = -\vec{g} - b\,\|\vec{r}\,'(t)\|\,\vec{r}\,'(t)\label{eq:quadresist}
\end{equation}
where $\vec{r}(t) = \begin{bmatrix} x(t) \\ y(t) \end{bmatrix}$ and
$\vec{g} = \begin{bmatrix} 0 \\ g \end{bmatrix}$.  Initial conditions
are $\vec{r}(0) = \begin{bmatrix} x(0) \\ y(0) \end{bmatrix}
= \begin{bmatrix} 0 \\ H \end{bmatrix}$, where $H$ is the height
at which the dart is fired, and $\vec{r}\,'(0) = \begin{bmatrix} x'(0)
  \\ y'(0) \end{bmatrix} = \begin{bmatrix} v_0\cos\theta
  \\ v_0\sin\theta \end{bmatrix}$, where $v_0$ is the initial speed of
the dart. We now have a coupled system of differential equations and
analytic solution is either impossible or very cumbersome, and we rely
on technology to solve and complete the analysis of this model.

First, we translate \eqref{eq:quadresist} in to a linear first-order system in components, where $R_1(t) = x(t)$ and $R_2(t) = \dfrac{dx}{dt}(t)$, $R_3(t) = y(t)$ and $R_4(t) = \dfrac{dy}{dt}(t)$:
\begin{equation}
\begin{aligned}
\dfrac{d R_1}{dt} &= R_2 \\
\dfrac{d R_2}{dt} &=-b\left[\sqrt{(R_2)^2 + (R_4)^2}\right] R_2 \\
\dfrac{d R_3}{dt} &= R_4  \\
\frac{d R_4}{dt} &= -g -b\left[\sqrt{(R_2)^2 + (R_4)^2}\right]R_4.\label{eq:Rquad}
\end{aligned}
\end{equation}
This system can be solved with a built-in initial value problem differential equations solver such as\\ {\tt SciPy.integrate.solve\_ivp()}, or MATLAB's {\tt ode45()}, with initial conditions\\ \
\begin{equation}
\vec{R}(0) = [0, v_0\cos\theta, H,v_0,\sin\theta]^T
\end{equation}  
The differential equations solver gives us a solution $\vec{R}(t) = \vec{R}(t, b, v_0, \theta, g, H)$, where the latter indicates the dependence on parameters and initial conditions. From this we can extract the $y$-component which is the third component of $\vec{R}(t)$: $y(t) = y(t,b,v_0, \theta, g, H) = R_3(t,b,v_0,\theta, g, H)$. First we use values of parameters from the scenario where the dart is dropped, and extract the value of $y(0.95, b, 0, \frac{\pi}{2}, g, 4.06)$. This allows us to use {\tt fsolve()} or {\tt fzero()} to estimate the value for $b$ for which $y(0.95, b, 0, \frac{\pi}{2}, g, 4.06)=0$.  We find $b = 0.06454$, and now we can use {\tt fsolve()} or {\tt fzero()} once again with the values of parameters from the scenario where the dart is fired straight up to find $v_0$ when $y(2.13, b, v_0, \frac{\pi}{2}, g, H) = 0$ to estimate $v_0=12.26$ m/s. 

A second method for finding parameters $b$ and $v_0$ is to solve the system of differential equations with a built-in boundary value solver such as {\tt SciPy.integrate.solve\_bvp()} or MATLAB's {\tt bvp4c()}.  These solvers must be provided with an initial guess for the solution that solve the boundary value problem and an initial guess for the unknown parameter, and will return the value of an unknown parameter along with a an array of values $t$ and the function values for $x(t)$, $x'(t)$, $y(t)$ and $y'(t)$ as the rows of $\vec{R} (t)$. 
By using the boundary value solver with the initial conditions from dropping the dart, we once again estimate $b=0.0645$, and from the initial conditions for firing the dart straight up in the air, we  once again estimate $v_0 = 12.3$ m/s.

Regardless of which method is used to estimate $b$ and $v_0$, once
estimates are obtained, we can use either MATLAB or Python's {\it
  event location} functionality with an initial value differential
equations solver like {\tt SciPy.integrate.solve\_ivp()} or MATLAB's
{\tt ode45()}.  Event location halts integration when an event, such
as the dart hitting the ground, occurs.  The initial value solver with
event location returns an array of values for $t$ and the function
values $x(t)$, $x'(t)$, $y(t)$ and $y'(t)$ are the first through
fourth rows of the solution, respectively. This allowed us to create
the plots in Figure 4. 

When $d(\theta)$ is maximized, $-d(\theta)$ is minimized.  Either MATLAB's {\tt fminsearch()} or Python's {\tt SciPy.optimize.fmin()} can be used to find the value of $\theta$ for which $-d(\theta)$ is the minimum.  Using an initial guess of $45^\circ$, we found $\theta = 40.45^\circ$ maximizes the distance traveled.  The MATLAB file {\tt nonlinair1.m} computes the angle that maximizes the distance traveled.

\begin{figure}[t]
  \begin{center}
\begin{tabular}{p{3 in} p{3 in}}
    \includegraphics[width=2.7 in]{Trajectory35.pdf} &\includegraphics[width=2.7 in]{DistAngleAll.pdf} 
\end{tabular}
\parbox{9. true cm}{\caption{Left: trajectory of dart for $\theta =35^\circ$ for all models. Right: distance vs. angle data and model results\label{fig:DistAngleAll}.}}
 \end{center}
  \end{figure}

We find the square root of the sum of squared errors for this model to be $E=1.709$, in comparison to $E=2.753$ for the model with linear air resistance and $E=5.099$ for the model with no air resistance.  This model offers substantial improvement from both of the prior models.

MATLAB code to find $b$ for physical air resistance using the initial value differential equations solver {\tt ode45()} and {\tt fzero()} is in {\tt bvalue.m}.  Code to find $b$ via the boundary value differential equations solver {\tt bvp4c()} is in {\tt bvalue2.m}.  Code to find $v_0$ also uses the boundary value differential equations solver {\tt bvp4c()} and is in {\tt v0value.m}.
Code to plot the 35$^\circ$ trajectory for all three models is in
{\tt ballistics\_trajectories.m}.  A comparison of model distances and
actual distances versus firing angle (again for all 3 models) is in
{\tt ballistics\_distances.m}, this code uses the function in {\tt
dist.m} to compute values for  the physical air resistance model. Code to find the value of $\theta$
that maximizes $d(\theta)$ for physical air resistance is in {\tt
nonlinair2.m}.



\subsection*{Regression Approach}

There are multiple ways to incorporate all of the data into the model.

In the Python code for this project, we use \eqref{eq:Rquad}
and an initial value differential equation solver with event location
to return an array of distances the dart travels for each firing
angle.  Then  Python's {\tt SciPy.Optimize.curve\_fit()} is used to find the values of
parameters $b$ and $v_0$ to minimize the least squared errors between
the model distances and the actual measured distances the dart travels
for the given firing angles.  We find $b = 0.08018$ and $v_0 = 13.27$.  This method can also be implemented in MATLAB using MATLAB's {\tt curve\_fit()} routine

In the MATLAB code {\tt distancefit2.m} (which uses {\tt distvec.m}), this same task is accomplished by minimizing the sum-of-squared errors for the model with {\tt fminsearch()}.  This can be accomplished with Python using {\tt SciPy.optimize.fmin()}.  The MATLAB code finds $b=0.0799$ and $v_0 = 13.25$ m/s.

We plot the data, the distance as predicted by the values of the parameters in the section on the Physical Model and the distances as predicted by the values of the parameters we obtained from using {\tt SciPy.Optimize.curve\_fit()}.  The two models are very close.  Results are in Figure \ref{fig:DistAngleAll2}

\begin{figure}[t]
  \begin{center}
    \includegraphics[width=4 in]{DistAngleAll2.pdf} 
 \parbox{8.1 true cm}{\caption{Distance vs. angle data and model results\label{fig:DistAngleAll2}. The two models and the data are all close.}}
 \end{center}
  \end{figure}

The square root of the sum of squared error for the model incorporating all of the data is $E = 1.551$, just slightly better than $E=1.7099$ from the physical air resistance model.  

\subsection*{Conclusions}

This project starts with a simple model for motion under gravity without air resistance, a model that can largely be analyzed with hand calculations, and proceeds to incorporate air resistance in a linear term and, after dimensional analysis, in a quadratic term that fully takes the vector properties of the model into account.  The contrast between the data and the models seen in Figure \ref{fig:DistAngleAll}, shows the importance of incorporating the correct physics (modeling assumptions) into the model.  Adding linear air resistance improves the model and the error, but this model is still far from the measured data.  It is only when we move to the physical air resistance model that the model truly tells the story of the data.  

\begin{thebibliography}{99}
\bibitem{anaconda} Anaconda, Inc. 2018. {\sl Anaconda Python Distribution}. \href{https://www.anaconda.com/download/}{https://www.anaconda.com/download/}.  Last checked 16 November 2018.
\bibitem{BoyceDiPrima} Boyce, William E. and Richard C. DiPrima. 2009. {\sl Elementary Differential Equations\/}, 9th edition. Hoboken, NJ. Wiley and Sons.
\bibitem{tracker} Brown, Douglas. 2018. {\sl Tracker: Video Analysis and Modeling Tool.\/} \href{http://physlets.org/tracker/}{http://physlets.org/tracker/}. Last checked 16 November 2018.
\bibitem{Buckingham} Buckingham, E. 1914. On Physically Similar Systems; Illustrations of the Use of Dimensional Equations. {\sl Physical Review Journal}. 4: 345--376. \href{https://journals.aps.org/pr/pdf/10.1103/PhysRev.4.345}{https://journals.aps.org/pr/pdf/10.1103/PhysRev.4.345}. Last checked 19 October 2018.
\bibitem{Edwards} Edwards, C. Henry and David E. Penney. 2008. {\sl Differential Equations and Boundary Value Problems\/} 4th Edition. Englewood Cliffs, NJ. Prentice Hall.
\bibitem{Galileo} Galilei, Galileo. 1638. {\sl Two New Sciences\/}. Leiden. Reprinted in Hawking, Stephen. 2004. {\sl On the Shoulders of Giants: the Great Works of Physics and Astronomy}. Philadelphia, PA. Running Press.
\bibitem{Giordano} Giordano, Nicholas J. and Hisao Nakanishi. 2005. {\sl Computational Physics\/} 2nd edition.  Boston, MA. Pearson. 
\bibitem{Malthe} Malthe-Sorenssen, Anders. 2015. {\sl Elementary Mechanics Using Python\/}. Heidelberg Germany: Springer.
\bibitem{Nagle}  Nagle, R. Kent, Edward B. Saff, and Arthur David Snider. 2000. {\sl Fundamentals of Differential Equations and Boundary Value Problems\/} 4th edition. Boston, MA. Pearson.
\bibitem{Noonburg} Noonburg, Virginia. 2014. {\sl Ordinary Differential Equations: From Calculus to Dynamical Systems}. Washington, DC. MAA Press.
\end{thebibliography}
}





\end{document}
