geodesics and more physics support (#63)

CHANGELOG.md

@@ -2,6 +2,19 @@
 
 ## [unreleased]
 
+- #63:
+
+  - Adds
+    `emmy.mathbox.physics.{lagrangian-curve,hamiltonian-curve,routhian-curve,geodesic}`
+    to the existing `ode-curve` function.
+
+  - Adds
+    `emmy.viewer.physics.{evolve-lagrangian,evolve-hamiltonian,evolve-routhian}`
+    to the existing `evolve` function.
+
+  - Updates the `emmy.viewer.components.physics/Evolve` component to use a
+    pre-allocated JS output array instead of forcing an allocation on each tick.
+
 ## [0.2.0]
 
 - #57:

dev/examples/manifold/klein.clj

@@ -5,6 +5,7 @@
    :exclude [+ - * / zero? compare divide numerator denominator
              infinite? abs ref partial =])
   (:require [emmy.clerk :as ec]
+            [emmy.mathbox.physics :as ph]
             [emmy.env :as e :refer :all]
             [emmy.leva :as leva]
             [emmy.mathbox.plot :as p]
@@ -17,17 +18,10 @@
 
 ;; # Fun with Klein Bottles
 
-(defn slider-surface [name {:keys [u v] :as opts}]
+(defn slider-surface [name {:keys [u v] :as opts} & children]
   (ev/with-let [!opts {:u (peek u) :v (peek v)}]
-    (p/scene
-     {:axes
-      {:x {:divisions 10
-           :label {:position 3}}
-       :y {:divisions 10
-           :ticks {:labels? false}}
-       :z {:divisions 10
-           :label-ticks? false}}
-      :grids []}
+    (apply
+     p/scene
      (leva/controls
       {:folder {:name name}
        :schema
@@ -37,7 +31,8 @@
      (p/parametric-surface
       (assoc opts
              :u [(first u) (ev/get !opts :u)]
-             :v [(first v) (ev/get !opts :v)])))))
+             :v [(first v) (ev/get !opts :v)]))
+     children)))
 
 ;; ## Mobius Strip
 
@@ -67,18 +62,31 @@
 
 ;; Can you find the Emmy logo hidden inside?
 
