sciencetalk

Code: Selecteer alles

# -*- coding: utf-8 -*-
"""
Created on Tue Dec 21 14:41:36 2021
@author: Ukster
"""
import math
import numpy as np
import matplotlib.pyplot as plt

#initiele datainvoer
r=0.05        #bolstraal [m]
rhoair=1.29   #luchtdichtheid [kg/m^3]
rhobol=1545   #boldichtheid [kg/m^3]  
g=9.81        #gravity (m/s^2]
Cd=0.47       #dragcoefficient [-] 
v=50          #initiele snelheid [m/s]
ang=30        #initiele hoek [degr] met de x-as

#berekening dragfactor en bolmassa
k=1/2*rhoair*np.pi*r**2*Cd #drag factor
M=4/3*np.pi*r**3*rhobol    #bolmassa (kg)

#hoek omzetten naar radialen
#horizontale en verticale snelheidscomponent
angr=math.radians(ang)
vx=[v*np.cos(angr)]
vy=[v*np.sin(angr)]

#drag Force      
F=k*v**2

#calculate initial acceleration component x and y
ax=[-(F*np.cos(angr))/M]
ay=[-g-(F*np.sin(angr))/M]
       
#calculate for time value on the object
t=[0]
counter=0
dt=0.01
x=[0]
y=[0]

while(y[counter]>=0):
    t.append(t[counter]+dt)
    vx.append(vx[counter]+dt*ax[counter])
    vy.append(vy[counter]+dt*ay[counter])
    x.append(x[counter]+dt*vx[counter])
    y.append(y[counter]+dt*vy[counter])
    angr=math.atan(vy[counter]/vx[counter])
       
        #calculatie magnitude of new velocity
    vel=np.sqrt(vx[counter+1]**2+vy[counter+1]**2)
    F=k*vel**2
    ax.append(-(F*np.cos(angr))/M)
    ay.append(-g-(F*np.sin(angr))/M)
    
    counter=counter+1
  
plt.plot(x,y,'r')
plt.ylabel("y (m)")
plt.xlabel("x (m)")
plt.show()