share/dynamics/plotdf.lisp

   1 ;; plotdf.mac - Adds a function plotdf() to Maxima, which draws a Direction
   2 ;;              Field for an ordinary 1st order differential equation,
   3 ;;              or for a system of two autonomous 1st order equations.
   4 ;;
   5 ;; Copyright (C) 2004, 2008, 2011 Jaime E. Villate <villate@fe.up.pt>
   6 ;;
   7 ;; This program is free software; you can redistribute it and/or modify
   8 ;; it under the terms of the GNU General Public License as published by
   9 ;; the Free Software Foundation; either version 2 of the License, or
  10 ;; (at your option) any later version.
  11 ;;
  12 ;; This program is distributed in the hope that it will be useful,
  13 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 ;; GNU General Public License for more details.
  16 ;;
  17 ;; You should have received a copy of the GNU General Public License
  18 ;; along with this program; if not, write to the Free Software
  19 ;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  20 ;; MA  02110-1301, USA.
  21 ;;
  22 ;; See plotdf.usg (which should come together with this program) for
  23 ;; a usage summary
  24 ;;
  25 ;; $Id: plotdf.lisp,v 1.12 2011-03-09 11:33:46 villate Exp $
  26
  27 (in-package :maxima)
  28
  29 ;; parses a plotdf option into a command-line option for tcl scripts
  30 (defun plotdf-option-to-tcl (value s1 s2)
  31   (let (v vv)
  32     (unless (and  ($listp value)
  33                   (symbolp (setq name (second value))))
  34       (merror
  35        (intl:gettext
  36         "plotdf-option-to-tcl: Expecting a symbol for the option name, found: \"~M\"") value))
  37     (case name
  38       (($xradius $yradius $xcenter $ycenter $tinitial $tstep)
  39        (setq v (check-option (cdr value) #'realp "a real number" 1))
  40        (setq value (list '(mlist) name v)))
  41       (($width $height $nsteps $versus_t)
  42        (setq v (check-option (cdr value) #'naturalp "a natural number" 1))
  43        (setq value (list '(mlist) name v)))
  44       ($trajectory_at
  45        (setq v (check-option (cdr value) #'realp "a real number" 2))
  46        (setq value (cons '(mlist) (cons name v))))
  47       ($bbox
  48        (setq v (check-option (cdr value) #'realp "a real number" 4))
  49        (setq value (cons '(mlist) (cons name v))))
  50       (($xfun $parameters $sliders $vectors $fieldlines $curves
  51               $windowtitle $xaxislabel $yaxislabel $psfile)
  52        (setq v (check-option (cdr value) #'stringp "a string" 1))
  53        (setq value (list '(mlist) name v)))
  54       ($axes
  55        (if (not (third value))
  56            (setq value '((mlist) $axes 0))
  57            (case (third value)
  58              ($x (setq value '((mlist) $axes "x")))
  59              ($y (setq value '((mlist) $axes "y")))
  60              (t (setq value '((mlist) $axes "xy"))))))
  61       ($box
  62        (if (not (third value))
  63            (setq value '((mlist) $nobox 1))
  64            (setq value '((mlist) $nobox 0))))
  65       ($direction
  66        (or
  67         (member (ensure-string (third value)) '("forward" "backward" "both") :test #'equal)
  68         (merror
  69          (intl:gettext
  70           "plotdf-option-to-tcl: direction should be forward, backward or both."))))
  71       (t (cond
  72            ((eql name s1)
  73             (setq value (cons '(mlist) (cons '$x (cddr (check-range value))))))
  74            ((eql name s2)
  75             (setq value (cons '(mlist) (cons '$y (cddr (check-range value))))))
  76            (t (merror (intl:gettext "plotdf-option-to-tcl: unknown option ~M") name)))))
  77     (setq vv (mapcar #'(lambda (a) (if (symbolp a) (ensure-string a) a)) (cdr value)))
  78     (with-output-to-string (st)
  79       (cond ((or (equal (first vv) "x") (equal (first vv) "y"))
  80              (format st "-~(~a~)center " (first vv))
  81              (format st "{~a} " (/ (+ (third vv) (second vv)) 2))
  82              (format st "-~(~a~)radius " (first vv))
  83              (format st "{~a}" (/ (- (third vv) (second vv)) 2)))
  84             (t
  85              (format st "-~(~a~) " (first vv))
  86              (format st "{~{~a~^ ~}}" (rest vv)))))))
  87
  88 ;; applies float(ev(expression, numer)) to an expression, and returns a string
  89
  90 (defun expr_to_str (fun)
  91   (mstring (mfuncall '$float (mfuncall '$ev fun '$numer))))
  92
  93 ;; plots the direction field for an ODE  dy/dx = f(x,y), or for an autonomous
  94 ;; system of 2 equations dx/dt = f(x,y), dy/dt = g(x,y)
  95 ;;
  96 (defun $plotdf (ode &rest options)
  97   (let (file cmd (opts " ") (s1 '$x) (s2 '$y) vars)
  98     (unless ($listp ode) (setf ode `((mlist) ,ode)))
  99
 100     ;; if the variables are not x and y, their names must be given right
 101     ;; after the expression for the ode's
 102     (when
 103       (and (listp (car options)) (= (length (car options)) 3)
 104             (or (symbolp (cadar options)) ($subvarp (cadar options)))
 105             (or (symbolp (caddar options)) ($subvarp (caddar options))))
 106       (setq s1 (cadar options))
 107       (setq s2 (caddar options))
 108       (setq options (cdr options)))
 109
 110     ;; parse options and copy them to string opts
 111     (cond (options
 112            (dolist (v options)
 113              (setq opts (concatenate 'string opts " "
 114                                   (plotdf-option-to-tcl v s1 s2))))))
 115
 116     (unless (search "-xaxislabel " opts)
 117       (setq opts (concatenate 'string opts " -xaxislabel " (ensure-string s1))))
 118     (unless (search "-yaxislabel " opts)
 119       (setq opts (concatenate 'string opts " -yaxislabel " (ensure-string s2))))
 120
 121     ;; check that the expressions given contain only the axis variables
 122     (setq vars (cdr (mfuncall '$listofvars ode)))
 123     (unless (search "-parameters " opts)
 124       (dolist (var vars)
 125         (unless (or (eq var s1) (eq var s2))
 126           (merror
 127            (intl:gettext "plotdf: expression(s) given can only depend on ~M and ~M~%Found extra variable ~M") s1 s2 var))))
 128
 129     ;; substitute $x by s1 and $y by s2
 130     (defun subxy (expr)
 131       (if (listp expr)
 132           (mapcar #'subxy expr)
 133         (cond ((eq expr s1) '$x) ((eq expr s2) '$y) (t expr))))
 134     (setf ode (mapcar #'subxy ode))
 135
 136     ;; parse the differential equation expressions
 137     (case (length ode)
 138           (3 (setq cmd (concatenate 'string " -dxdt \""
 139                                     (expr_to_str (second ode)) "\" -dydt \""
 140                                     (expr_to_str (third ode)) "\"")))
 141           (2 (setq cmd (concatenate 'string " -dydx \""
 142                                     (expr_to_str (second ode)) "\"")))
 143           (t (merror
 144               (intl:gettext "plotdf: first argument must be either an expression or a list with two expressions."))))
 145     (let (data)
 146       (setq file (plot-file-path (format nil "~a.xmaxima" (random-name 16))))
 147       (cond ($show_openplot
 148              (setq data (format nil "plotdf ~a ~a~%" cmd opts))
 149              (with-open-file
 150               (fl
 151                #+sbcl (sb-ext:native-namestring file)
 152                #-sbcl file
 153                :direction :output :if-exists :supersede)
 154               (princ data fl))
 155              ($system (concatenate 'string *maxima-prefix*
 156                                    (if (string= *autoconf-windows* "true")
 157                                        "\\bin\\" "/bin/")
 158                                    $xmaxima_plot_command)
 159                       #-(or (and sbcl win32) (and sbcl win64) (and ccl windows))
 160                       (format nil " ~s &" file)
 161                       #+(or (and sbcl win32) (and sbcl win64) (and ccl windows))
 162                       file)
 163              (list '(mlist) file))
 164             (t
 165              (setq data (format nil "{plotdf ~a ~a}" cmd opts))
 166              (princ data) "")))))
 167
 168 ;; plot equipotential curves for a scalar field f(x,y)
 169 (defun $ploteq (fun &rest options)
 170
 171   (let (file cmd mfun (opts " ") (s1 '$x) (s2 '$y))
 172     (setf mfun `((mtimes) -1 ,fun))
 173     ;; if the variables are not x and y, their names must be given right
 174     ;; after the expression for the ode's
 175     (when
 176       (and (listp (car options)) (= (length (car options)) 3)
 177             (or (symbolp (cadar options)) ($subvarp (cadar options)))
 178             (or (symbolp (caddar options)) ($subvarp (caddar options))))
 179       (setq s1 (cadar options))
 180       (setq s2 (caddar options))
 181       (setq options (cdr options)))
 182     (defun subxy (expr)
 183       (if (listp expr)
 184           (mapcar #'subxy expr)
 185           (cond ((eq expr s1) '$x) ((eq expr s2) '$y) (t expr))))
 186     (setf mfun (mapcar #'subxy mfun))
 187 ;; the next two lines should take into account parameters given in the options
 188 ;;    (if (delete '$y (delete '$x (rest (mfuncall '$listofvars ode))))
 189 ;;        (merror "The equation(s) can depend only on 2 variable which must be specified!"))
 190     (setq cmd (concatenate 'string " -dxdt \""
 191                            (expr_to_str (mfuncall '$diff mfun '$x))
 192                            "\" -dydt \""
 193                            (expr_to_str (mfuncall '$diff mfun '$y))
 194                            "\" "))
 195
 196     ;; parse options and copy them to string opts
 197     (cond (options
 198            (dolist (v options)
 199              (setq opts (concatenate 'string opts " "
 200                                   (plotdf-option-to-tcl v s1 s2))))))
 201
 202     (unless (search "-vectors " opts)
 203       (setq opts (concatenate 'string opts " -vectors {}")))
 204     (unless (search "-fieldlines " opts)
 205       (setq opts (concatenate 'string opts " -fieldlines {}")))
 206     (unless (search "-curves " opts)
 207       (setq opts (concatenate 'string opts " -curves {red}")))
 208     (unless (search "-windowtitle " opts)
 209       (setq opts (concatenate 'string opts " -windowtitle {Ploteq}")))
 210     (unless (search "-xaxislabel " opts)
 211       (setq opts (concatenate 'string opts " -xaxislabel " (ensure-string s1))))
 212     (unless (search "-yaxislabel " opts)
 213       (setq opts (concatenate 'string opts " -yaxislabel " (ensure-string s2))))
 214
 215     (let (data)
 216       (setq file (plot-file-path (format nil "~a.xmaxima" (random-name 16))))
 217       (cond ($show_openplot
 218              (setq data (format nil "plotdf ~a ~a~%" cmd opts))
 219              (with-open-file
 220               (fl
 221                #+sbcl (sb-ext:native-namestring file)
 222                #-sbcl file
 223                :direction :output :if-exists :supersede)
 224               (princ data fl))
 225              ($system (concatenate 'string *maxima-prefix*
 226                                    (if (string= *autoconf-windows* "true")
 227                                        "\\bin\\" "/bin/")
 228                                    $xmaxima_plot_command)
 229                       #-(or (and sbcl win32) (and sbcl win64) (and ccl windows))
 230                       (format nil " ~s &" file)
 231                       #+(or (and sbcl win32) (and sbcl win64) (and ccl windows))
 232                       file)
 233              (list '(mlist) file))
 234             (t
 235              (setq data (format nil "{plotdf ~a ~a}" cmd opts))
 236              (princ data) "")))))