--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/intensity.lisp	Sun Oct 13 14:00:24 2019 +0200
@@ -0,0 +1,34 @@
+#!/usr/bin/env -S sbcl --script
+;; Calculate intensity of polarized light
+
+(defun norm (v) (sqrt (reduce (lambda (a e) (+ a (expt e 2))) v :initial-value 0)))
+
+(defun make-ray (x y) (vector x y))
+
+(defun deg-to-rad (deg) (* PI (/ deg 180)))
+
+(defun random-normed-ray (&optional (max 100))
+  (flet ((random-component () (/ (1+ (random (1- max))) max)))
+    (let* ((a (random-component))
+	   (b (random-component))
+	   (n (norm (list a b))))
+      (vector (/ a n) (/ b n)))))
+
+(defun polarize-ray (ray angle)
+  (vector (* (cos angle) (aref ray 0)) (* (sin angle) (aref ray 1))))
+
+(defun polarize-rays (rays angle)
+  (map 'vector #'(lambda (r) (polarize-ray r angle)) rays))
+
+(defun many-normed-rays (n)
+  (loop for i from 1 to n
+     collect (random-normed-ray) into rays
+     finally (return (make-array (length rays) :initial-contents rays))))
+
+(defun average-intensity (rays)
+  (/ (reduce #'+ (map 'vector #'(lambda (r) (expt (norm r) 2)) rays) :initial-value 0) (length rays)))
+
+(let* ((angle 45)
+       (sample 1000)
+       (result (average-intensity (polarize-rays (many-normed-rays sample) angle))))
+  (format t "The intensity after polarizing with 45 degree filter is: ~a~%" result))

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/macroplay.lisp	Sun Oct 13 14:00:24 2019 +0200
@@ -0,0 +1,57 @@
+(defun nth-fib (n)
+  (do ((c 0 (+ c 1))
+       (last 0 current)
+       (current 1 (+ last current)))
+      ((= n c) current)))
+
+;; Allow defining dynamic unique symbols in other macros
+(defmacro with-gensyms (symbols &body body)
+  `(let ,(loop for sym in symbols collect `(,sym (gensym)))
+     ,@body))
+
+(defmacro do-something-times (times &body body)
+  (with-gensyms (loop-var-name)
+    `(loop for ,loop-var-name from 1 to ,times do ,@body)))
+
+(defmacro once-only ((&rest names) &body body)
+  (let ((gensyms (loop for n in names collect (gensym))))
+    `(let (,@(loop for g in gensyms collect `(,g (gensym))))
+       `(let (,,@(loop for g in gensyms for n in names collect ``(,,g ,,n)))
+	  ,(let (,@(loop for n in names for g in gensyms collect `(,n ,g)))
+	     ,@body)))))
+
+(defmacro do-something-times-primitive (times &body body)
+  (oo (times)
+    `(loop for loop-var from 1 to ,times do ,@body)))
+
+(defmacro oo ((&rest names) &body body)
+  ;; For each name, generate a symbol that is used
+  ;; in the invoking macro to store the unique symbol;
+  ;; for the emitted code
+  (let ((gensyms (loop for n in names collect (gensym))))
+    ;; in other defmacro (level-1)
+    ;; generate symbols to store evaluated names in. Assign them
+    ;; to unique (l-0) symbols generated above.
+    `(let (,@(loop for g in gensyms collect `(,g (gensym))))
+       ;; in final code (level-2)
+       ;; evaluate names. Assign them to the locally-generated (l-1) symbols.
+       `(let (,,@(loop for g in gensyms for n in names collect ``(,,g ,,n)))
+	  ;; in other defmacro (level-1 again):
+	  ;; Bind evaluated names back to names
+	  ,(let (,@(loop for g in gensyms for n in names collect `(,n ,g)))
+	     ,@body)))))
+
+;; stuff
+(defmacro my-list (&rest vals)
+  (if (not vals) nil
+      ``(,,(car vals) . ,(my-list ,@(cdr vals)))))
+
+;; Anaphoric macros??
+
+(defmacro my-if (test if-so if-else)
+  `(let ((it ,test))
+     (if it ,if-so ,if-else)))
+
+(defmacro my-lambda (arg-list &body body)
+  `(labels ((self ,arg-list ,@body))
+     #'self))

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test3.lisp	Sun Oct 13 14:00:24 2019 +0200
@@ -0,0 +1,111 @@
+;; syntax test
+
+(load "~/dev/lisp/lists.lisp")
+(load "~/dev/lisp/macroplay.lisp")
+
+;; heap?
+
+(defun make-heap (&optional (max-size 64))
+  (list 0 (make-array max-size :initial-element nil)))
+
+(defun heap-get-pos (h pos) (aref (cadr h) pos))
+(defun heap-set-pos (h pos el) (setf (aref (cadr h) pos) el))
+
+(defun swap-heap-elems (heap ix1 ix2)
+  (let ((tmp (heap-get-pos heap ix1)))
+    (heap-set-pos heap ix1 (heap-get-pos heap ix2))
+    (heap-set-pos heap ix2 tmp)))
+(defun left-child (ix) (+ 1 (* 2 ix)))
+(defun right-child (ix) (+ 2 (* 2 ix)))
+(defun parent (ix) (floor (/ ix 2)))
+
+(defun bubble-heap-up (heap pos)
+  "Take the element at pos and move elements as long as the invariant is not fulfilled."
+  (loop while (and (> pos 0) (< (heap-get-pos heap (parent pos)) (heap-get-pos heap pos)))
+     do
+       (swap-heap-elems heap (parent pos) pos)
+       (setf pos (parent pos))))
+
+(defun bubble-heap-fill (heap)
+  "After popping the top element, move elements until the invariant is fulfilled"
+  (let ((pos 0))
+	(loop while (heap-get-pos heap pos) do
+	     (let ((left-elem (heap-get-pos heap (left-child pos)))
+		   (right-elem (heap-get-pos heap (right-child pos))))
+	       (cond ((and (eql nil left-elem) (not (eql nil right-elem)))
+		      (progn
+			(swap-heap-elems heap (right-child pos) pos)
+			(setf pos (right-child pos))))
+		     ((and (eql nil right-elem) (not (eql nil left-elem)))
+		      (progn
+			(swap-heap-elems heap (left-child pos) pos)
+			(setf pos (left-child pos))))
+		     ((and (eql nil left-elem) (eql nil right-elem)) (return))
+		     ('t
+		      (if (> left-elem right-elem)
+			  (progn
+			    (swap-heap-elems heap (left-child pos) pos)
+			    (setf pos (left-child pos)))
+			  (progn
+			    (swap-heap-elems heap (right-child pos) pos)
+			    (setf pos (right-child pos))))))))
+	(heap-set-pos heap pos nil)))
+
+(defvar *debug-heap-level* 0)
+(defun debug-heap-at-level (heap target-level pos)
+  (defun or_ (a b) (or a b))
+  (if (/= target-level *debug-heap-level*)
+      (let ((*debug-heap-level* (1+ *debug-heap-level*)))
+	(or_ (debug-heap-at-level heap target-level (left-child pos))
+	     (debug-heap-at-level heap target-level (right-child pos))))
+      (if (< pos (length heap))
+	  (progn (format t "~a " (heap-get-pos heap pos)) 't)
+	  nil)))
+
+(defun debug-heap (heap)
+  (loop
+     for level from 0
+     with ok = 't
+     while ok
+     do
+       (format t "level ~a:~%" level)
+       (setf ok (debug-heap-at-level heap level 0))
+       (format t "~%")))
+
+(defun put-heap (heap elem)
+  (let ((elems (if (listp elem) elem (list elem)))
+	(h (cadr heap)))
+    (loop for e in elems do
+	 (let ((newpos (car heap)))
+	   (heap-set-pos heap newpos e)
+	   (incf (car heap))
+	   (bubble-heap-up heap newpos)))))
+
+(defun pop-heap (heap)
+  (let ((elem (heap-get-pos heap 0)))
+    (if (eql 0 (car heap))
+	nil
+	(progn (decf (car heap))
+	       (bubble-heap-fill heap)
+	       elem))))
+
+(defun len-heap (heap)
+  (car heap))
+
+(defparameter *test-suite* 'heap-self-test)
+(defmacro check-true (form)
+  (let ((result-name (gensym)))
+    `(let ((,result-name ,form))
+       (format t "~a: ~:[FAIL~;pass~] ~a ~%" *test-suite* ,result-name ',form)
+       ,result-name)))
+
+(defun heap-self-test ()
+  (let ((heap (make-heap)))
+    (or (and (check-true (eql nil (put-heap heap '(5 2 1 8 3 0))))
+	     (check-true (eql 8 (pop-heap heap)))
+	     (check-true (eql 5 (pop-heap heap)))
+	     (check-true (eql 4 (len-heap heap))))
+	(progn (format t "~a~%" heap) (debug-heap heap)))
+    ))
+
+(heap-self-test)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testfw.lisp	Sun Oct 13 14:00:24 2019 +0200
@@ -0,0 +1,19 @@
+
+(defun report-result (name result form)
+  (format t "~a : ~:[FAIL~;pass~] ~a~%" name result form)
+  result)
+
+(defmacro strict-and (&rest cases)
+  `(let ((evaluated (list ,@cases)))
+     (and evaluated)))
+
+(defmacro run-tests (name &rest cases)
+  `(strict-and ,@(loop for case in cases collect
+	     `(report-result ,name ,case (quote ,case)))))
+
+(format t "~%")
+(run-tests '+-tests
+	   (eq 2 (+ 1 2))
+	   (eq 5 (+ 2 3))
+	   )
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/threebody.lisp	Sun Oct 13 14:00:24 2019 +0200
@@ -0,0 +1,110 @@
+;; Calculate 3-body problem (later: generalize to n-body)
+
+;; General helpers
+
+(defun unique-pairs (up-to)
+  "Generate all unique combinations of numbers from 0 to up-to - 1 without (n,n) pairs. I.e., (1,2), (1, 3), ..., (2, 3), (2, 4), ..."
+  (let ((result '()))
+    (loop for i from 0 to (1- up-to) do
+	 (loop for j from (1+ i) to (1- up-to) do
+	      (push (list i j) result)))
+    result))
+
+;; VEC impl
+
+(defstruct vec
+  (x 0 :type number)
+  (y 0 :type number)
+  (z 0 :type number))
+
+(defmacro for-all-vec-components ((vec var) &body body)
+  `(progn ,@(loop for C in (list 'vec-x 'vec-y 'vec-z) collect
+       `(let ((,var (,C ,vec))) ,@body))))
+
+;(for-all-vec-components ((body-pos *tbody1*) comp) (format t "component: ~a~%" comp))
+
+(defun norm-of (vec)
+  (sqrt (+ (expt (vec-x vec) 2) (expt (vec-y vec) 2) (expt (vec-z vec) 2))))
+
+(defun real-product (vec real)
+  (make-vec :x (* real (vec-x vec))
+	    :y (* real (vec-y vec))
+	    :z (* real (vec-z vec))))
+
+(defun vec-sum (&rest vecs)
+  (flet ((sum (list) (reduce #'+ list)))
+    (make-vec :x (sum (mapcar #'vec-x vecs))
+	      :y (sum (mapcar #'vec-y vecs))
+	      :z (sum (mapcar #'vec-z vecs)))))
+
+;; BODY impl
+
+(defstruct body
+  (mass 0 :type number)
+  (pos (make-vec) :type vec)
+  (vel (make-vec) :type vec))
+
+(defparameter *tbody1*
+  (make-body :mass 20000000000
+	     :pos (make-vec :x 0l0 :y 0l0 :z 0l0)
+	     :vel (make-vec :x 0.005 :y 0 :z 0)))
+(defparameter *tbody3*
+  (make-body :mass 10000000
+	     :pos (make-vec :x 2500l0 :y 3000l0 :z -200l0)
+	     :vel (make-vec :x 0.01 :y 0 :z 0)))
+(defparameter *tbody2*
+  (make-body :mass 10000000
+	     :pos (make-vec :x -2500l0 :y -3000l0 :z 500l0)
+	     :vel (make-vec :x -0.01 :y 0 :z 0)))
+
+(defconstant GRAVCONST 6.67428l-11)
+
+(defun direction-between (body1 body2)
+  (let* ((p1 (body-pos body1))
+	(p2 (body-pos body2))
+	(dx (- (vec-x p2) (vec-x p1)))
+	(dy (- (vec-y p2) (vec-y p1)))
+	(dz (- (vec-z p2) (vec-z p1))))
+    (make-vec :x dx :y dy :z dz)))
+
+(defun distance-between (body1 body2) (norm-of (direction-between body1 body2)))
+
+(defun gravitational-force (body1 body2)
+  (let* ((direction (direction-between body1 body2))
+	 (dist (norm-of direction))
+	 (unity-direction (real-product direction (/ 1 dist)))
+	 (force (/ (* GRAVCONST (body-mass body1) (body-mass body2))
+		   (expt dist 2)))
+	 (directed-force (real-product unity-direction force)))
+	 directed-force))
+
+(defun accelerate-body (body force delta-t)
+  (let* ((accel (real-product force (/ 1 (body-mass body))))
+	 (delta-s-by-vel (real-product (body-vel body) delta-t))
+	 (delta-s-by-accel (real-product accel (* 1/2 (expt delta-t 2))))
+	 (delta-v-by-accel (real-product accel delta-t)))
+    (setf (body-pos body) (vec-sum (body-pos body) delta-s-by-vel delta-s-by-accel))
+    (setf (body-vel body) (vec-sum (body-vel body) delta-v-by-accel))
+    body))
+
+(defun calculate-body-step (bodies-vec delta-t)
+  (loop for pair in (unique-pairs (length bodies-vec)) do
+       (let* ((b1 (first pair))
+	      (b2 (second pair))
+	      (force-2-on-1 (gravitational-force (aref bodies-vec b1) (aref bodies-vec b2)))
+	      (force-1-on-2 (real-product force-2-on-1 -1)))
+	 (accelerate-body (aref bodies-vec b1) force-2-on-1 delta-t)
+	 (accelerate-body (aref bodies-vec b2) force-1-on-2 delta-t))))
+
+(defparameter *bodies-state-output* t)
+(defun output-bodies-state (bodies-vec)
+  (flet ((body-state (body)
+	   (format *bodies-state-output* "~,3f, ~,3f, ~,3f"
+		   (vec-x (body-pos body))
+		   (vec-y (body-pos body))
+		   (vec-z (body-pos body)))))
+    (loop for body across bodies-vec do
+	 (body-state body)
+	 (format *bodies-state-output* ", "))
+    (format *bodies-state-output* "~%")))
+