;; ODE system setup ;; Liam Healy, Sun Apr 15 2007 - 14:19 ;; Time-stamp: <2009-01-24 18:23:33EST ode-system.lisp> ;; \$Id\$ (in-package :gsl) (cffi:defcstruct ode-system ; gsl_odeiv_system ;; See /usr/include/gsl/gsl_odeiv.h "The definition of an ordinary differential equation system for GSL." (function :pointer) (jacobian :pointer) (dimension sizet) (parameters :pointer)) (export '(make-ode-functions with-ode-integration)) (defmacro make-ode-functions (name jacobian dimension) "Setup functions for ODE integrators. The variable name is used as the name of the The CL functions name and jacobian should be defined previously with defuns. The function and Jacobian arguments are the same: time and dimension dependent variables as scalars. The function returns dimension multiple values corresponding to the derivatives of the dependent variables, f(t). The Jacobian returns dimension^2 values corresponding to the partial derivatives of each of the functions f with respect to each of the dependent variables y, and dimension values corresponding to the derivatives of f with respect to t." ;; set return values ;; Possible future improvements: take/set arrays easily, allow lambdas instead of named functions. (with-unique-names (solverfn solverdf) `(progn (defmcallback ,solverfn :success-failure (:double (:double ,dimension) (:set :double ,dimension)) nil nil ,name) (defmcallback ,solverdf :success-failure (:double (:double ,dimension) (:set :double ,(expt dimension 2)) (:set :double ,dimension)) nil nil ,jacobian) (defcbstruct (,solverfn function ,solverdf jacobian) ode-system ((dimension ,dimension)))))) (defmacro with-ode-integration ((time step-size dependent dimensions &optional (stepper '*step-rk8pd*) (absolute-error 1.0d-6) (relative-error 0.0d0)) &body body) "Environment for integration of ordinary differential equations. The variables time and step-size will become C doubles in the body; to convert back, use #'dcref The dependent variable may be specified as a list being the same as the first argument to with-c-double." (let ((ctime (make-symbol "CTIME")) (cstep (make-symbol "CSTEP"))) `(let ((stepper (make-ode-stepper ,stepper ,dimensions)) (control (new-y-control ,absolute-error ,relative-error)) (evolve (make-ode-evolution ,dimensions))) (unwind-protect (cffi:with-foreign-objects ((,(if (listp dependent) (first dependent) dependent) :double ,dimensions) (,ctime :double) (,cstep :double)) (setf (dcref ,cstep) ,step-size (dcref ,ctime) ,time) ,(if (listp dependent) `(with-c-double ,dependent (symbol-macrolet ((,time ,ctime) (,step-size ,cstep)) ,@body)) `(symbol-macrolet ((,time ,ctime) (,step-size ,cstep)) ,@body))) (free-control control)))))