import math
def deg_to_rad(x):
return x*math.pi/180.0
def arsinh(z):
return math.log(z+math.sqrt(z*z+1))
# The quantity that Shouryya Ray says is constant, according to my interpreteation:
def ray(u,v,g,alpha):
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)) )
# Numerically simulate the motion, printing out Ray's invariant at a certain interval.
# For testing purposes, we also tabulate the KE+PE and show the range.
def projectile(v0,theta,k,t,n1,n2):
# All units are mks.
# v0 = initial speed
# theta = angle of launch, in radians
# k Fair = kv^2
# t = amount of time to simulate
# n1 = number of tests to print out
# n2 = number of iterations between printouts
m = .146 # mass of a baseball
g=9.8 # gravitational field
alpha = k/m
vx = v0*math.cos(theta)
vy = v0*math.sin(theta)
x = 0.0
y = 0.0
dt = t/(n1*n2)
r_th = (v0*v0/g)*math.sin(2.0*theta)
heat = 0
for i in range(n1):
e = .5*m*(vx*vx+vy*vy)+m*g*y+heat
print("u=",vx," v=",vy," Ray=",ray(vx,vy,g,alpha)," KE+PE+heat=",e)
for j in range(n2):
theta = math.atan2(vy,vx)
f_air = k*(vx*vx+vy*vy)
#f_air = k*math.sqrt(vx*vx+vy*vy)
f_grav = m*g
fx = -f_air*math.cos(theta)
if fx>0.0:
print("fx>0, vx=",vx," vy=",vy," theta=",theta)
return
fy_air = -f_air*math.sin(theta)
fy = fy_air-f_grav
ax = fx/m
ay = fy/m
vx = vx + ax*dt
vy = vy + ay*dt
heat = heat - (fx*vx+fy_air*vy)*dt
if y>0.0 and y + vy*dt<0.0:
print("range=",x,", vs ",r_th," from theory without air resistance")
x = x + vx*dt
y = y + vy*dt
v0 = 45.0 # m/s, a good home-run hit
theta = deg_to_rad(35.0) # degrees
k = 7.0e-4 # realistic value for a baseball
t = 10.0
n1 = 10
n2 = 10000
projectile(v0,theta,k,t,n1,n2)
aW1wb3J0IG1hdGgKCmRlZiBkZWdfdG9fcmFkKHgpOgogIHJldHVybiB4Km1hdGgucGkvMTgwLjAKCmRlZiBhcnNpbmgoeik6CiAgcmV0dXJuIG1hdGgubG9nKHorbWF0aC5zcXJ0KHoqeisxKSkKCiMgVGhlIHF1YW50aXR5IHRoYXQgU2hvdXJ5eWEgUmF5IHNheXMgaXMgY29uc3RhbnQsIGFjY29yZGluZyB0byBteSBpbnRlcnByZXRlYXRpb246CmRlZiByYXkodSx2LGcsYWxwaGEpOgogIHJldHVybiAwLjUqZyoqMi8odSoqMikrMC41KmFscGhhKmcqKCB2Km1hdGguc3FydCh1KnUrdip2KSAvICh1KnUpICsgYXJzaW5oKG1hdGguZmFicyh2L3UpKSApCgojIE51bWVyaWNhbGx5IHNpbXVsYXRlIHRoZSBtb3Rpb24sIHByaW50aW5nIG91dCBSYXkncyBpbnZhcmlhbnQgYXQgYSBjZXJ0YWluIGludGVydmFsLgojIEZvciB0ZXN0aW5nIHB1cnBvc2VzLCB3ZSBhbHNvIHRhYnVsYXRlIHRoZSBLRStQRSBhbmQgc2hvdyB0aGUgcmFuZ2UuCmRlZiBwcm9qZWN0aWxlKHYwLHRoZXRhLGssdCxuMSxuMik6CiAgIyBBbGwgdW5pdHMgYXJlIG1rcy4KICAjIHYwICAgID0gaW5pdGlhbCBzcGVlZAogICMgdGhldGEgPSBhbmdsZSBvZiBsYXVuY2gsIGluIHJhZGlhbnMKICAjIGsgICAgIEZhaXIgPSBrdl4yCiAgIyB0ICAgICA9IGFtb3VudCBvZiB0aW1lIHRvIHNpbXVsYXRlCiAgIyBuMSAgICA9IG51bWJlciBvZiB0ZXN0cyB0byBwcmludCBvdXQKICAjIG4yICAgID0gbnVtYmVyIG9mIGl0ZXJhdGlvbnMgYmV0d2VlbiBwcmludG91dHMKICBtID0gLjE0NiAjIG1hc3Mgb2YgYSBiYXNlYmFsbAogIGc9OS44ICMgZ3Jhdml0YXRpb25hbCBmaWVsZAogIGFscGhhID0gay9tCiAgdnggPSB2MCptYXRoLmNvcyh0aGV0YSkKICB2eSA9IHYwKm1hdGguc2luKHRoZXRhKQogIHggPSAwLjAKICB5ID0gMC4wCiAgZHQgPSB0LyhuMSpuMikKICByX3RoID0gKHYwKnYwL2cpKm1hdGguc2luKDIuMCp0aGV0YSkKICBoZWF0ID0gMAogIGZvciBpIGluIHJhbmdlKG4xKToKICAgIGUgPSAuNSptKih2eCp2eCt2eSp2eSkrbSpnKnkraGVhdAogICAgcHJpbnQoInU9Iix2eCwiIHY9Iix2eSwiIFJheT0iLHJheSh2eCx2eSxnLGFscGhhKSwiIEtFK1BFK2hlYXQ9IixlKQogICAgZm9yIGogaW4gcmFuZ2UobjIpOgogICAgICB0aGV0YSA9IG1hdGguYXRhbjIodnksdngpCiAgICAgIGZfYWlyID0gayoodngqdngrdnkqdnkpCiAgICAgICNmX2FpciA9IGsqbWF0aC5zcXJ0KHZ4KnZ4K3Z5KnZ5KQogICAgICBmX2dyYXYgPSBtKmcKICAgICAgZnggPSAtZl9haXIqbWF0aC5jb3ModGhldGEpCiAgICAgIGlmIGZ4PjAuMDoKICAgICAgICBwcmludCgiZng+MCwgdng9Iix2eCwiIHZ5PSIsdnksIiB0aGV0YT0iLHRoZXRhKQogICAgICAgIHJldHVybgogICAgICBmeV9haXIgPSAtZl9haXIqbWF0aC5zaW4odGhldGEpCiAgICAgIGZ5ID0gZnlfYWlyLWZfZ3JhdgogICAgICBheCA9IGZ4L20KICAgICAgYXkgPSBmeS9tCiAgICAgIHZ4ID0gdnggKyBheCpkdAogICAgICB2eSA9IHZ5ICsgYXkqZHQKICAgICAgaGVhdCA9IGhlYXQgLSAoZngqdngrZnlfYWlyKnZ5KSpkdAogICAgICBpZiB5PjAuMCBhbmQgeSArIHZ5KmR0PDAuMDoKICAgICAgICBwcmludCgicmFuZ2U9Iix4LCIsIHZzICIscl90aCwiIGZyb20gdGhlb3J5IHdpdGhvdXQgYWlyIHJlc2lzdGFuY2UiKQogICAgICB4ID0geCArIHZ4KmR0CiAgICAgIHkgPSB5ICsgdnkqZHQKCnYwID0gNDUuMCAjIG0vcywgYSBnb29kIGhvbWUtcnVuIGhpdAp0aGV0YSA9IGRlZ190b19yYWQoMzUuMCkgIyBkZWdyZWVzCmsgPSA3LjBlLTQgIyByZWFsaXN0aWMgdmFsdWUgZm9yIGEgYmFzZWJhbGwKdCA9IDEwLjAKbjEgPSAxMApuMiA9IDEwMDAwCnByb2plY3RpbGUodjAsdGhldGEsayx0LG4xLG4yKQo=