1 files changed, 78 insertions, 0 deletions
diff --git a/test/damped-harmonic.scm b/test/damped-harmonic.scm
new file mode 100644
index 0000000..a47be34
--- /dev/null
+++ b/test/damped-harmonic.scm
@@ -0,0 +1,78 @@
+
+; This is the "Example" from R7RS-little
+
+(define (integrate-system system-derivative
+                          initial-state
+                          h)
+  (let ((next (runge-kutta-4 system-derivative h)))
+    (letrec ((states
+               (cons initial-state
+                     (delay (map-streams next
+                                         states)))))
+      states)))
+
+(define (runge-kutta-4 f h)
+  (let ((*h (scale-vector h))
+        (*2 (scale-vector 2))
+        (*1/2 (scale-vector (/ 1 2)))
+        (*1/6 (scale-vector (/ 1 6))))
+    (lambda (y)
+      ;; y is a system state
+      (let* ((k0 (*h (f y)))
+             (k1 (*h (f (add-vectors y (*1/2 k0)))))
+             (k2 (*h (f (add-vectors y (*1/2 k1)))))
+             (k3 (*h (f (add-vectors y k2)))))
+        (add-vectors y
+            (*1/6 (add-vectors k0
+                               (*2 k1)
+                               (*2 k2)
+                               k3)))))))
+
+(define (elementwise f)
+  (lambda vectors
+    (generate-vector
+      (vector-length (car vectors))
+      (lambda (i)
+        (apply f
+               (map (lambda (v) (vector-ref v i))
+                    vectors))))))
+
+(define (generate-vector size proc)
+  (let ((ans (make-vector size)))
+    (letrec ((loop
+               (lambda (i)
+                 (cond ((= i size) ans)
+                       (else
+                         (vector-set! ans i (proc i))
+                         (loop (+ 1 i)))))))
+      (loop 0))))
+
+(define add-vectors (elementwise +))
+
+(define (scale-vector s)
+  (elementwise (lambda (x) (* x s))))
+
+(define (map-streams f s)
+  (cons (f (head s))
+        (delay (map-streams f (tail s)))))
+
+(define head car)
+(define (tail stream)
+  (force (cdr stream)))
+
+; Ok, now the actual example
+
+(define (damped-oscillator R L C)
+  (lambda (state)
+    (let ((Vc (vector-ref state 0))
+          (Il (vector-ref state 1)))
+      (vector (- 0 (+ (/ Vc (* R C)) (/ Il C)))
+              (/ Vc L)))))
+
+(define the-states
+  (integrate-system
+    (damped-oscillator 10000 1000 .001)
+    '#(1 0)
+    .01))
+
+(display the-states)