-(ev/with-let [!r {:r 2}]
+(ev/with-let [!r {:r 2 :theta_0 1 :alpha_0 0 :steps 20}]
   [:<>
-   (leva/controls {:atom !r
-                   :folder {:name "Klein Bagel"}
-                   :schema {:r {:min 2 :max 6 :step 0.01}}})
+   (leva/controls
+    {:atom !r
+     :folder {:name "Klein Bagel"}
+     :schema {:r {:min 2 :max 6 :step 0.01}
+              :theta_0 {:label (->infix 'theta_0) :min 0 :max Math/PI :step 0.02}
+              :alpha_0 {:label (->infix 'alpha_0) :min 0 :max Math/PI :step 0.02}
+              :steps   {:min 500 :max 9000 :step 50}}})
    (slider-surface
     "Klein Bagel"
     {:f (ev/with-params {:atom !r :params [:r]}
           klein-bagel)
-     :opacity 0.75
      :u [0 (* 2 Math/PI)]
-     :v [0 (* 2 Math/PI)]})])
+     :v [0 (* 2 Math/PI)]}
+
+    (ph/geodesic
+     {:x0 [(ev/get !r :theta_0) 0]
+      :v0 [(list 'Math/cos (ev/get !r :alpha_0))
+           (list 'Math/sin (ev/get !r :alpha_0))]
+      :xform (ev/with-params {:atom !r :params [:r]}
+               klein-bagel)
+      :steps (ev/get !r :steps)
+      :width 2
+      :end? true}))])
 
 ;; ## Klein bottle:
 
@@ -102,11 +110,30 @@
    (* 2/15 (sin v)
       (+ 3 (* 5 (cos u) (sin u))))])
 
-(slider-surface
- "Klein"
- {:f klein-bottle
-  :u [0 Math/PI]
-  :v [0 (* 2 Math/PI)]})
+(ev/with-let [!r {:theta_0 1 :alpha_0 0 :steps 20}]
+  [:<>
+   (leva/controls
+    {:atom !r
+     :folder {:name "Klein"}
+     :schema {:theta_0 {:label (->infix 'theta_0) :min 0 :max Math/PI :step 0.02}
+              :alpha_0 {:label (->infix 'alpha_0) :min 0 :max Math/PI :step 0.02}
+              :steps   {:min 10 :max 9000 :step 50}}})
+
+   (slider-surface
+    "Klein"
+    {:f klein-bottle
+     :opacity 0.2
+     :u [0 Math/PI]
+     :v [0 (* 2 Math/PI)]}
+
+    (ph/geodesic
+     {:x0 [(ev/get !r :theta_0) 0]
+      :v0 [(list 'Math/cos (ev/get !r :alpha_0))
+           (list 'Math/sin (ev/get !r :alpha_0))]
+      :xform klein-bottle
+      :steps (ev/get !r :steps)
+      :width 2
+      :end? true}))])
 
 ;; ## Plucker's Conoid
 

dev/examples/simulation/double_ellipsoid.clj

@@ -52,12 +52,11 @@
      {:threestrap {:plugins ["core" "controls" "cursor" "stats"]}}
      (leva/controls {:atom !opts})
 
-     (emmy.viewer.physics/evolve
+     (emmy.viewer.physics/evolve-lagrangian
       {:atom !state
        :initial-state initial-state
-       :f' (ev/with-params {:atom !opts :params [:m :k :x0 :a :b :c]}
-             (comp Lagrangian->state-derivative
-                   L-central-triaxial))})
+       :L (ev/with-params {:atom !opts :params [:m :k :x0 :a :b :c]}
+            L-central-triaxial)})
 
      (emmy.mathbox.physics/comet
       {:length 16

dev/examples/simulation/ellipsoid.clj

@@ -117,12 +117,11 @@
         :b {:min 1 :max 5 :step 0.01}
         :c {:min 1 :max 5 :step 0.01}}})
 
-     (emmy.viewer.physics/evolve
+     (emmy.viewer.physics/evolve-lagrangian
       {:atom !state
        :initial-state initial-state
-       :f' (ev/with-params {:atom !opts :params [:m :g :a :b :c]}
-             (comp Lagrangian->state-derivative
-                   L-central-triaxial))})
+       :L (ev/with-params {:atom !opts :params [:m :g :a :b :c]}
+            L-central-triaxial)})
 
      (plot/parametric-surface
       {:opacity 0.2

dev/examples/simulation/oscillator.clj

@@ -36,14 +36,14 @@
     (plot/scene
      (leva/controls {:atom !opts})
 
-     (emmy.viewer.physics/evolve
+     (emmy.viewer.physics/evolve-lagrangian
       {:atom !state
        :initial-state initial-state
-       :f' (ev/with-params {:atom !opts :params [:g :m :k]}
-             (comp e/Lagrangian->state-derivative L-harmonic))})
+       :L (ev/with-params {:atom !opts :params [:g :m :k]}
+            L-harmonic)})
 
      (emmy.mathbox.physics/comet
-      {:length        1
+      {:length        10
        :state->xyz    coordinate
        :initial-state initial-state
        :atom          !state}))))

dev/user.clj

@@ -23,7 +23,7 @@
    :browse? true
    :watch-paths ["src" "dev"]
    ;; Enable this when working on new components.
-   ;; :cljs-namespaces '[emmy-viewers.sci-extensions]
+   :cljs-namespaces '[emmy-viewers.sci-extensions]
    })
 
 (def static-defaults

src/emmy/mathbox/physics.cljc

@@ -4,6 +4,9 @@
   namespace."
   (:refer-clojure :exclude [vector])
   (:require [emmy.mathbox.plot :as plot]
+            [emmy.mechanics.lagrange :as l]
+            [emmy.mechanics.hamilton :as h]
+            [emmy.mechanics.routhian :as r]
             [emmy.viewer :as ev]
             [emmy.viewer.compile :as vc]
             [emmy.viewer.physics :as ph]))
@@ -26,6 +29,9 @@
     responsible for transforming each integrated state into a 3D point. Defaults
     to `identity`.
 
+  - `:simplify?` if `true`, the compiler will attempt to simplify the body of
+    `:state->xyz` and `:f'`. `false` by default.
+
   - `:steps`: the number of `:dt`-spaced steps for the integrator to take.
     Defaults to 1000.
 
@@ -63,6 +69,109 @@
                  ['emmy.mathbox.components.physics/ODECurve opts])
         (ev/fragment plot/scene))))
 
