;;; Lindenmayer System. Implements this commands:
;;; http://www.biologie.uni-hamburg.de/b-online/e28_3/lsys.html
;;;
;;; See http://algorithmicbotany.org/ and http://www.ii.uib.no/~knute/lsystems/llisp.html
;;; for more advanced applications

(defun starts-with (string search start)
  (let ((string-len (length string))
        (search-len (length search)))
    (if (> (+ start search-len) string-len)
        nil
        (string= string search :start1 start :end1 (+ start search-len)))))

(defun search-longest (rules string start)
  "Returns the successor and predecessor length, if a predecessor
   is found in 'string' at 'start'. Returns the longest match."
  (let ((found-successor nil)
        (found-predecessor-length 0))
    (loop for (predecessor . successor) in rules do
          (when (starts-with string predecessor start)
            (let ((predecessor-length (length predecessor)))
              (when (>= predecessor-length found-predecessor-length)
                (setq found-predecessor-length predecessor-length)
                (setq found-successor successor)))))
    (values found-successor found-predecessor-length)))

(defun get-max-successor-length (rules)
  (let ((max 0))
    (loop for (predecessor . successor) in rules do
          (let ((len (length successor)))
            (when (> len max) (setq max len))))
    max))
  
(defun apply-rules (rules axiom)
  "Returns the next iteration, starting with the string in 'axiom'."
  (let ((result (make-string (max (length axiom) (* (length axiom) (get-max-successor-length rules)))))
        (result-length 0))
    (loop for i from 0 to (1- (length axiom)) do
          (multiple-value-bind (successor predecessor-length)
              (search-longest rules axiom i)
            (if successor
                (progn
                  (replace result successor :start1 result-length)
                  (incf result-length (length successor))
                  (setq i (1- (+ i predecessor-length))))
                (progn
                  (replace result (subseq axiom i (1+ i)) :start1 result-length)
                  (incf result-length)))))
    (subseq result 0 result-length)))

(defun l-system (rules axiom depth)
  "Returns 'depth' iterations, starting with 'axiom' and applying the 'rules'."
  (let ((result axiom))
    (loop repeat depth do
          (setq result (apply-rules rules result)))
    result))

(defun forward (point len angle)
  "Returns a new point by starting from 'point' and adding the polar coordinates 'len' and 'angles'."
  (let ((x (car point))
        (y (cdr point))
        (rad (* (/ (coerce pi 'single-float) 180.0) angle)))
    (cons (+ x (* (sin rad) len)) (+ y (* (cos rad) len)))))

(defun do-l-system (commands len angle fun)
  "Calls 'fun x0 y0 x1 y1' for every command in the 'commands' string."
  (let ((point-stack '())
        (angle-stack '())
        (current-point '(0e0 . 0e0))
        (current-angle 0e0))
    (loop for i from 0 to (1- (length commands)) do
          (let ((command (elt commands i)))
            (cond ((eq command #\f) ;; move only
                   (setq current-point (forward current-point len current-angle)))
                  ((eq command #\F) ;; forward with draw
                   (let ((next-point (forward current-point len current-angle)))
                     (funcall fun (car current-point) (cdr current-point) (car next-point) (cdr next-point))
                     (setq current-point next-point)))
                  ((eq command #\+)
                   (setq current-angle (+ current-angle angle)))
                  ((eq command #\-)
                   (setq current-angle (- current-angle angle)))
                  ((eq command #\[)
                   (push current-point point-stack)
                   (push current-angle angle-stack))
                  ((eq command #\])
                   (setq current-point (pop point-stack))
                   (setq current-angle (pop angle-stack))))))))

(defun postscript-l-system (rules axiom angle depth)
  "Calculates and prints a Lindenmayer System as postscript."
  (format t "500 500 scale~%")
  (format t ".1 .1 translate~%")
  (format t "0 setlinewidth~%")
  (let ((commands (l-system rules axiom depth))
        (min-x 1e30)
        (min-y 1e30)
        (max-x -1e30)
        (max-y -1e30))
    (do-l-system commands 1 angle
                 (lambda (x0 y0 x1 y1)
                   (when (< x0 min-x) (setq min-x x0))
                   (when (< y0 min-y) (setq min-y y0))
                   (when (< x1 min-x) (setq min-x x1))
                   (when (< y1 min-y) (setq min-y y1))
                   (when (> x0 max-x) (setq max-x x0))
                   (when (> y0 max-y) (setq max-y y0))
                   (when (> x1 max-x) (setq max-x x1))
                   (when (> y1 max-y) (setq max-y y1))))
    (let* ((len (float (/ (max (- max-y min-y) (- max-x min-x)))))
           (min-x (* min-x len))
           (min-y (* min-y len)))
      (do-l-system commands len angle
                   (lambda (x0 y0 x1 y1)
                     (format t "~D ~D moveto ~D ~D lineto~%" (- x0 min-x) (- y0 min-y) (- x1 min-x) (- y1 min-y))))))
  (format t "stroke~%")
  (format t "showpage~%"))

;;; dragon curve
; (postscript-l-system '(("FL"."FL+FR+") ("FR"."-FL-FR")) "FL" 90 16)

;;; snowflake
(defun test ()
  (ghostscript-l-system '(("FL"."FL+FR+") ("FR"."-FL-FR")) "FL" 90 16))

(defun ghostscript-l-system (rules axiom angle depth)
  (with-open-file (s "c:/tmp/test.ps" :direction :output :if-exists :supersede)
    (let ((*standard-output* s))
      (postscript-l-system rules axiom angle depth)))
  (ext:shell "c:/gs/gs7.05/bin/gswin32.exe -g800x800 c:/tmp/test.ps"))

;;; bush
; (ghostscript-l-system '(("F"."FF-[-F+F+F]+[+F-F-F]")) "+F" 23 5)


;;; quadratic Koch island
;(ghostscript-l-system '(("F"."F+F-F-FF+F+F-F")) "F+F+F+F" 90 3)

;;; fern
;(ghostscript-l-system '(("F"."F[-F]F[+F][F]")) "F" 30 5)

;;; Sierpinski triangle
;(ghostscript-l-system '(("F"."FF") ("X"."--FXF++FXF++FXF--")) "X" 60 7)

;;; Sierpienski carpet
;(ghostscript-l-system
; '(("F"."F+F-F-F-f+F+F+F-F")
;   ("f"."fff"))
; "F" 90 5)

;;; Hilbert curve (L -> +RF-LFL-FR+, R -> -LF+RFR+FL-):
;(ghostscript-l-system '(("L"."+RF-LFL-FR+") ("R"."-LF+RFR+FL-")) "L" 90 7)

;;; Penrose tiling
(ghostscript-l-system
 '(("F"."")
   ("W"."YF++ZF----XF[-YF----WF]++")
   ("X"."+YF--ZF[---WF--XF]+")
   ("Y"."-WF++XF[+++YF++ZF]-")
   ("Z"."--YF++++WF[+ZF++++XF]--XF"))
 "+WF--XF---YF--ZF" 36 5)

;;; hexagon tiling
;(ghostscript-l-system '(("F"."F+F+F+F+F+f+f-F-F-F-F-f-f")) "F" 60 4)

;;; square tiling
;(ghostscript-l-system '(("F"."F+F+F+F+f-F-F-F-f")) "F" 90 5)