# ballbound.jl --- 投げたボールのバウンド, 空気抵抗ありのシミュレーション

# Runge-Kutta 法の1ステップ
function rungekutta(f,t,x,dt)
  k1=dt*f(t,x)
  k2=dt*f(t+dt/2, x+k1/2)
  k3=dt*f(t+dt/2, x+k2/2)
  k4=dt*f(t+dt, x+k3)
  x + (k1 + 2 * k2 + 2 * k3 + k4) / 6
end

function f(t,x)
  Gamma=1.0
  m=100.0
  g=9.8
  y=similar(x)
  y[1] = x[3]
  y[2] = x[4]
  y[3] = -Gamma/m*x[3]
  y[4] = -g-Gamma/m*x[4]
  y
end

function ballbound(n=1000) # 1000等分くらい
  t=0.0
  v0=50.0
  theta=50.0
  x=[0.0,0.0,v0*cos(theta*pi/180),v0*sin(theta*pi/180)]
  println("t x")
  println("$t $x")
  Tmax=20.0
  dt=Tmax/n
  println("Tmax=$Tmax, dt=$dt")
  s=@sprintf "%f %f %f %f %f\n" t x[1] x[2] x[3] x[4]
  print(s)
  of = open("ballbound.dat","w")
  print(of,s)
  for i=1:n
    # x=euler(f,t,x,dt)
    x=rungekutta(f,t,x,dt)
    if x[2]<0
      x[2] = - x[2]
      x[4] = - x[4]
    end
    t=i*dt
    s=@sprintf "%f %f %f %f %f\n" t x[1] x[2] x[3] x[4]
    print(s)
    print(of,s)
  end
  close(of)
  # 以下は工事中
  of = open("ballbound.gp","w")
  println(of,"plot \"ballbound.dat\" using 2:3 with lp")
  close(of)
  run(`gnuplot ballbound.gp`)
end