;;; Lisplab, level3-fft-blas.lisp
;;; Runge-Kuttas fourth order methods for solving 
;;; differential equations.

;;; Copyright (C) 2009 Joern Inge Vestgaarden
;;;
;;; 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 2 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.

(in-package :lisplab)

(defun rk4 (func y x 
	    &key 
	    (step)
	    (update)
	    (test (constantly nil)))
  "The fourth order Runge-Kutta"
  (flet ((rk4-step (func y x h)
           (let* ((k1 (.* h (funcall func y x)))
                  (k2 (.* h (funcall func (.+ y (.* 0.5 k1)) (+ x (* 0.5 h)))))
                  (k3 (.* h (funcall func (.+ y (.* 0.5 k2)) (+ x (* 0.5 h)))))
                  (k4 (.* h (funcall func (.+ y k3) (+ x h)))))
	     (.+ y (.* 1/6 (.+ k1 (.+ (.* 2 k2) (.+ (.* 2 k3) k4))))))))
    (do* ((h (funcall step y x) (funcall step y x)))
	 ((not (funcall test y x)) (values y x))
      	 (setf y (rk4-step func y x h)
	       x (+ x h))
      (when update
        (funcall update y x)))))