+(defn lagrangian-curve
+  "Returns a fragment that plots a curve by integrating a system of ordinary
+  differential equations generated from the Lagrangian passed via `:L` forward
+  from the initial input state `initial-state` for `steps` steps of `dt` each.
+
+  Required arguments:
+
+  - `:L'`: function of the form `(fn [[t q qdot]] <energy>)`
+
+  - `:initial-state`: the initial input vector to `:L`.
+
+  See [[ode-curve]] for a description of all supported optional options."
+  [{:keys [L] :as opts}]
+  (let [f' (if (ev/param-f? L)
+             (update L :f #(comp l/Lagrangian->state-derivative %))
+             (l/Lagrangian->state-derivative L))]
+    (ode-curve
+     (-> (dissoc opts :L)
+         (assoc :f' f')))))
+
+(defn hamiltonian-curve
+  "Returns a fragment that plots a curve by integrating a system of ordinary
+  differential equations generated from the Hamiltonian passed via `:H` forward
+  from the initial input state `initial-state` for `steps` steps of `dt` each.
+
+  Required arguments:
+
+  - `:H`: function of the form `(fn [[t q p]] <energy>)`
+
+  - `:initial-state`: the initial input vector to `:H`.
+
+  See [[ode-curve]] for a description of all supported optional options."
+  [{:keys [H] :as opts}]
+  (let [f' (if (ev/param-f? H)
+             (update H :f #(comp h/Hamiltonian->state-derivative %))
+             (h/Hamiltonian->state-derivative H))]
+    (ode-curve
+     (-> (dissoc opts :H)
+         (assoc :f' f')))))
+
+(defn routhian-curve
+  "Returns a fragment that plots a curve by integrating a system of ordinary
+  differential equations generated from the Routhian passed via `:R` forward
+  from the initial input state `initial-state` for `steps` steps of `dt` each.
+
+  Required arguments:
+
+  - `:R`: function of the form `(fn [[t q qdot-or-p]] <energy>)`
+
+  - `:initial-state`: the initial input vector to `:R`.
+
+  See [[ode-curve]] for a description of all supported optional options."
+  [{:keys [R] :as opts}]
+  (let [f' (if (ev/param-f? R)
+             (update R :f #(comp r/Routhian->state-derivative %))
+             (r/Routhian->state-derivative R))]
+    (ode-curve
+     (-> (dissoc opts :R)
+         (assoc :f' f')))))
+
+(defn- L-params [xform]
+  (fn [& params]
+    (comp (l/L-free-particle 1)
+          (l/F->C
+           (comp (apply xform params)
+                 l/coordinate)))))
+
+(defn- L-bare [xform]
+  ((L-params (fn [] xform))))
+
+(defn geodesic
+  "Returns a fragment that plots a geodesic curve over the parametric surface
+  defined by `:xform`, starting from initial position `:x0` with initial
+  velocity `:v0`. The geodesic curve is generated by integrating a system of
+  ordinary differential equations forward for `steps` steps of `dt` each.
+
+  Required arguments:
+
+  - `:xform`: function of the form `(fn [[u v]] [<x> <y> <z>])`
+
+  - `:x0`: starting position of the geodesic in `[u v]` coordinates.
+
+  - `:v0`: starting velocity of the geodesic in `[udot vdot]`.
+
+  See [[ode-curve]] for a description of all supported optional options."
+  [{:keys [xform x0 v0] :as opts}]
+  (let [L (if (ev/param-f? xform)
+            (update xform :f L-params)
+            (L-bare xform))
+        state->xyz (if (ev/param-f? xform)
+                     (update xform :f (fn [f]
+                                        (fn [& params]
+                                          (comp (apply f params)
+                                                l/coordinate))))
+                     (comp xform l/coordinate))]
+    (lagrangian-curve
+     (-> (dissoc opts :xform :x0 :v0)
+         (assoc :L L
+                :initial-state [0 x0 v0]
+                :state->xyz state->xyz)))))
+
+;; ## Visual Elements
+
 (defn comet
   "Returns a fragment that renders a point that trails its historical positions
   out behind it in a glowing tail.

src/emmy/viewer/components/physics.cljs

@@ -136,13 +136,13 @@
     :as opts}]
   ;; TODO make this in a hook and call the no-arity function on the hook.
   ;; exit.
-  (let [update (monotonic-integrator
-                f'
-                (:state (.-state !state))
-                (select-keys opts [:epsilon :max-steps]))]
+  (let [initial (->arr (:state (.-state !state)))
+        update  (monotonic-integrator
+                 f'
+                 initial
+                 (select-keys opts [:epsilon :max-steps]))]
     [sw/Stopwatch
      {:onTick
-      (fn [t]
-        ;; TODO can we keep the output here mutable and provide an out to
-        ;; update?
-        (swap! !state assoc :state (update t)))}]))
+      (let [out (js/Array. (count initial))]
+        (fn [t]
+          (swap! !state assoc :state (update t out))))}]))