def newton alpha if alpha or alpha pi print Bad angle return def fx re

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def newton(v, alpha):
if alpha <= 0 or alpha >= pi/2:
print("Bad angle")
return .0
def fx(t, v):
return -(beta*v**2)
def fy(t, v):
return -copysign(1, v)*beta*v**2-g
dvx, dvy = RK4(fx), RK4(fy)
dt = 0.0001
vx, vy = v*cos(alpha), v*sin(alpha)
t = dt
x_old = y_old = 0
x = vx*dt
y = vy*dt
while (y > 0):
t += dt
vx_new = vx + dvx(t, vx, dt)
vy_new = vy + dvy(t, vy, dt)
x_old = x
x += (vx_new + vx) * dt / 2
y_old = y
y += (vy_new + vy) * dt / 2
vx = vx_new
vy = vy_new
if y_old == 0:
return x_old
return (x_old*(y_old - y) - y_old*(x_old-x))/(y_old - y)