#!/usr/bin/env python
# coding: utf-8

# <style>
# @import url(https://www.numfys.net/static/css/nbstyle.css);
# </style>
# <a href="https://www.numfys.net"><img class="logo" /></a>
# 
# # Shooting a projectile into the air
# 
# ## Examples - Mechanics
# <section class="post-meta">
# By Niels Henrik Aase, Eilif Sommer Øyre, Jaakko Akola, and Jon Andreas Støvneng
# </section>
# Last edited: February 20th 2019
# 
# ___
# This notebook is based on a project assignment [[1]](#rsc) given in the course TFY4245 Classical Mechanics at NTNU, fall 2018.
# ___

# ## Introduction
# The trajectory of a projectile in a uniform gravity field is a classic and dear physics problem. However, in most introductory physics courses, one works under the assumption that the frictional forces can be neglected. This is of course not always the case, especially not at high speeds. If frictional forces is included, the path of the projectile can no longer be described analytically, and the trajectory must be calculated numerically. To do this, one must solve a coupled set of differential equations, which provides us with an excellent opportunity to implement an ODE-solver on vector form. This turns out to be an elegant and easy implementation, that can be used for solving both differential equations of higher order and coupled differential equations. 

# ## Theory
# ### Equations of motion
# The governing equations for a projectile is derived from Newtons second law, $\vec{F} = m\vec{a}$, where $\vec{F}$ is the force acting on the projectile, $m$ is the mass, and $\vec{a}$ is the acceleration of the ball. In this example there are two forces acting on the ball. The gravitational force 
# 
# $$\vec{F_G} = m\vec{g},$$
# 
# where $\vec{g}$ denotes the the standard gravity. The other force is the drag force
# 
# $$\vec{F_D} = -\frac{1}{2}\rho ACv_r^2\hat{v_r}.$$ 
# 
# $A$ is the cross-sectional area, $C$ is the drag coefficient, and $v_r$ is the velocity of the projectile relative to the fluid it is traveling through, $\rho$ is the air density and $\hat{v_r}$ is the direction of the projectile with length unity. A more practical way of writing this equation is
# 
# $$\vec{F_D} = - \frac{B}{\rho_0}\rho v_r\left(\vec{v}-\vec{V}\right),$$
# 
# where $\rho_0$ is the air density at sea level altitude, $B$ is some constant, $\vec{v}$ is the velocity of the projectile, and $\vec{V}$ is the velocity of the wind. 

# In this example we will use the shells from the German cannon "the Paris gun" as the projectiles, and the numerical constants mentioned above will be based on historical data of the cannon. The Paris gun operated during World War I and bombarded Paris from far away, the distances often exceeded 100 kilometers, with the velocities of the shells being correspondingly high.
# 
# The complete equation for the drag force includes an extra term, which is not mentioned above, namely the term originating from Stokes' law. For a spherical projectile the term takes the form 
# 
# $$\vec{F_D'} = -6 \pi \eta \vec{v},$$
# 
# where $\vec{F_D'}$ denotes the extra term, and $\eta$ denotes the air's viscosity. At high velocities this term is dwarfed by $\vec{F_D}$ as $\vec{F_D}$ is proportional to $v^2$, thus making the drag force $\vec{F_D}$ an excellent approximation for the total drag force.  

# The last part of the theory needed before we can find the equations of motion is an expression for the air density. As the shells of the Paris gun will reach the stratosphere, it is not reasonable to treat $\rho$ as a constant. In this notebook we will use an adiabatic air density model. The air density is then given by 
# 
# $$\rho = \rho_0 \left(1 - \frac{ay}{T_0}\right)^\alpha,$$
# 
# where $\rho$ and $T_0$ denotes the air density and temperature at sea level, $y$ denotes the altitude, $a$ and $\alpha$ are constants $a= 6.5 \cdot 10^{-3}$ K/m, $\alpha$ = 2.5. At high altitudes the expression in the parentheses will become negative, thus making the air density an imaginary number. To avoid this, the air density is set equal to zero when reaching these altitudes. This critical height is around 44 kilometers, well into the stratosphere. 

