SVG, STL file exporter [draft]

src/emmy/numerical/ode.cljc

@@ -278,11 +278,10 @@
   (make-integrator state-derivative state-derivative-args))
 
 (defn integrate-state-derivative
-  "A wrapper for evolve, which is more convenient when you just
-  want a vector of (time, state) pairs over the integration interval
-  instead of having to deal with a callback. Integrates the supplied
-  state derivative (and its argument package) from [0 to t1] in steps
-  of size dt"
+  "A wrapper for evolve, which is more convenient when you just want a vector of
+  states over the integration interval instead of having to deal with a
+  callback. Integrates the supplied state derivative (and its argument package)
+  from [0 to t1] in steps of size dt"
   [state-derivative state-derivative-args initial-state t1 dt]
   (let [f (make-integrator* state-derivative state-derivative-args initial-state
                             {:epsilon 1e-6

src/emmy/stl.clj

@@ -0,0 +1,104 @@
+(ns emmy.stl
+  (:require [clojure.java.io :as io])
+  (:import (java.nio ByteBuffer ByteOrder)))
+
+(defn- float->bytes [f]
+  (let [bb (ByteBuffer/allocate 4)]
+    (.order bb ByteOrder/LITTLE_ENDIAN)
+    (.putFloat bb f)
+    (.array bb)))
+
+(defn- vec3->bytes [[x y z]]
+  (byte-array (concat (float->bytes (float x))
+                      (float->bytes (float y))
+                      (float->bytes (float z)))))
+
+(defn- normal-vector [[x1 y1 z1] [x2 y2 z2] [x3 y3 z3]]
+  (let [ux (- x2 x1)
+        uy (- y2 y1)
+        uz (- z2 z1)
+        vx (- x3 x1)
+        vy (- y3 y1)
+        vz (- z3 z1)
+        nx (- (* uy vz) (* uz vy))
+        ny (- (* uz vx) (* ux vz))
+        nz (- (* ux vy) (* uy vx))
+        magnitude (Math/sqrt (+ (* nx nx) (* ny ny) (* nz nz)))]
+    [(/ nx magnitude) (/ ny magnitude) (/ nz magnitude)]))
+
+(defn- triangle->bytes [v1 v2 v3]
+  (let [normal (normal-vector v1 v2 v3)]
+    (byte-array (concat (vec3->bytes normal)
+                        (vec3->bytes v1)
+                        (vec3->bytes v2)
+                        (vec3->bytes v3)
+                        [0 0]))))  ; Attribute byte count (unused)
+
+(defn- generate-triangles [points rows cols]
+  (for [i (range rows)
+        j (range cols)
+        :let [next-i (mod (inc i) rows)
+              next-j (mod (inc j) cols)
+              v1 (nth points (+ (* i cols) j))
+              v2 (nth points (+ (* i cols) next-j))
+              v3 (nth points (+ (* next-i cols) j))
+              v4 (nth points (+ (* next-i cols) next-j))]]
+    [(triangle->bytes v1 v2 v3)
+     (triangle->bytes v2 v4 v3)]))
+
+(defn parametric-to-stl
+  "Convert a parametric function to an STL file"
+  [parametric-fn filename & {:keys [u-start u-end u-steps
+                                    v-start v-end v-steps]
+                             :or {u-start 0 u-end (* 2 Math/PI) u-steps 100
+                                  v-start 0 v-end Math/PI v-steps 50}}]
+  (let [u-range (range u-start u-end (/ (- u-end u-start) u-steps))
+        v-range (range v-start v-end (/ (- v-end v-start) v-steps))
+        points (vec (for [v v-range
+                          u u-range]
+                      (parametric-fn u v)))
+        triangles (generate-triangles points v-steps u-steps)
+        triangle-count (* 2 u-steps v-steps)]
+    (with-open [out (io/output-stream filename)]
+      ;; Write STL header
+      (.write out (byte-array 80))
+      ;; Write triangle count
+      (.write out (.array (-> (.order (ByteBuffer/allocate 4) ByteOrder/LITTLE_ENDIAN)
+                              (.putInt triangle-count))))
+      ;; Write triangles
+      (doseq [triangle (flatten triangles)]
+        (.write out triangle)))
+    filename))
+
+;; Example usage:
+(comment
+  (require '[emmy.env :refer :all])
+
+  (defn toroidal->rect [R r]
+    (fn [theta phi]
+      (*
+       (rotate-z-matrix phi)
+       (up (+ R (* r (cos theta)))
+           0
+           (* r (sin theta))))))
+
+  (parametric-to-stl (toroidal->rect 50 10) "torus.stl"
+                     :u-start 0
+                     :u-end (* 2 Math/PI)
+                     :v-start 0
+                     :v-end (* 2 Math/PI)
+                     :u-steps 200
+                     :v-steps 100)
+
+  (defn example-parametric-fn [u v]
+    [(* 50 (Math/cos u) (Math/sin v))
+     (* 50 (Math/sin u) (Math/sin v))
+     (* 50 (Math/cos v))])
+
+  (parametric-to-stl example-parametric-fn "output.stl"
+                     :u-start 0
+                     :u-end (* 2 Math/PI)
+                     :v-start 0
+                     :v-end Math/PI
+                     :u-steps 20
+                     :v-steps 20))

src/emmy/svg.clj

@@ -0,0 +1,83 @@
+(ns emmy.svg
+  (:require [clojure.java.io :as io]
+            [clojure.string :as str]
+            [emmy.env :as e :refer :all]))
+
+(defn scale-points
+  "Scale points to fit within the SVG viewbox"
+  [points width height]
+  (let [x-values (map first points)
+        y-values (map second points)
+        x-min (apply min x-values)
+        x-max (apply max x-values)
+        y-min (apply min y-values)
+        y-max (apply max y-values)
+        x-range (- x-max x-min)
+        y-range (- y-max y-min)
+        scale (min (/ width x-range) (/ height y-range))]
+    (for [[x y] points]
+      [(* scale (- x x-min))
+       (* scale (- y y-min))])))
+
+(defn points-to-path
+  "Convert a sequence of points to an SVG path string"
+  [points]
+  (str "M " (str/join " L " (map #(str/join "," %) points))))
+
+(defn generate-svg
+  "Generate an SVG string from a sequence of points"
+  [points width height]
+  (let [scaled-points (scale-points points width height)
+        path-data (points-to-path scaled-points)]
+    (format
+     "<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 %d %d\">
+        <path d=\"%s\" fill=\"none\" stroke=\"black\" stroke-width=\"1\"/>
+      </svg>"
+     width height path-data)))
+
+(defn coordinates-to-svg
+  "Convert a vector of [x y] coordinates to an SVG file"
+  [coordinates filename & {:keys [width height]
+                           :or {width 1000 height 1000}}]
+  (let [svg-content (generate-svg coordinates width height)]
+    (with-open [writer (io/writer filename)]
+      (.write writer svg-content))
+    filename))
+
+(defn lorenz-attractor
+  "Generate the derivatives (dx, dy, dz) for the Lorenz attractor given the current point (x, y, z) and parameters"
+  [sigma beta rho]
+  (fn [[x y z]]
+    (let [dx (* sigma (- y x))
+          dy (- (* x (- rho z)) y)
+          dz (- (* x y) (* beta z))]
+      [dx dy dz])))
+
+;; Example usage:
+(comment
+  (def example-coordinates
+    [[0 0] [100 100] [200 50] [300 200] [400 0]])
+
+  (coordinates-to-svg example-coordinates "output.svg")
+
+  (let [sigma 10
+        beta 8/3
+        rho 28
+        t1 20
+        dt 0.005 
+        init-state [-8.0 8.0 27.0]]
+    (-> (for [[x y] (integrate-state-derivative lorenz-attractor [sigma beta rho] init-state t1 dt)]
+          [x y])
+        (coordinates-to-svg "lorenz.svg"))))
+
+
+
+(defn project-lorenz-xy
+  "Project the Lorenz attractor onto the XY plane"
+  [sigma beta rho dt steps]
+  (let [initial-state [0.1 0 0]
+        integrator (state-advancer lorenz-attractor sigma beta rho)]
+    (mapv (fn [t]
+            (let [[x y _] (integrator initial-state t)]
+              [x y]))
+          (range 0 (* dt steps) dt))))