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
      k1 = ax*dt
      k2 = (ax+0.5*k1)*dt
      k3 = (ax+0.5*k2)*dt
      k4 = (ax+k3)*dt
      vx = vx + 1/6*(k1+2*k2+2*k3+k4)
      k1 = ay*dt
      k2 = (ay+0.5*k1)*dt
      k3 = (ay+0.5*k2)*dt
      k4 = (ay+k3)*dt
      vy = vy + 1/6*(k1+2*k2+2*k3+k4)
      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")
      k1 = vx*dt
      k2 = (vx+0.5*k1)*dt
      k3 = (vx+0.5*k2)*dt
      k4 = (vx+k3)*dt
      x = x + 1/6*(k1+2*k2+2*k3+k4)
      k1 = vy*dt
      k2 = (vy+0.5*k1)*dt
      k3 = (vy+0.5*k2)*dt
      k4 = (vy+k3)*dt
      y = y + 1/6*(k1+2*k2+2*k3+k4)
 
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+MCwgdng9Iix2eCwiIHZ5PSIsdnksIiB0aGV0YT0iLHRoZXRhKQogICAgICAgIHJldHVybgogICAgICBmeV9haXIgPSAtZl9haXIqbWF0aC5zaW4odGhldGEpCiAgICAgIGZ5ID0gZnlfYWlyLWZfZ3JhdgogICAgICBheCA9IGZ4L20KICAgICAgYXkgPSBmeS9tCiAgICAgIGsxID0gYXgqZHQKICAgICAgazIgPSAoYXgrMC41KmsxKSpkdAogICAgICBrMyA9IChheCswLjUqazIpKmR0CiAgICAgIGs0ID0gKGF4K2szKSpkdAogICAgICB2eCA9IHZ4ICsgMS82KihrMSsyKmsyKzIqazMrazQpCiAgICAgIGsxID0gYXkqZHQKICAgICAgazIgPSAoYXkrMC41KmsxKSpkdAogICAgICBrMyA9IChheSswLjUqazIpKmR0CiAgICAgIGs0ID0gKGF5K2szKSpkdAogICAgICB2eSA9IHZ5ICsgMS82KihrMSsyKmsyKzIqazMrazQpCiAgICAgIGhlYXQgPSBoZWF0IC0gKGZ4KnZ4K2Z5X2Fpcip2eSkqZHQKICAgICAgaWYgeT4wLjAgYW5kIHkgKyB2eSpkdDwwLjA6CiAgICAgICAgcHJpbnQoInJhbmdlPSIseCwiLCB2cyAiLHJfdGgsIiBmcm9tIHRoZW9yeSB3aXRob3V0IGFpciByZXNpc3RhbmNlIikKICAgICAgazEgPSB2eCpkdAogICAgICBrMiA9ICh2eCswLjUqazEpKmR0CiAgICAgIGszID0gKHZ4KzAuNSprMikqZHQKICAgICAgazQgPSAodngrazMpKmR0CiAgICAgIHggPSB4ICsgMS82KihrMSsyKmsyKzIqazMrazQpCiAgICAgIGsxID0gdnkqZHQKICAgICAgazIgPSAodnkrMC41KmsxKSpkdAogICAgICBrMyA9ICh2eSswLjUqazIpKmR0CiAgICAgIGs0ID0gKHZ5K2szKSpkdAogICAgICB5ID0geSArIDEvNiooazErMiprMisyKmszK2s0KQoKdjAgPSA0NS4wICMgbS9zLCBhIGdvb2QgaG9tZS1ydW4gaGl0CnRoZXRhID0gZGVnX3RvX3JhZCgzNS4wKSAjIGRlZ3JlZXMKayA9IDcuMGUtNCAjIHJlYWxpc3RpYyB2YWx1ZSBmb3IgYSBiYXNlYmFsbAp0ID0gMTAuMApuMSA9IDEwCm4yID0gMTAwMDAKcHJvamVjdGlsZSh2MCx0aGV0YSxrLHQsbjEsbjIpCg==