# ballbound-v3.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 ballbound3(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")
  p=open(pipeline(`gnuplot --persist`; stderr=Pipe()), "r+")
  println(p.in, "plot '-' with lines title \"n=$n\"");

  println(p.in, "$(x[1]) $(x[2])")
  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
    # @printf(p.in, "%f %f\n", x[1], x[2])
    println(p.in, "$(x[1]) $(x[2])")
  end
  println(p.in, "e")
  flush(p.in)
  println(p.in, "set term png");
  println(p.in, "set output \"ballbound3.png\"");
  println(p.in, "replot");
  println(p.in, "quit")
  flush(p.in)
  close(p)
end