#!/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))