fork download
  1. import math
  2.  
  3. def deg_to_rad(x):
  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