import
  io(open),
  ipl.functional(_a, lambda)

$encoding UTF-8

procedure main ()
  local f, plot, ty
  local data2, data5, data10

  # Newton's cooling law, f(t,Temp) = -0.07*(Temp-20)
  f := lambda { -0.07 * (_a[2] - 20.0) }

  data2 := euler_method (f, 100, 0, 100, 2)
  data5 := euler_method (f, 100, 0, 100, 5)
  data10 := euler_method (f, 100, 0, 100, 10)

  plot := open ("gnuplot", "pw")
  plot.write ("set encoding utf8")
  plot.write ("set term png size 1000,750 font 'Fanwood Text,18'")
  plot.write ("set output 'newton-cooling-OI.png'")
  plot.write ("set grid")
  plot.write (u"set title 'Newton’s Law of Cooling'")
  plot.write ("set xlabel 'Elapsed time (seconds)'")
  plot.write ("set ylabel 'Temperature (Celsius)'")
  plot.write ("set xrange [0:100]")
  plot.write ("set yrange [15:100]")
  plot.write ("y(x) = 20.0 + (80.0 * exp (-0.07 * x))")
  plot.write ("plot y(x) with lines title 'Analytic solution', \\")
  plot.write ("     '-' with linespoints title 'Euler method, step size 2s', \\")
  plot.write ("     '-' with linespoints title 'Step size 5s', \\")
  plot.write ("     '-' with linespoints title 'Step size 10s'")
  every plot.write (ty := !data2 & ty[1] || " " || ty[2])
  plot.write ("e")
  every plot.write (ty := !data5 & ty[1] || " " || ty[2])
  plot.write ("e")
  every plot.write (ty := !data10 & ty[1] || " " || ty[2])
  plot.write ("e")
  plot.close()
end

# Approximate y(t) in dy/dt=f(t,y), y(a)=y0, t going from a to b with
# positive step size h.
procedure euler_method (f, y0, a, b, h)
  local t, y, results

  t := a
  y := y0
  results := [[t, y]]
  while t + h <= b do
  {
    y +:= h * f(t, y)
    t +:= h
    put (results, [t, y])
  }
  return results
end