# In a two-dimensional coordinate system, where $x$ and $y$ are the displacement from the origin in horizontal and vertical direction respectively, the equations of motion become
# $$
# \begin{align}
# m\ddot{x} &= -\frac{B}{\rho_0}\rho\left(\dot{x}-V_x\right)\left|\vec{v}-\vec{V}\right|,\\
# m\ddot{y} &= -\frac{B}{\rho_0}\rho\left(\dot{y}-V_y\right)\left|\vec{v}-\vec{V}\right|-mg,
# \end{align}
# $$
# 

# where $V_x$ and $V_y$ denotes the velocity of the wind in the horizontal and vertical direction. Simplifying even further, and introducing $u$ and $v$ as the velocities in the horizontal and vertical direction yields four coupled differential equations
# 
# \begin{equation}
# \begin{aligned}
# \dot{x} &= u, \\
# \dot{u} &= -\frac{B}{m}\Big[ 1 - \frac{ay}{T_0} \Big]^\alpha \left(u-V_x\right)\Big[ (u-V_x)^2+(v-V_y)^2 \Big]^{1/2},\\
# \dot{y} &= w, \\
# \dot{w} &= -\frac{B}{m}\Big[ 1 - \frac{ay}{T_0} \Big]^\alpha \left(v-V_y\right)\Big[ (u-V_x)^2+(v-V_y)^2 \Big]^{1/2}-g.\\
# \end{aligned}
# \label{ODEs}
# \end{equation}

# ### Discretization: Runge-Kutta 4th order
# 
# First we define a set of $N$ discrete time values $$
# t_n=t_0+n\cdot h~~~~~\mathrm{with}~~~~n=0,1,2,3,...,N.
# $$
# 
# $h$ is the time step length and is calculated by the formula $h = \frac{t_N-t_0}{N}$, where $t_0$ and $t_N$ denotes the initial and final time values of the projectile's flight. These $N$ time values correspond to $N$ number of $x$, $y$, $u$, and $w$ values, each denoted $x_n$, $y_n$, $u_n$ and $w_n$. A more practical way of denoting this collection of values is on vector form
# $$
# \vec{w_n} = \begin{bmatrix}x_n \\ y_n \\ u_n \\ v_n\end{bmatrix}.
# $$
# 
# This notation makes it possible to implement the Runge-Kutta method in an elegant way. It is always possible to reduce an ordinary differential equation of order $q$ to a problem of solving a set of two $q-1$ order coupled differential equations. We have a coupled set of two second order differential equations, thus we need to solve a set of four first order coupled differential equations. By introducing a function $f$ that transforms $\vec{w_n}$ in the way described by equation \eqref{ODEs}, one can implement the Runge-Kutta method directly. In our case, $f$ is given by
# 
# $$
# f(\vec{w_n} ) = f \begin{bmatrix}x_n \\ y_n \\ u_n \\ v_n\end{bmatrix} = \begin{bmatrix}u_n \\ v_n \\ -\frac{B}{m}\Big[ 1 - \frac{ay_n}{T_0} \Big]^\alpha \left(u_n-V_x\right)\Big[ (u_n-V_x)^2+(v_n-V_y)^2 \Big]^{1/2} \\ -\frac{B}{m}\Big[ 1 - \frac{ay_n}{T_0} \Big]^\alpha \left(v_n-V_y\right)\Big[ (u_n-V_x)^2+(v_n-V_y)^2 \Big]^{1/2}-mg\end{bmatrix},
# $$
# with no time dependence, such that 
# $$
# \dot{\vec{w_n}}=f(\vec{w_n})
# $$
# 
# describes the complete system of differential equations. For further reading, check out the [notebook on the Runge-Kutta method](https://nbviewer.jupyter.org/urls/www.numfys.net/media/notebooks/runge_kutta_method.ipynb). The following code implements this, and can easily be expanded to solve either an even higher order ODE or an ODE describing a three-dimensional problem. 
# 

