1. import math
2.
4. return x*math.pi/180.0
5.
6. def arsinh(z):
7. return math.log(z+math.sqrt(z*z+1))
8.
9. # The quantity that Shouryya Ray says is constant, according to my interpreteation:
10. def ray(u,v,g,alpha):
11. return 0.5*g**2/(u**2)+0.5*alpha*g*( v*math.sqrt(u*u+v*v) / (u*u) + arsinh(math.fabs(v/u)) )
12.
13. # Numerically simulate the motion, printing out Ray's invariant at a certain interval.
14. # For testing purposes, we also tabulate the KE+PE and show the range.
15. def projectile(v0,theta,k,t,n1,n2):
16. # All units are mks.
17. # v0 = initial speed
18. # theta = angle of launch, in radians
19. # k Fair = kv^2
20. # t = amount of time to simulate
21. # n1 = number of tests to print out
22. # n2 = number of iterations between printouts
23. m = .146 # mass of a baseball
24. g=9.8 # gravitational field
25. alpha = k/m
26. vx = v0*math.cos(theta)
27. vy = v0*math.sin(theta)
28. x = 0.0
29. y = 0.0
30. dt = t/(n1*n2)
31. r_th = (v0*v0/g)*math.sin(2.0*theta)
32. heat = 0
33. for i in range(n1):
34. e = .5*m*(vx*vx+vy*vy)+m*g*y+heat
35. print("u=",vx," v=",vy," Ray=",ray(vx,vy,g,alpha)," KE+PE+heat=",e)
36. for j in range(n2):
37. theta = math.atan2(vy,vx)
38. f_air = k*(vx*vx+vy*vy)
39. #f_air = k*math.sqrt(vx*vx+vy*vy)
40. f_grav = m*g
41. fx = -f_air*math.cos(theta)
42. if fx>0.0:
43. print("fx>0, vx=",vx," vy=",vy," theta=",theta)
44. return
45. fy_air = -f_air*math.sin(theta)
46. fy = fy_air-f_grav
47. ax = fx/m
48. ay = fy/m
49. vx = vx + ax*dt
50. vy = vy + ay*dt
51. heat = heat - (fx*vx+fy_air*vy)*dt
52. if y>0.0 and y + vy*dt<0.0:
53. print("range=",x,", vs ",r_th," from theory without air resistance")
54. x = x + vx*dt
55. y = y + vy*dt
56.
57. v0 = 45.0 # m/s, a good home-run hit
58. theta = deg_to_rad(35.0) # degrees
59. k = 7.0e-4 # realistic value for a baseball
60. t = 10.0
61. n1 = 10
62. n2 = 10000
63. projectile(v0,theta,k,t,n1,n2)
64.
Success #stdin #stdout 0.48s 5852KB
stdin
Standard input is empty
stdout
u= 36.861841993  v= 25.8109396358  Ray= 0.0707591549121  KE+PE+heat= 147.825
u= 30.6462536241  v= 12.4678301768  Ray= 0.0707593974484  KE+PE+heat= 147.82340481
u= 26.6088469837  v= 1.66254890559  Ray= 0.0707595771036  KE+PE+heat= 147.822430352
u= 23.5594201658  v= -7.76184161898  Ray= 0.0859724743979  KE+PE+heat= 147.821713101
u= 20.8616382626  v= -16.0948023952  Ray= 0.104132074212  KE+PE+heat= 147.821152309
range= 120.889365695 , vs  194.171179295  from theory without air resistance
u= 18.251416064  v= -23.2425061291  Ray= 0.120666666457  KE+PE+heat= 147.820728647
u= 15.7067336398  v= -29.0889765847  Ray= 0.135385086927  KE+PE+heat= 147.820430788
u= 13.3014376668  v= -33.6520004806  Ray= 0.148676037201  KE+PE+heat= 147.82023562
u= 11.1126791141  v= -37.0751711583  Ray= 0.16096016949  KE+PE+heat= 147.820114408
u= 9.18620095669  v= -39.5647637  Ray= 0.172558265671  KE+PE+heat= 147.82004161