Skip to content
Snippets Groups Projects
Commit 779f30fb authored by cprutean's avatar cprutean
Browse files

Update file Checkpoint5.ipynb

parent 1ebaf281
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Checkpoints/Checkpoint5.ipynb
+ 5
− 5
View file @ 779f30fb
%% Cell type:markdown id:abaa6e4a tags: %% Cell type:markdown id:abaa6e4a tags:
## Checkpoint 5 ## Checkpoint 5
### Aim ### Aim
To write two Python programs to numerically analyse the path of a projectile under gravity and drag force. There are two tasks, one to display the trajectory and the second to display the relation between final kinetic energy and launch angle. To write two Python programs to numerically analyse the path of a projectile under gravity and drag force. There are two tasks, one to display the trajectory and the second to display the relation between final kinetic energy and launch angle.
### Numerical Integration ### Numerical Integration
In this checkpoint we will use the simplest numerical integration technique called First Order Euler. This involves iteratively calculating the position and velocity at a set of times $t_i$=iΔt for i=0,1,2,… where we assume that the position and velocity at time $t_{i+1}$ depend only on the position and velocity at time $t_i$. In particular if at time $t_i$ the particle has In this checkpoint we will use the simplest numerical integration technique called First Order Euler. This involves iteratively calculating the position and velocity at a set of times $t_i$=iΔt for i=0,1,2,… where we assume that the position and velocity at time $t_{i+1}$ depend only on the position and velocity at time $t_i$. In particular if at time $t_i$ the particle has
$$\text{Position} = \vec{x}_i$$ $$\text{Position} = \vec{x}_i$$
$$\text{Velocity} = \vec{v}_i$$ $$\text{Velocity} = \vec{v}_i$$
$$\text{Acceleration} = \vec{a}_i$$ $$\text{Acceleration} = \vec{a}_i$$
then at time $t_{i+1}=t_i+Δt$ the position and velocity of the particle can be approximated by then at time $t_{i+1}=t_i+Δt$ the position and velocity of the particle can be approximated by
$$\vec{x}_{i+1} = \vec{x}_i+Δt\vec{v}_i$$ $$\vec{x}_{i+1} = \vec{x}_i+Δt\vec{v}_i$$
$$\vec{v}_{i+1} = \vec{v}_i+Δt\vec{a}_i$$ $$\vec{v}_{i+1} = \vec{v}_i+Δt\vec{a}_i$$
So if we start at time $t_0=0$ at position $\vec{x}_0$ and velocity $\vec{v}_0$, then provided we can calculate the force on the particle, and hence the acceleration $\vec{a}_i$, and use small enough time steps Δt, we can iteratively trace out the path by repeated use of the above relations. So if we start at time $t_0=0$ at position $\vec{x}_0$ and velocity $\vec{v}_0$, then provided we can calculate the force on the particle, and hence the acceleration $\vec{a}_i$, and use small enough time steps Δt, we can iteratively trace out the path by repeated use of the above relations.
### Forces on Particle ### Forces on Particle
![ch5.png](ch5.png) ![ch5.png](ch5.png)
If we consider a particle of mass m travelling with velocity $\vec{v}$ subject to gravity and drag from a fluid as shown below, then the drag force is given by If we consider a particle of mass m travelling with velocity $\vec{v}$ subject to gravity and drag from a fluid as shown below, then the drag force is given by
$$ \vec{F}_d=−\frac{1}{2}\rho AC_Dv^2\hat{v}$$ $$ \vec{F}_d=−\frac{1}{2}\rho AC_Dv^2\hat{v}$$
where $\rho$ is the density of the fluid, A is the cross-sectional area of the particle, $C_D$ is the drag coefficient and $\hat{v}=\vec{v}/|v|$, being the unit vector in the direction of the velocity. where $\rho$ is the density of the fluid, A is the cross-sectional area of the particle, $C_D$ is the drag coefficient and $\hat{v}=\vec{v}/|v|$, being the unit vector in the direction of the velocity.
Thus the acceleration of of the particle is given by $\vec{a}=\frac{1}{m}\vec{F}_d−g\hat{j}$ where g is acceleration due to gravity, being 9.81 $ms^{−2}$. We can write this as: Thus the acceleration of of the particle is given by $\vec{a}=\frac{1}{m}\vec{F}_d−g\hat{j}$ where g is acceleration due to gravity, being 9.81 $ms^{−2}$. We can write this as:
$$ \vec{a} =−\beta v^2\hat{v}−g\hat{j}$$ $$ \vec{a} =−\beta v^2\hat{v}−g\hat{j}$$
where $\beta=\frac{\rho AC_D}{2m}$, where $\beta$ is the normalised drag coefficient for the system. where $\beta=\frac{\rho AC_D}{2m}$, where $\beta$ is the normalised drag coefficient for the system.
In this problem the position, velocity and acceleration are all two dimensional vectors so they have a horizontal and vertical component that you have to treat separately in the calculation. In this problem the position, velocity and acceleration are all two dimensional vectors so they have a horizontal and vertical component that you have to treat separately in the calculation.
So in terms of components in x (horizontal) and y (vertical) we can write: So in terms of components in x (horizontal) and y (vertical) we can write:
$\vec{x}=(xx,xy); \vec{v}=(vx,vy); \vec{a}=(ax,ay)$ $\vec{x}=(xx,xy); \vec{v}=(vx,vy); \vec{a}=(ax,ay)$
where the components of acceleration from the above equations are: where the components of acceleration from the above equations are:
$a_x=−\beta vv_x$ and $a_y=−\beta vv_y−g$ $a_x=−\beta vv_x$ and $a_y=−\beta vv_y−g$
where v is the magnitude of the velocity, so being $v=\sqrt{v^2_x+v^2_y}$ where v is the magnitude of the velocity, so being $v=\sqrt{v^2_x+v^2_y}$
### Task ### Task
**Task 1:** Write a Python program to plot out the trajectory of a projectile starting at position $\vec{x}_0=(0,0)$ for specified initial velocity and normalised drag coefficient. You program should prompt for **Task 1:** Write a Python program to plot out the trajectory of a projectile starting at position $\vec{x}_0=(0,0)$ for specified initial velocity and normalised drag coefficient. You program should prompt for
- Magnitude of the initial velocity $v_0$ in $ms^{−1}$. - Magnitude of the initial velocity $v_0$ in $ms^{−1}$.
- Angle θ from the horizontal of initial velocity in degrees. - Angle θ from the horizontal of initial velocity in degrees.
- β the normalised drag coefficient. - β the normalised drag coefficient.
- Δt the step interval in seconds. - Δt the step interval in seconds.
The path should be traced and displayed using Matplotlib from the inital starting position until the particle re-crosses the x-axis. Your graph should have a suitable title and axis labels. Your program should also print the range of the projectile to the terminal. The path should be traced and displayed using Matplotlib from the inital starting position until the particle re-crosses the x-axis. Your graph should have a suitable title and axis labels. Your program should also print the range of the projectile to the terminal.
Test your program with: Test your program with:
- Case 1: $v_0 = 10 ms^{−1}$, $\theta=45^{\circ}$, $\beta$=0.0, $\Delta t$=0.001 - Case 1: $v_0 = 10 ms^{−1}$, $\theta=45^{\circ}$, $\beta$=0.02, $\Delta t$=0.001
- Case 2: $v_0 = 10 ms^{−1}$, $\theta=60^{\circ}$, $\beta$=0.3, $\Delta t$=0.001 - Case 2: $v_0 = 10 ms^{−1}$, $\theta=60^{\circ}$, $\beta$=0.3, $\Delta t$=0.001
**Task 2:** Modify the above program to produce a plot of the ratio of final kinetic energy to initial kinetic energy against lauch angle θ in the range $0 \to 90^{\circ}$. Your program should prompt for : **Task 2:** Modify the above program to produce a plot of the ratio of final kinetic energy to initial kinetic energy against lauch angle θ in the range $0 \to 90^{\circ}$. Your program should prompt for :
- Initial velocity $v_0$. - Initial velocity $v_0$.
- Normalised drag coefficient $\beta$. - Normalised drag coefficient $\beta$.
- Timestep $\Delta t$. - Timestep $\Delta t$.
and plot out, using Matplotlib, a graph of $K_f / K_i$ against $\theta$. and plot out, using Matplotlib, a graph of $K_f / K_i$ against $\theta$.
### Checkpoint ### Checkpoint
- To help you with this checkpoint we've already provided you with a program skeletong structure below. - To help you with this checkpoint we've already provided you with a program skeleton structure below.
- Test your program with the following conditions: - Test your program with the following conditions:
- $v_0 = 10 ms^{−1}$, $\beta$=0.02, $\Delta t$=0.001 - Case 1: $v_0 = 10 ms^{−1}$, $\beta$=0.02, $\Delta t$=0.001
- $v_0 = 10 ms^{−1}$, $\beta$=0.3, $\Delta t$=0.001 - Case 2: $v_0 = 10 ms^{−1}$, $\beta$=0.3, $\Delta t$=0.001
- Call a demonstrator, show them your code, and run your programs for the two different conditions. - Call a demonstrator, show them your code, and run your programs for the two different conditions.
### Assessment ### Assessment
To pass this checkpoint: To pass this checkpoint:
- At least one of the progammes works fully with the supplied test examples. - At least one of the progammes works fully with the supplied test examples.
- A good attempt has been made for the second program. - A good attempt has been made for the second program.
- The output graphs have the correct title and labels and the code is well structured with sensible use of functions called from the main() function. - The output graphs have the correct title and labels and the code is well structured with sensible use of functions called from the main() function.
- There are sensible header comments and sufficient in-line comments to explain the operation. - There are sensible header comments and sufficient in-line comments to explain the operation.
### Next steps ### Next steps
Once you have completed this checkpoint, please continue with the [Week5](../WeeklyTasks/Week5.ipynb) tasklist. Once you have completed this checkpoint, please continue with the [Week5](../WeeklyTasks/Week5.ipynb) tasklist.
%% Cell type:code id:2539ac4d tags: %% Cell type:code id:2539ac4d tags:
``` python ``` python
### Task 1 - plotting the trajectory ### Task 1 - plotting the trajectory
### Define here your imports ### Define here your imports
g = 9.81 #m/s**2 g = 9.81 #m/s**2
def set_initial(v_initial, theta): def set_initial(v_initial, theta):
""" Set inital conditions """ Set inital conditions
""" """
def acceleration(vx, vy, beta): def acceleration(vx, vy, beta):
""" Calculate the acceleration """ Calculate the acceleration
""" """
def step_forward(x, y, vx, vy, beta, delta_t): def step_forward(x, y, vx, vy, beta, delta_t):
""" Do a forward step """ Do a forward step
""" """
def Checkpoint5_1(v_initial, theta, beta, delta_t): def Checkpoint5_1(v_initial, theta, beta, delta_t):
""" Main program to read in from terminal, do iteration, and plot out graph """ Main program to read in from terminal, do iteration, and plot out graph
""" """
print("Trajectory for (v0, theta, beta, delta_t) = ({0}, {1}, {2}, {3})".format(v_initial, theta, beta, delta_t)) print("Trajectory for (v0, theta, beta, delta_t) = ({0}, {1}, {2}, {3})".format(v_initial, theta, beta, delta_t))
... ...
t = 0.0 t = 0.0
x, y, vx, vy = set_initial(v_initial, theta) x, y, vx, vy = set_initial(v_initial, theta)
... ...
x, y, vx, vy = step_forward(x, y, vx, vy, beta, delta_t) x, y, vx, vy = step_forward(x, y, vx, vy, beta, delta_t)
while y>0.0: while y>0.0:
#Propagate motion #Propagate motion
plt.plot(x_values,y_values,"r") # Plot in red plt.plot(x_values,y_values,"r") # Plot in red
plt.title("Trajectory of particle") # Add title/lables plt.title("Trajectory of particle") # Add title/lables
plt.xlabel("Distance (m)") plt.xlabel("Distance (m)")
plt.ylabel("Height") plt.ylabel("Height")
plt.grid() plt.grid()
plt.show() plt.show()
#Run interactively #Run interactively
v_initial = float(input("v0 (initial velocity) : ")) v_initial = float(input("v0 (initial velocity) : "))
theta = float(input("theta (angle from the horizontal) : ")) theta = float(input("theta (angle from the horizontal) : "))
beta = float(input("beta (drag coefficient) : ")) beta = float(input("beta (drag coefficient) : "))
delta_t = float(input("timestep (delta_t) : ")) delta_t = float(input("timestep (delta_t) : "))
Checkpoint5_1(v_initial, theta, beta, delta_t) Checkpoint5_1(v_initial, theta, beta, delta_t)
``` ```
%% Cell type:code id:a60f660b tags: %% Cell type:code id:a60f660b tags:
``` python ``` python
# Task 2 - ratio of kinetic energies # Task 2 - ratio of kinetic energies
def kinetic_energy(v): def kinetic_energy(v):
### Calculate kinetic energy ### Calculate kinetic energy
def Checkpoint5_2(v_initial, beta, delta_t): def Checkpoint5_2(v_initial, beta, delta_t):
print("K_f/K_i ratios for (v0, beta, delta_t) = ({0}, {1}, {2})".format(v_initial, beta, delta_t)) print("K_f/K_i ratios for (v0, beta, delta_t) = ({0}, {1}, {2})".format(v_initial, beta, delta_t))
... ...
#Control resolution for angle input (total # calculations = 90 * angle_res) #Control resolution for angle input (total # calculations = 90 * angle_res)
angle_res = 10 angle_res = 10
#Scan all launch angles theta from 0 to 90 #Scan all launch angles theta from 0 to 90
for theta in range(90*angle_res+1): for theta in range(90*angle_res+1):
... ...
x, y, vx, vy = set_initial(v_initial, theta/angle_res) x, y, vx, vy = set_initial(v_initial, theta/angle_res)
x, y, vx, vy = step_forward(x, y, vx, vy, beta, delta_t) x, y, vx, vy = step_forward(x, y, vx, vy, beta, delta_t)
while y>0.0: while y>0.0:
#Propagate motion #Propagate motion
#Plot graph #Plot graph
plt.plot(thetas, ke_ratios,"r") # Plot in red plt.plot(thetas, ke_ratios,"r") # Plot in red
plt.title("Ratio of Final-to-initial Kinetic Energies") # Add title/lables plt.title("Ratio of Final-to-initial Kinetic Energies") # Add title/lables
plt.xlabel(r'Launch angles $\theta$ (degrees)') plt.xlabel(r'Launch angles $\theta$ (degrees)')
plt.ylabel(r'K$_f$/K$_i$') plt.ylabel(r'K$_f$/K$_i$')
plt.xlim(0,90) plt.xlim(0,90)
plt.grid() plt.grid()
plt.show() plt.show()
#Run interactively #Run interactively
v_initial = float(input("v0 (initial velocity) : ")) v_initial = float(input("v0 (initial velocity) : "))
beta = float(input("beta (drag coefficient) : ")) beta = float(input("beta (drag coefficient) : "))
delta_t = float(input("timestep (delta_t) : ")) delta_t = float(input("timestep (delta_t) : "))
Checkpoint5_2(v_initial, beta, delta_t) Checkpoint5_2(v_initial, beta, delta_t)
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment