;; Function $rk implements the 4th order Runge-Kutta numerical method.
;;  exprs:   an expression or maxima list with n expressions
;;  vars:    name of (or list of names of) the independent variable(s)
;;  initial: initial value for the variable (or list of initial values)
;;  domain:  maxima list with four elements (name of the independent
;;           variable, its initial and final values and its increment)
;;; derived from complex-dynamics.lisp which was
;;;;; Copyright (C) 2006-2013 Jaime E. Villate <villate@fe.up.pt>

;;; trying to put double-float declarations in every place? 
;;; hard -- work must be done by coerce-float-fun defined in plot.lisp.
;;; presumably in concert with proper declaration by user's  mode_declare().

;;; Mostly I changed constants here to double-floats, and compiled it all (in SBCL)
;;; I renamed the program rkf   from rk.

;;; last(rk(t-x^2,x,1.0,[t,0,8,0.1d0])); (in complex-dynamics) took 0.313 seconds and 24 megabytes
;;; last(rkf(t-x^2,x,1.0,[t,0,8,0.1d0]));  took 0.016 seconds and 0.698 megabytes

(defun $rkf (exprs vars initial domain ;; float version
            &aux d u fun k1 k2 k3 k4 r1 r2 r3 traj r
              (it (mapcar #'coerce-float (cddr domain))))
  (unless ($listp exprs) (setq exprs `((mlist) ,exprs)))
  (unless ($listp initial) (setq initial `((mlist) ,initial)))
  (unless ($listp vars) (setq vars `((mlist) ,vars)))
  (dolist (var (cdr vars))
    (unless (symbolp var)
      (merror (intl:gettext "rk: variable name expected; found: ~M") var)))
  (unless (symbolp (cadr domain))
    (merror (intl:gettext "rk: variable name expected; found: ~M")
            (cadr domain)))
  
  (setq vars (concatenate 'list '((mlist)) (list (cadr domain)) (cdr vars)))
  
  (setq r (concatenate 'list `(,(car it)) (mapcar #'coerce-float (cdr initial))));; now r contains double-float numbers 
  (setq fun (mapcar #'(lambda (x) (coerce-float-fun x vars)) (cdr exprs)))
  ;;now fun is a list of functions, one for each expression in expr
  ;; these functions are, I think, not compiled, nor are they
  ;; fully declared as float inputs etc.
  
  (setq d (/ (- (cadr it) (car it)) (caddr it)))
  (setq traj (list (cons '(mlist) r)))
  (do ((m 1 (1+ m))) ((> m d))
    (ignore-errors
      (setq k1 (mapcar #'(lambda (x) (apply x r)) fun))
      (setq r1 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* 0.5d0 (caddr it) x)) k1)))
      (push (+ (car r) (/ (caddr it) 2)) r1)
      (setq k2 (mapcar #'(lambda (x) (apply x r1)) fun))
      (setq r2 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* 0.5d0 (caddr it) x)) k2)))
      (push (+ (car r) (/ (caddr it) 2)) r2)
      (setq k3 (mapcar #'(lambda (x) (apply x r2)) fun))
      (setq r3 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (caddr it) x)) k3)))
      (push (+ (car r) (caddr it)) r3)
      (setq k4 (mapcar #'(lambda (x) (apply x r3)) fun))
      (setq u (map 'list #'+
                   (mapcar #'(lambda (x) (* #.(/ 1.0d0 6.0d0) x)) k1)
                   (mapcar #'(lambda (x) (* #.(/ 1.0d0 3.0d0) x)) k2)
                   (mapcar #'(lambda (x) (* #.(/ 1.0d0 3.0d0) x)) k3)
                   (mapcar #'(lambda (x) (* #.(/ 1.0d0 6.0d0) x)) k4)))
      (setq r 
            (concatenate 'list
                         `(,(+ (car it) (* m (caddr it))))
                         (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (caddr it) x)) u))))
      (push (cons '(mlist) r) traj)))
  (when (< (car r) (cadr it))
    (let ((s (- (cadr it) (car r))))
      (ignore-errors
        (setq k1 (mapcar #'(lambda (x) (apply x r)) fun))
        (setq r1 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* 0.5d0 s x)) k1)))
        (push (+ (car r) (* 0.5d0 s)) r1)
        (setq k2 (mapcar #'(lambda (x) (apply x r1)) fun))
        (setq r2 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* 0.5d0 s x)) k2)))
        (push (+ (car r)(* 0.5d0 s)) r2)
        (setq k3 (mapcar #'(lambda (x) (apply x r2)) fun))
        (setq r3 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* s x)) k3)))
        (push (+ (car r) s) r3)
        (setq k4 (mapcar #'(lambda (x) (apply x r3)) fun))
        (setq u (map 'list #'+
		     (mapcar #'(lambda (x) (* #.(/ 1.0d0 6.0d0) x)) k1)
                   (mapcar #'(lambda (x) (* #.(/ 1.0d0 3.0d0) x)) k2)
                   (mapcar #'(lambda (x) (* #.(/ 1.0d0 3.0d0) x)) k3)
                   (mapcar #'(lambda (x) (* #.(/ 1.0d0 6.0d0) x)) k4)))

        (setq r 
              (concatenate 'list
                           `(,(cadr it))
                           (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* s x)) u)))))
        (push (cons '(mlist) r) traj)))
  (cons '(mlist) (nreverse traj)))