Write a README

.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,42 @@
+				 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
+
+Before the ASDF system can be loaded, you must build the Maxima
+sources locally and install them into a directory (I couldn't get
+Maxima's ASDF system to load successfully otherwise, though I'm not
+entirely sure why this fixes it).  This can be done with:
+
+	cd maxima
+	mkdir $PWD-install
+	./bootstrap
+	./configure --prefix=$PWD-install
+	make -j
+	make install
+
+I'll try to figure out if this step itself can be done with ASDF at
+some point ;)
+
+Once that's done you should be able to load the system with
+(asdf:load-system :ham), given that the sources are in your ASDF
+registry.  The system also uses :sdl2, so make sure that's available
+-- you can get it with (ql:quickload :sdl2) if you're so inclined.

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,23 @@
 
 (in-package :asdf-user)
 
+(eval-when (:compile-toplevel :load-toplevel :execute)
+  (let ((base-dir
+		  (make-pathname :directory (pathname-directory *load-pathname*))))
+	(pushnew (merge-pathnames "maxima/src/" base-dir)
+			 asdf:*central-registry*
+			 :test #'equal)
+	(pushnew (merge-pathnames "maxima-interface/" base-dir)
+			 asdf:*central-registry*
+			 :test #'equal)))
+
 (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)));

physics-compiler.lisp

@@ -0,0 +1,81 @@
+;; 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 momentum-symbol (q)
+  (intern (format nil "P_~a" (symbol-name q))))
+
+(defun flatten (l)
+  (apply #'append l))
+
+(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 update-logic (result coords)
+  (defun coord-updates (exprs q)
+	(let ((p (momentum-symbol q))
+		  (q-dot-expr (cadr exprs))
+		  (p-dot-expr (caddr exprs)))
+	  `((incf ,p (* dt ,(unmaxima-expression p-dot-expr)))
+		(incf ,q (* dt ,(unmaxima-expression q-dot-expr))))))
+  (flatten (loop for coord-expressions
+				   in (cdr (maxima-interface::maxima-to-lisp result))
+				 for q in coords
+				 collect (coord-updates coord-expressions q))))
+
+(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))
+	   ,@(update-logic result coords)
+	   (list ,@(flatten
+				(mapcar (lambda (s) (list `',s s))
+						vars))))))

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)))))