/*
 * Copyright (c) Camden Dixie O'Brien
 * SPDX-License-Identifier: AGPL-3.0-only
 */

ham(L, qs) := block([ps, qdots, H],
	ps: maplist(lambda([q], concat(p_, q)), qs),
	qdots: maplist(lambda([q], 'diff(q, t)), qs),
	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))));