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)