# ## Implementation

# In[2]:


# Importing necessary packages
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
get_ipython().run_line_magic('matplotlib', 'inline')

newparams = {'figure.figsize': (15, 7), 'axes.grid': False,
             'lines.markersize': 10, 'lines.linewidth': 2,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral', 'figure.dpi': 200}
plt.rcParams.update(newparams)

# Constants

g = 9.81             # Force of gravity per kilo on earth's surface
alpha = 2.5          # Parameter in the adiabatic air density model
a = 6.5 * 10 ** (-3) # Parameter in the adiabatic air density model
T_0 = 288            # Temperature at sea level
m = 50               # Mass of the projectile
B = 2 * 10 ** (-3)   # Constant based on the Paris gun
v_0 = 1640           # Initial velocity


# Wind strength. This is a constant that never varies 
V = np.zeros(2) # Starting of with no wind


def f(w):
    """ A function describing the right hand side of the equations
    of motion/the set of ODEs.
    Parameter:
        w       vector containg the needed coordinates and their velocities 
    """
    # V : the vector describing the wind strength and direction
    temp_vector = np.zeros(4)
    temp_vector[0] = w[2]
    temp_vector[1] = w[3]
    # Saving the next parameters in order to optimize run time
    k = (1 - a * w[1] / T_0)   # Air density factor
    s = np.sqrt((w[2]-V[0]) ** 2 + (w[3]-V[1]) ** 2 )
    if k>0:
        temp_vector[2] = - B / m * k ** (alpha) * (w[2]-V[0]) * s
        temp_vector[3] = - B / m * k ** (alpha) * (w[3]-V[1]) * s - g
    else:
        temp_vector[2] = 0 
        temp_vector[3] = - g
    return temp_vector 


def RK4_step(f, w, h):
    """Performing a single step of the Runge-Kutta fourth order method.
    Parameters:
        f      RHS of ODEs to be integrated
        w      numerical approximation of w at time t
        h      unit/integration step length
    Returns:
        numerical approximation of w at time t+h
    """
    s1 = f(w)
    s2 = f(w + (h / 2) * s1)
    s3 = f(w + (h / 2) * s2)
    s4 = f(w + h * s3)
    return w + (h / 6) * (s1 + (2 * s2) + (2 * s3) + s4)


# Now we can finally start to explore the trajectories of the Paris gun. The functions below will be helpful going forward.

# In[3]:


def shoot(theta, v0):
    """ Initializes the vector w (x and y position and velocity of the
    projectile) given a initial shooting angle, theta, and 
    absolute velocity, v0. 
    """
    w = np.zeros(4)
    w[2] = v0 * np.cos(np.deg2rad(theta))
    w[3] = v0 * np.sin(np.deg2rad(theta))
    return w


def projectile_motion(h, theta):
    """ Calculates the motion of the projectile using the functions
    defined above. While the projectile is in the air (w[1] >=0) the 
    position and velocity is updated using a single step of RK4.
    Parameters:
        h        unit/integration step length
        theta    initial shooting angle
    Returns:
        X_list   array of the projectile's x-position
        Y_list   array of the porjectile's y-position
    """
    w = shoot(theta, v_0)
    X_list = np.zeros(0)
    Y_list = np.zeros(0)
    while  w[1] >= 0:
        w = RK4_step(f, w, h)
        X_list = np.append(X_list, w[0])
        Y_list = np.append(Y_list, w[1])
    return X_list, Y_list


# ### Finding the optimal launch angle
# With no drag force, it is trivial to show that the angle maximizing the range of a projectile is $\theta = 45^\circ$, where $\theta$ is the angle between the horizontal axis and the direction of the initial velocity. It is interesting to se how the air drag affects this maximum angle, both with and without wind. The following code finds the optimal angle to one decimal precision, and shows how the range depends on the firing angle. It is also instructive to see some of the trajectories of the projectile, so some of the trajectories are included. 

# In[4]:


def find_optimal_angle(h):
    """ Given an integration time step, this function calculates the optimal initial 
    shooting angle for the projectile to obtain maximum range, in x-direction. The 
    set of angles tested, with their corresponding range, along with the optimal
    angle are returned.
    """
    
    record = 0        # Placeholder variable that holds the maximum range
    optimal_angle = 0 # Placeholder variable that holds the angle yielding the maximum range
    # Lists containing the initial angle and its corresponding range
    theta_list = np.zeros(0) 
    range_list = np.zeros(0)
    for theta in range (1,90,2):
        x_list, y_list = projectile_motion(h, theta)
        # Using linear interpolation do determine the landing point more precisely
        m = (y_list[-1] - y_list[-2]) / (x_list[-1] - x_list[-2])   # The landing point
        x_range = - y_list[-1] / m + x_list[-1] 
        theta_list = np.append(theta_list, theta)
        range_list = np.append(range_list, x_range)
        # Update records
        if x_range >= record:
            record = x_range
            optimal_angle = theta

    # Rerunning the same code on a smaller interval in order to approximate the optimal angle
    # more precicely
    theta_list_smaller = np.linspace(optimal_angle - 2, optimal_angle + 2, 41)
    for theta_small in theta_list_smaller:
        x_list, y_list = projectile_motion(h, theta)
        # Again, using linear interpolation do determine the landing point more precisely
        m = (y_list[-1] - y_list[-2]) / (x_list[-1] - x_list[-2])
        x_range = - y_list[-1] / m + x_list[-1]
        if x_range >= record:
            record = x_range
            optimal_angle = theta_small
            
    return theta_list, range_list, optimal_angle

theta, x , best = find_optimal_angle(0.1)

print("The optimal angle is: ", best, " degrees")


# In[5]:


plt.plot(theta, x/1000)
plt.title(r"Projectile range as a function of shooting angle, $\theta$")
plt.xlabel(r"$\theta $ [$\degree$]")
plt.ylabel(r"range [km]")

def trajectories(h):
    plt.figure()
    plt.title("Projectile trajectories by alternating shooting angle")
    plt.xlabel(r"x [km]")
    plt.ylabel(r"y [km]")
    theta_list = np.arange(30.0,75,5)
    for angle in theta_list:
        x_list, y_list = projectile_motion(h, angle)
        plt.plot(x_list/1000, y_list/1000, '--', label=r'$\theta = $%.i $\degree$'%(angle))
    plt.legend(loc='best')
    plt.gca().set_ylim(bottom=0)
    plt.show()
trajectories(0.1)


# In contrast to projectile motion with no friction, maximum range is not obtained with firing angle $\theta = 45^\circ$. Because  the drag force is lower at higher altitudes, the range increases when the projectile spends the majority of it's flight at higher altitudes. This explains why the maximum range is obtained at $\theta = 55^\circ$.

# ### Numerical validity
# 
# In order to determine the validity of the numerical implementation, we can use conservation of energy. During the projectile's trajectory the drag force $\vec{F_D}$ will do work on the projectile. The cumulative work done on the projectile is given by 
# 
# $$
# W = \int \vec{F_D} \cdot \vec{dr}.
# $$
# 
# In our discrete model, we can easily calculate this integral by summing over all the infinitesimal contributions, and using the fact that $\vec{F_D}$ always is antiparallel with $\vec{dr}$ such that the cumulative work at time $t_M$, $M<N$, is equal to
# 
# $$
# W = \sum_{n=1}^{M} |\vec{F_{D_{n-1}}}| \cdot |\vec{r_n} - \vec{r_{n-1}}|,
# $$
# where $\vec{F_{D_{n-1}}}$ denotes the drag force at position $\vec{r_{n-1}}$. However, this can again be written to the form
# 
# $$
# W = h \sum_{n=1}^{M} |\vec{F_{D_{n-1}}}| \cdot |\vec{v_{n-1}}|,
# $$
# which is the most convenient for our purposes. It must be noted that the choice of  index that is used for the force and speed respectively is non-trivial, as they both will vary during the interval $dr$. Here both the speed of the projectile and the force at point $\vec{r_{n-1}}$ is used. This is mainly a pragmatic choice, as this combination yields the lowest, and thus the most realistic energy loss. Thus we can finally express the energy of the projectile (which of course is equal to it's kinetic and potential energy) and the work done on the air at time $t_M$
# 
# $$
# E(t_M) = h \sum_{n=1}^{M} |\vec{F_{D_{n-1}}}| \cdot |\vec{v_{n-1}}| + \frac{1}{2} m v_m^2 + mgy_m,
# $$
# 
# where $v_m$ is the absolute velcoity of the projectile, and $y_m$ is the height of the projectile over the ground (zero potential reference point).
# 
# Now all that remains is to check if this expression for the energy indeed remains constant while it is in the air. 
# 

# In this notebook, the Runge-Kutta method is implemented without much difficulties. However, it is interesting to see whether or not it was necessary to use a 4th order differential equation solver, instead of using the most primitive ODE solver, Euler's method. The following code will use both Runge-Kutta and Euler's method, and compare the energy loss during the motion for both of them.

# In[6]:


def euler_step(f ,w ,h):
    """Simple implementation of Euler's method on vector form
    """
    return w + h * f(w)

def energy_conservation(h, theta, func):
    """ This function performs the same steps as "projectile_motion()" and "f()"
    but extracts the value for the drag, F_D, at every step. This enables the 
    calculation of work and storage of energy at each time step.
    Parameters:
        h           integration time step
        theta       initial shooting angle
        func        a function defining the integration step (euler or RK4)
    Returns:
        x_list      array of the trajectory x-coordinates
        y_list      array of the trajectory y-coordinates
        time_steps  array of the elapsed time at each step
        energy      array of the energy at each step
        work        the total work by the drag force
    """
    w = shoot(theta, v_0)
    
    x_list = np.zeros(1)
    y_list = np.zeros(1)
    work = 0
    
    # Initial total energy
    v = np.sqrt(w[2] ** 2 + w[3] ** 2)
    energy = np.array([0.5 * m * v ** 2])
    
    
    while w[1] >= 0:
        w = func(f, w, h)
        
        # Updating lists
        x_list = np.append(x_list, w[0])
        y_list = np.append(y_list, w[1])
        
        # Temporary parameters for checking if high altitude
        k = (1 - a * y_list[-2] / T_0)
        s = np.sqrt((w[2] - V[0]) ** 2 + (w[3] - V[1]) ** 2)
        if k > 0:
            F_D = B * k ** (alpha) * s ** 2
        else:
            F_D = 0
        
        # Add the amount of work done by drag this time interval
        work +=  v * F_D * h
        # Updating the speed in order to calculate the total energy
        v = np.sqrt(w[2] ** 2 + w[3] ** 2)
        # Use np.sum on the work_list to get the cumulative work
        energy = np.append(energy, 0.5 * m * v ** 2 + m * g * w[1] + work)
    
    n = len(x_list)
    time_steps = np.linspace(0, n * h, n)
    
    return x_list, y_list, time_steps, energy, work
 

# Setting initial values
theta = 45
h = 0.1
# Getting the results
x_list_euler, y_list_euler, time_steps_euler, energy_euler, work = energy_conservation(h, theta, euler_step)
x_list_RK4, y_list_RK4, time_steps_RK4, energy_RK4, work = energy_conservation(h, theta, RK4_step)

print("At firing angle %.1f degrees and with time step h = %.2f seconds:" %(theta, h))
print("Energy lost with RK4: %.3E = %.2f percent of total energy" %(energy_RK4[0] - energy_RK4[-1],
                                          (energy_RK4[0] - energy_RK4[-1]) / energy_RK4[0] * 100))
print("Energy lost with Euler's method: %.3E = %.2f percent of total energy" %(energy_euler[0] - energy_euler[-1],
                                        (energy_euler[0] - energy_euler[-1]) / energy_euler[0] * 100))

# Plotting the separate trajectories 
plt.plot(x_list_euler/1000, y_list_euler/1000, label ="Euler")
plt.plot(x_list_RK4/1000, y_list_RK4/1000, label ="RK4")
plt.title("Firing angle = " + str(theta) + "$\degree$")
plt.xlabel(r"$x$ [km]")
plt.ylabel(r"$y$ [km]")
plt.gca().set_ylim(bottom=0)
plt.legend()
plt.show()


# From the figure above, it is clear that Euler's method work surprisingly well, even though it us just ODE-solver of first order. It produces nearly the same trajectory as Runge Kutta 4th method. The most likely explanation of the small difference between Runge Kutta 4th order and Euler's method is that the coupled set of ODEs we are solving are quite linear, in other words, the derivatives vary quite slowly. Thus Euler's method produce satisfactory results, even though it is only of first order.
# 
# It might seem like loosing around $10^3$ joules is a lot, but one needs to remember the enormous dimensions of the Paris gun, and that this is a small fraction of the original energy of the projectile.
# 

# ## The physics behind the numerics
# Even though the amount of disappeared energy is relatively small, it is still interesting to study the nature of the loss. In our case it turns out that it has an explanation that has roots in both the physics and the numerics of the problem. 

# In[7]:


# When does the energy disappear during the flight?
starting_energy_RK4 = np.repeat(energy_RK4[0], len(energy_RK4))
starting_energy_euler = np.repeat(energy_euler[0], len(energy_euler))

plt.plot(time_steps_RK4, starting_energy_RK4 - energy_RK4, label = "RK4")
plt.plot(time_steps_euler, starting_energy_euler - energy_euler, label = "Euler's method")
plt.legend()
plt.xlabel(r"$t$ [s]")
plt.ylabel(r"$\Delta $ E [J]")

plt.show()


# Exploring the numerical validity as a function of the firing angle
energy_loss_euler = np.zeros(0)
energy_loss_RK4 = np.zeros(0)
angles = np.zeros(0)
for i in range(90):
    energy_euler = energy_conservation(0.1, i, euler_step)[3] # only returns the fourth return value
    energy_loss_euler = np.append(energy_loss_euler, energy_euler[0] - energy_euler[-1])
    energy_RK4 = energy_conservation(0.1, i, RK4_step)[3]
    energy_loss_RK4 = np.append(energy_loss_RK4, energy_RK4[0] - energy_RK4[-1])
    angles = np.append(angles, i)

plt.plot(angles, energy_loss_RK4, label ="RK4")
plt.plot(angles, energy_loss_euler, label ="Euler's method")
plt.legend()
plt.xlabel(r"$\theta $ [$\degree$]")
plt.ylabel(r"$\Delta $ E [J]")
plt.show()


# The first figure illustrates how much the energy differs from the original energy of the projectile during the the trajectory. Both Euler's method and Runge-Kutta demonstrate the same numerical behavior, the main contribution to the error in energy occurs at lower altitudes, at the start and end of the trajectory. As the projectile quickly reaches high altitudes where the air density is lower, the drag force will also decrease. When the drag force is negligible the motion of the projectile can be described analytically by a second order expression, an easy task, especially for Runge-Kutta as it is a fourth order ODE solver. The large energy loss at the start of the trajectory with Euler's method is expected, the drag force is proportional to $v^3$ and because the speed is much higher during the start of the trajectory, the derivatives of $u$ and $v$ will vary quicker and Euler's method will in a sense fail to keep up. 
# 
# The same explanations are valid for the second figure. The less time the projectile spends in lower altitudes, the smaller the energy loss. If the firing angle is low, the projectile will not reach the altitudes where the air density and energy loss is lower. 

# In[8]:


# Exploring the numerical validity as a function of the step length h
theta = 45

h_list = np.linspace(0.01, 0.1, 30)

energy_loss_euler = np.zeros(0)
energy_loss_RK4 = np.zeros(0)

for h in h_list:
    energy_euler = energy_conservation(h, theta, euler_step)[3]
    energy_loss_euler = np.append(energy_loss_euler, energy_euler[0] - energy_euler[-1])
    energy_RK4 = energy_conservation(h, theta, RK4_step)[3]
    energy_loss_RK4 = np.append(energy_loss_RK4, energy_RK4[0] - energy_RK4[-1])
plt.plot(h_list, energy_loss_RK4, label ="RK4")
plt.plot(h_list, energy_loss_euler, label ="Euler's method")
plt.legend()
plt.xlabel(r" h [s]")
plt.ylabel(r"$\Delta $ E [J]")
plt.show()


# The last numerical check is to see how the energy deviance depends on the time step length $h$. Unsurprisingly the deviance shrinks as $h$ becomes smaller, and the error seems to be linear, both for Euler's method and Runge-Kutta fourth order.
# 
# As for which ODE solver that is the superior for this project, one can certainly argue that the more realistic results produced by Runge-Kutta fourth order makes it the preferable candidate. Especially considering that the Runge-Kutta implementation only requires four extra lines of code, it is hard to se any major drawbacks of using Runge-Kutta instead of Euler's method. However it must be noted that they produce nearly the same trajectory, so as a simple demonstration of the trajectory, Euler's method will work just fine.

# ### Exploring the effects of the wind - Small Bertha
# 
# In this section we will see how the wind affects the motion of the projectile. As the dimensions of the Paris gun is far to big to be affected by any realistic wind strength, dimensions from a fictitious cannon "Small Bertha" is used. The name "Small Bertha" originates from the series of guns used by Germany in World War 1, guns of this series were commonly referred to as "Big Bertha". Because the air density will be approximately constant, we now expect that the maximum range is obtained at an angle lower than previously. We will study how the range and trajectory of the projectile is affected.

# In[9]:


# Re-initializing variables in order to scale down the dimensions from the Paris gun
m = 10             # Mass of the projectile
B = 5 * 10 ** (-3) # Constant based on the Paris gun
v_0 = 164          # Initial velocity


# Wind strength
V = np.zeros(2)

print("The optimal angle is: ", find_optimal_angle(0.1)[2], 
      " degrees with no wind. \nThe trajectories with no wind:")            
trajectories(0.1)

V[0] = 20 # Strong tailwind in the horizontal direction
print("The optimal angle is: ", find_optimal_angle(0.1)[2], 
      " degrees with strong tailwind.", 
      "\nThe trajectories with strong tailwind:")            
trajectories(0.1)

V[0] = -20 # Strong headwind in the horizontal direction
print("The optimal angle is: ", find_optimal_angle(0.1)[2], 
      " degrees with strong headwind", 
      "\nThe trajectories with strong headwind:")                     
trajectories(0.1)


# The wind affects the trajectories as expected. When the projectile experiences strong tailwind, the drag force decreases, making the motion more similar to the analytical trajectory of a projectile without friction. Thus the optimal firing angle will be closer to $45^\circ$, the optimal angle for projectile motion without friction. The opposite case is true when the projectile experiences strong headwind. The drag force increases, which makes the optimal firing angle lower, as the projectile is best served spending as little time as possible in the air. 

# ___
# <a id="rsc"></a>
# ## Resources and Further Reading
# <a>[1]</a>: Jaakko Akola, *Computational assignment 1*, NTNU, 2018. <br />
# <a>[2]</a>: Jaakko Akola, *Appendix for computational assignment 1*, NTNU, 2018. <br />