Update package export list

.gitignore

@@ -0,0 +1 @@
+maxima-install

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "maxima"]
+	path = maxima
+	url = https://git.code.sf.net/p/maxima/code
+[submodule "maxima-interface"]
+	path = maxima-interface
+	url = https://git.sr.ht/~jmbr/maxima-interface

README

@@ -0,0 +1,45 @@
+				 HAM: A PHYSICS SIMULATION FRAMEWORK
+
+Ham is a physics simulation framework that takes a Lagrangian of a
+system and a description of how to render it and produces a real-time,
+Hamiltonian-based simulation of the system.  It uses Maxima for the
+symbolic mathematics and SDL2 for the graphics.
+
+For example, for a simple harmonic oscillator:
+
+	(define-simulation simple-harmonic-oscillator
+	  :lagrangian (- (* 1/2 m (^ (dot q) 2)) (* 1/2 k (^ q 2)))
+	  :coords (q)
+	  :render (translate (* q 100) 0 (circle 20 0 0))
+	  :start (q 1.0 p_q 0.0)
+	  :params (m 1.0 k 2.0))
+
+To then run the simulation:
+
+	(make-instance 'simple-harmonic-oscillator)
+
+
+								SETUP
+
+If you have :sdl2 available and the sources are in ASDF's registry,
+loading the system should just work™:
+
+	(asdf:load-system :ham)
+
+Otherwise, the full process fetching :sdl2 with Quicklisp and adding
+the current directory to ASDF's registry is
+
+	(ql:quickload :sdl2)
+	(push *default-pathname-defaults* asdf:*central-registry*)
+	(asdf:load-system :ham)
+
+As it stands the system will only work with SBCL, though I don't think
+it would be *that* much work to get it working with other Common Lisp
+implementations.
+
+Maxima ended up being a bit of a pain to build fully within ASDF so
+the ASDF config first compiles it with its autotools set-up (so you'll
+need autoconf et al) and uses the FASLs from that build when possible.
+It's a little cursed but it does the job.  I did have to call
+asdf::mark-operation-done directly so it might break with future
+versions of ASDF (I'm using version 3.3.1 at time of writing).

drawing.lisp

@@ -92,6 +92,12 @@
     (sdl2:destroy-window window)
     (sdl2:quit)))
 
