
# -- num/rkf45.lua.in
# -- 
# -- Copyright (C) 2009-2011 Francesco Abbate
# -- 
# -- This program is free software; you can redistribute it and/or modify
# -- it under the terms of the GNU General Public License as published by
# -- the Free Software Foundation; either version 3 of the License, or (at
# -- your option) any later version.
# -- 
# -- This program is distributed in the hope that it will be useful, but
# -- WITHOUT ANY WARRANTY; without even the implied warranty of
# -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# -- General Public License for more details.
# -- 
# -- You should have received a copy of the GNU General Public License
# -- along with this program; if not, write to the Free Software
# -- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# --

# -- Adapted from the GSL Library, version 1.14

# -- Runge-Kutta-Fehlberg 4(5)

# -- Reference eg. Hairer, E., Norsett S.P., Wanner, G. Solving ordinary
# -- differential equations I, Nonstiff Problems, 2nd revised edition,
# -- Springer, 2000.

local abs, max, min = math.abs, math.max, math.min

# order = 5

$(include 'ode-defs')

# AH = { '1/4', '3/8', '12/13', '1', '1/2' }

# B = { 0 }
# B[2] = { '1/4' }
# B[3] = { '3/32', '9/32' }
# B[4] = { '1932/2197', '-7200/2197', '7296/2197'}
# B[5] = { '8341/4104', '-32832/4104', '29440/4104', '-845/4104'}
# B[6] = { '-6080/20520', '41040/20520', '-28352/20520', '9295/20520', '-5643/20520'}

# C = {'902880/7618050', 
#      '0', 
#      '3953664/7618050', 
#      '3855735/7618050', 
#      '-1371249/7618050', 
#      '277020/7618050' }

# -- These are the differences of fifth and fourth order coefficients
# -- for error estimation

# EC = { '1/360', '0', '-128/4275', '-2197/75240', '1/50', '2/55' }

# y_err_only = (a_dydt == 0)

local function rkf45_step(s, t1)
   local t, h, f = s.t, s.h, s.f
   local s_y, s_dydt = s.y, s.dydt
   local hadj, inc

   local $(VL'y')
   local $(VL'k1') = $(AL's_dydt.data')

   if t < t1 and t + h > t1 then h = t1 - t end

   while h > 0 do
      $(VL'y') = $(AL's_y.data')
      local rmax = 0

      do
      local $(VL'ytmp')

#     for S = 2, 6 do
#        for i = 0, N-1 do
            ytmp_$(i) = y_$(i) + h * ($(KCONV(B[S], S-1, i)))
#        end
         local $(VLI('k', S)) = f(t + $(AH[S-1]) * h, $(VL'ytmp'))
#     end

      local di
#     for i = 0, N-1 do
         di = $(KCONV(C, 6, i))
         y_$(i) = y_$(i) + h * di
#     end

#     if not y_err_only then
         local $(VL'dydt') = f(t + h, $(VL'y'))
#        for i = 0, N-1 do
            s_dydt.data[$(i)] = dydt_$(i)
#        end
#     end

      do
      local yerr, r, d0

#     for i = 0, N-1 do
         yerr = h * ($(KCONV(EC, 6, i)))
#     if y_err_only then
         d0 = $(eps_rel) * ($(a_y) * abs(y_$(i))) + $(eps_abs)
#     else
         d0 = $(eps_rel) * ($(a_y) * abs(y_$(i)) + $(a_dydt) * abs(h * dydt_$(i))) + $(eps_abs)
#     end
         r = abs(yerr) / abs(d0)
         rmax = max(r, rmax)
#     end
      end
      end

      hadj, inc = hadjust(rmax, h)
      if inc >= 0 then break end
      h = hadj
   end

#  if y_err_only then
      local $(VL'dydt') = f(t + h, $(VL'y'))
#     for i = 0, N-1 do
         s_dydt.data[$(i)] = dydt_$(i)
#     end
#  end

#  for i = 0, N-1 do
      s_y.data[$(i)] = y_$(i)
#  end
   s.t = t + h
   s.h = hadj
end

return {new= ode_new, init= ode_init, evolve= ode_evolve, step= rkf45_step}