+(defmacro with-graphics-context ((var &rest args) &body body)
+  `(let ((,var (make-instance 'graphics-context ,@args)))
+     (unwind-protect
+		  (progn ,@body)
+       (cleanup ,var))))
+
 (define-transform project-to-screen (screen-height)
                   :x-transform (round x)
                   :y-transform (round (- screen-height y)))

ham.asd

@@ -3,11 +3,71 @@
 
 (in-package :asdf-user)
 
+(eval-when (:compile-toplevel :load-toplevel :execute)
+  (let* ((base-dir
+           (make-pathname :directory (pathname-directory *load-pathname*)))
+         (maxima-dir (merge-pathnames "maxima/" base-dir))
+         (maxima-src-dir (merge-pathnames "src/" maxima-dir))
+         (maxima-install-dir (merge-pathnames "maxima-install/" base-dir)))
+    (unless (probe-file (merge-pathnames "bin/maxima" maxima-install-dir))
+      (uiop:with-current-directory (maxima-dir)
+        (uiop:run-program "./bootstrap"
+                          :output *standard-output*
+                          :error-output *error-output*)
+        (uiop:run-program (format nil "./configure --prefix=~A"
+                                  (namestring maxima-install-dir))
+                          :output *standard-output*
+                          :error-output *error-output*)
+        (uiop:run-program "make -j"
+                          :output *standard-output*
+                          :error-output *error-output*)
+        (uiop:run-program "make install"
+                          :output *standard-output*
+                          :error-output *error-output*)))
+    (pushnew maxima-src-dir
+             asdf:*central-registry*
+             :test #'equal)
+    (pushnew (merge-pathnames "maxima-interface/" base-dir)
+             asdf:*central-registry*
+             :test #'equal)
+    (defparameter *maxima-src-dir* maxima-src-dir)))
+
+(defun prebuilt-maxima-fasl (source-file)
+  (let ((source-path (asdf:component-pathname source-file)))
+    (when (uiop:subpathp source-path *maxima-src-dir*)
+      (let ((fasl-file
+              (make-pathname
+               :directory (substitute "src/binary-sbcl" "src"
+                                      (pathname-directory source-path)
+                                      :test #'string=)
+               :name (pathname-name source-path)
+               :type "fasl")))
+        (when (probe-file fasl-file)
+          (list fasl-file))))))
+
+(defmethod asdf:output-files :around ((op asdf:compile-op)
+                                      (c asdf:cl-source-file))
+  (or (prebuilt-maxima-fasl c)
+      (call-next-method)))
+
+(defmethod asdf:perform :around ((op asdf:compile-op)
+                                 (c asdf:cl-source-file))
+  (if (prebuilt-maxima-fasl c)
+      (asdf::mark-operation-done op c)
+      (call-next-method)))
+
+(defmethod asdf:input-files :around ((op asdf:load-op)
+                                     (c asdf:cl-source-file))
+  (or (prebuilt-maxima-fasl c)
+      (call-next-method)))
+
 (defsystem :ham
   :version "0.1.0"
   :author "Camden Dixie O'Brien"
   :license "AGPL-3.0-only"
   :description "A framework for creating physics simulations"
-  :depends-on (:sdl2)
+  :depends-on (:sdl2 :maxima :maxima-interface)
   :components ((:file "package")
-			   (:file "drawing" :depends-on ("package"))))
+               (:file "drawing" :depends-on ("package"))
+               (:file "physics-compiler" :depends-on ("package"))
+               (:file "simulation" :depends-on ("package" "drawing"))))

ham.mac

@@ -9,7 +9,7 @@ ham(L, qs) := block([ps, qdots, H],
 	H: subst(solve(maplist(lambda([p,qdot], p = diff(L, qdot)), ps, qdots),
 	               qdots),
 	         apply("+", maplist("*", ps, qdots)) - L),
-	maplist(trigsimp, flatten(maplist(lambda([q], [
-		'diff(q,t) = diff(H, concat(p_, q)),
-		'diff(concat(p_, q),t) = -diff(H, q)
-	]), qs))));
+	maplist(trigsimp, maplist(lambda([q], [
+		diff(H, concat(p_, q)),
+        -diff(H, q)
+	]), qs)));

package.lisp

@@ -7,4 +7,4 @@
   (:use :common-lisp :sdl2)
   (:export :line :circle :rectangle :centre-rectangle
 		   :define-transform :translate :scale :rotate
-		   :overlay :graphics-context :cleanup :display-frame))
+		   :overlay :simulation :define-simulation))

physics-compiler.lisp

@@ -0,0 +1,82 @@
+;; Copyright (c) Camden Dixie O'Brien
+;; SPDX-License-Identifier: AGPL-3.0-only
+
+(in-package :ham)
+
+(eval-when (:load-toplevel :execute)
+  (maxima-interface:maxima-eval
+   `((maxima::$load)
+     ,(namestring (asdf:system-relative-pathname :ham "ham.mac")))))
+
+(defun enmaxima-symbol (sym)
+  (intern (format nil "$~a" (symbol-name sym))
+		  (find-package :maxima)))
+
+(defun enmaxima-lagrangian (expr)
+  (defun prepare (expr)
+	(cond
+	  ((symbolp expr) (enmaxima-symbol expr))
+	  ((and (listp expr) (eq 'dot (car expr)))
+	   (append '(maxima::%derivative)
+			   (mapcar #'prepare (cdr expr))
+			   '(maxima::$t)))
+	  ((and (listp expr) (eq '^ (car expr)))
+	   (cons 'expt
+			 (mapcar #'prepare (cdr expr))))
+	  ((listp expr)
+	   (cons (car expr)
+			 (mapcar #'prepare (cdr expr))))
+	  (t expr)))
+  (maxima-interface::lisp-to-maxima (prepare expr)))
+
+(defun enmaxima-coords (coords)
+  (cons '(maxima::mlist) (mapcar #'enmaxima-symbol coords)))
+
+(defun hamilton-eqns (lagrangian coords)
+  (maxima-interface:maxima-eval
+   `((maxima::$ham)
+	 ,(enmaxima-lagrangian lagrangian)
+	 ,(enmaxima-coords coords))))
+
+(defun unmaxima-expression (expr)
+  (defun unmaxima-symbol (sym)
+	(let ((name (symbol-name sym)))
+	  (if (and (eq (symbol-package sym) (find-package :maxima))
+			   (> (length name) 1)
+			   (char= (char name 0) #\$))
+		  (intern (subseq name 1)) sym)))
+  (cond
+	((symbolp expr) (unmaxima-symbol expr))
+	((listp expr) (mapcar #'unmaxima-expression expr))
+	(t expr)))
+
+(defun momentum-symbol (q)
+  (intern (format nil "P_~a" (symbol-name q))))
+
+(defun next-symbol (s)
+  (intern (format nil "NEXT-~a" (symbol-name s))))
+
+(defun update-logic (result coords)
+  (defun coord-update (q exprs)
+	`(,(next-symbol q)
+	  (+ ,q (* dt ,(unmaxima-expression (cadr exprs))))))
+  (defun momentum-update (q exprs)
+	(let ((p (momentum-symbol q)))
+	  `(,(next-symbol p)
+		(+ ,p (* dt ,(unmaxima-expression (caddr exprs)))))))
+  (loop for exprs in (cdr (maxima-interface::maxima-to-lisp result))
+        for q in coords
+		collect (coord-update q exprs)
+		collect (momentum-update q exprs)))
+
+(defun keys (plist)
+  (loop for key in plist by #'cddr collect key))
+
+(defun update-body (result coords params)
+  (let ((vars (append coords (mapcar #'momentum-symbol coords))))
+    `(let (,@(mapcar (lambda (s) `(,s (getf params ',s))) (keys params))
+           ,@(mapcar (lambda (s) `(,s (getf state ',s))) vars))
+	   (let ,(update-logic result coords)
+		 (list ,@(loop for s in vars
+					   collect `',s
+					   collect (next-symbol s)))))))

simulation.lisp

@@ -0,0 +1,54 @@
+;; Copyright (c) Camden Dixie O'Brien
+;; SPDX-License-Identifier: AGPL-3.0-only
+
+(in-package :ham)
+
+(defclass simulation ()
+  ((state :initarg :start)
+   (params :initarg :params)
+   (width :initarg :width :initform 800)
+   (height :initarg :height :initform 600)
+   (running :initform nil)
+   simulation-thread))
+
+(defgeneric update (simulation dt))
+(defgeneric render (simulation))
+
+(defmethod initialize-instance :after ((sim simulation) &key)
+  (start sim))
+
+(defmethod run ((sim simulation))
+  (sdl2:make-this-thread-main
+   (lambda ()
+     (with-slots (width height running) sim
+       (with-graphics-context (ctx :width width :height height)
+         (loop while running do
+           (update sim (/ 1.0 60.0))    ; assume 60 Hz
+           (let ((frame (render sim)))
+             (display-frame ctx frame))))))))
+
+(defmethod start ((sim simulation))
+  (with-slots (running simulation-thread) sim
+    (setf running t)
+    (setf simulation-thread
+          (sb-thread:make-thread (lambda () (run sim))
+                                 :name "ham-simulation-thread"))))
+
+(defmethod stop ((sim simulation))
+  (with-slots (running simulation-thread) sim
+    (setf running nil)
+    (sb-thread:join-thread simulation-thread)))
+
+(defmacro define-simulation (name &key lagrangian coords render start params)
+  `(progn
+     (defclass ,name (simulation)
+       ()
+       (:default-initargs :start ',start :params ',params))
+     (defmethod update ((sim ,name) dt)
+       (with-slots (state params) sim
+         (setf state ,(let ((eqns (hamilton-eqns lagrangian coords)))
+                        (update-body eqns coords params)))))
+     (defmethod render ((sim ,name))
+       (with-slots (state) sim
+         (let ,(mapcar (lambda (q) `(,q (getf state ',q))) coords)
+           ,render)))))