; console program
; Lisp / Common Lisp / SBCL 
; based on : 
; http://www.mostlymaths.net/2009/08/lavaurs-algorithm.html
; lisp code by R Berenguel
; reducible angle is not angle of lower period 
; for example 9/15 = 3/5 has period 4 

;
; Rules by Lavaurs :
; no arc crosses any other arc
; arc connects 2 angles with the same period : first and second
; first angle is the smallest angles not yet connected, and second angle is the next smallest angle not yet connected
;
; orthogonal circles  (x1,y1,r1) and (x2,y2,r2)
; r1^2 + r2^2 = (x2-x1)^2 +(y2-y1)^2
; http://planetmath.org/encyclopedia/OrthogonalCircle.html
; http://classes.yale.edu/fractals/Labs/NonLinTessLab/BasicConstr3.html
; 
; example of use : 
; 
; sbcl 
; (load "lsv.lisp")
; ; look for lavaurs-5.svg file in your home directory
; ; after loading file you can :
; (draw-lavaurs "lavaurs-6.svg" 2000 6)
;
; Adam Majewski
; fraktal.republika.pl
; 2010.09.04- 11.17
;
;
;;  This program is free software: you can redistribute it and/or
;;  modify it under the terms of the GNU General Public License as
;;  published by the Free Software Foundation, either version 3 of the
;;  License, or (at your option) any later version.

;;  This program is distributed in the hope that it will be useful,
;;  but WITHOUT ANY WARRANTY; without even the implied warranty of
;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
;;  General Public License for more details.

;;  You should have received a copy of the GNU General Public License
;;  along with this program. If not, see
;;  <http://www.gnu.org/licenses/>.

(defun doubling-map (ratio-angle)
" period doubling map =  The dyadic transformation (also known as the dyadic map, 
 bit shift map, 2x mod 1 map, Bernoulli map, doubling map or sawtooth map "
(let* ((n (numerator ratio-angle))
       (d (denominator ratio-angle)))
  (setq n  (mod (* n 2) d)) ; (2 x n) modulo d = doubling
  (/ n d)))

(defun give-period (ratio-angle)
	"gives period of angle in turns (ratio) under doubling map"
	(let* ((n (numerator ratio-angle))
	       (d (denominator ratio-angle))
	       (temp n)) ; temporary numerator
	  
	  (loop for p from 1 to 100 do 
		(setq temp  (mod (* temp 2) d)) ; (2 x n) modulo d = doubling)
		when ( or (= temp n) (= temp 0)) return p )))

(defun give-list (period)
  "Returns a list of all  angles of  given period 
   without angles with lower periods, which divide period.
   period is integer > 0 "
  (let* ((angles-list '())
	 (den (- (expt 2 period) 1 )) ; denominator of angle = (2^period)-1
	  a ) ; angle
    (when (> period 0) 
      (dotimes  (num (+ den 1)) ; from 0 to den
	(setq a (/ num den ))
	(when (= period (give-period a)) ; 
         (setq angles-list (append angles-list (list a)))))) ; 
      angles-list))

(defun not-crosses (new-pair old-pair)
"checks if new arc ( = between new-angle-1 and new-angle-2)
crosses old arc ( = between (first old-pair) and (second old-pair)).
It is a part of Lavaurs algorithm.
Angles are external angles mesured in turns 
angle-1 < angle-2 
test : 
(setq old-p '(1/3 2/3))
(setq new-p '(4/15 6/15))
(not-crosses new-p old-p) 
(not-crosses (list 1/7 2/7) old-p) ; t
(not-crosses (list 1/7 3/7) old-p) ; nil
"
(let ((old-angle-1 (first old-pair))
	(old-angle-2 (second old-pair))
	(new-angle-1 (first new-pair))
	(new-angle-2 (second new-pair)))

; check the order of input angles
(when (< new-angle-2 new-angle-1) (rotatef new-angle-1 new-angle-2))
(when (< old-angle-2 old-angle-1) (rotatef old-angle-1 old-angle-2))
(cond
	((< new-angle-1 new-angle-2 old-angle-1 old-angle-2) t) 
	((< old-angle-1 old-angle-2 new-angle-1 new-angle-2) t) 
	((< new-angle-1 old-angle-1 old-angle-2 new-angle-2) t) 
	((< old-angle-1 new-angle-1 new-angle-2 old-angle-2) t)
	(t nil))))

(defun not-crosses-all (new-pair old-list)
"checks if new pair of rays do not crosses any of pairs from old list
test : 
(setq old-pairs '((1/3 2/3) (1/7 2/7)))
(not-crosses-all (list 4/15 6/15) old-pairs) ; nil
"
(let ((b-result T)) 
(loop for old-pair in old-list do (setq b-result (and b-result (not-crosses new-pair old-pair))))
b-result ))

(defun give-pairs-up-to (period)
"gives list of external angles  (pairs landing on the same root point)
for periods up to input period
period >= 3 
examples of use : 
(give-pairs-old-old 3)
(give-pairs-old-old 7)
"
 (let* ((pairs (list (list 1/3 2/3))) ; list for period 2 
	angles-list ; temporary list
	i	
	new-first-angle
	new-second-angle) 

	( loop for p from 3 to period do 
  		(setq angles-list (give-list p))
		(loop for j from  0 to (/ (- (length angles-list) 1) 2)  do 

  			(setq new-first-angle (pop angles-list)) ; pop removes angle from list
			; find second ( conjugate) angle
			(setq i 0)
			(loop  do ;for i from 0 to (- (length angles-list) 1) do  ; first = nth 0 
 		
				(setq new-second-angle (nth i angles-list)) ; nth do not removes angle from list
				(setq i (+ i 1)) 
  			until (not-crosses-all (list new-first-angle new-second-angle) pairs))
		(setq angles-list (remove new-second-angle angles-list))
	(push (list new-first-angle new-second-angle)  pairs))) ; save new pair to pairs list

(reverse pairs)))

; it should be the same as number of components
;(loop for p from 3 to 10 do (print (length (give-pairs p))))
;3 
;6 
;15 
;27 
;63 
;120 
;252 
;495 

; 
(defun give-pairs (period-max)
"gives list of external angles  (pairs landing on the same root point)
for period = pariod-max
period >= 2 
examples of use : 
(give-pairs 3)
(give-pairs 7) ; 
(time (give-pairs 16))
(compile 'give-pairs)

"
 (let* ((pairs (list (list 1/3 2/3))) ; list for period 2 
	angles-list ; temporary list
	(previous-list pairs) ; list for "previous period" = (period -1)
	i	
	new-first-angle
	new-second-angle) 

	( loop for period from 3 to period-max do 
  		(setq angles-list (give-list period)) ; find all angles for given period
		(setq previous-list pairs) ; update previous list
		; match pairs of angles 
		(loop for j from  0 to (/ (- (length angles-list) 1) 2)  do 

  			(setq new-first-angle (pop angles-list)) ; pop removes angle from list
			; find second ( conjugate) angle
			(setq i 0)
			(loop  do ;for i from 0 to (- (length angles-list) 1) do  ; first = nth 0 
 		
				(setq new-second-angle (nth i angles-list)) ; nth do not removes angle from list
				(setq i (+ i 1)) 
  			until (not-crosses-all (list new-first-angle new-second-angle) pairs))
			(setq angles-list (remove new-second-angle angles-list))
			(push (list new-first-angle new-second-angle)  pairs))); save new pair to pairs list
		
	 
	(setq pairs (set-difference pairs previous-list  :test 'equal))	; remove previous angles
	(reverse pairs)))

; --------------------------------------  drawing code ------------------------------------------------

(defun ttr (turn)           
" Turns to Radians"
(* turn  (* 2 pi) ))

(defun give-arc-list (circle-list angle-list)
  "
  Copyright 2009 Rubén Berenguel
  ruben /at/ maia /dot/ ub /dot/ es

  Find the ortogonal circle to the main circle, given the angles in
  it. 
  Input : 
  R: radius of the main circle 
  angle1, angle2 :  angles of main circles (in turns)
  (a, b) , (ba, bb) : points of main circle and new ortogonal circle
  Output is a list for svg path procedure
  thru draw-arc procedure

  http://classes.yale.edu/fractals/Labs/NonLinTessLab/BasicConstr3.html 

  With minor changes by Adam Majewski   " 
  (let* ((x0 (first circle-list))
	 (y0 (second circle-list))
	 (r0 (third circle-list))
	 (alpha (ttr ( first angle-list))) ; convert units from turns to radians
	 (balpha (ttr (second angle-list)))
	 (gamma (+ alpha (/ (- balpha alpha) 2))) ; angle between alpha and balpha
         (ca (cos alpha))
	 (cg (cos gamma))
	 (sa (sin alpha))
	 (sg (sin gamma))
	 (temp (/ r0 (+ (* ca cg) (* sa sg))))
         ; first common point 
	 (a (+ x0 (* r0 ca))) ; a = x0 + r0 * cos(alpha)
	 (b (+ y0 (* r0 sa))) ; b = y0 + r0 * sin(alpha)
	 ; second common point 
	 (ba (+ x0 (* r0 (cos balpha)))) ; ba = x0 + r0 * cos(balpha)
	 (bb (+ y0 (* r0 (sin balpha)))) ; bb = y0 + r0 * sin(balpha)	
	 ; center of ortogonal circle
	 (x (+ x0 (* temp cg)))
	 (y (+ y0 (* temp sg)))
	 ; center of middle circle 
	 (xma (- x a))
	 (ymb (- y b))
	 ; radius of ortogonal circle
	 (r (sqrt (+ (* xma xma) (* ymb ymb))))
	 ; where write labals of arcs
	 (rt (+ r0 35))
	 ; first common point fot label
	 (at (+ x0 (* rt ca))) 
	 (bt (+ y0 (* rt sa))) 
	 ; second common point for label
	 (bat(+ x0 (* rt (cos balpha)))) 
	 (bbt (+ y0 (* rt (sin balpha)))))
	 ; result
	 (list 	a b ; first point of arc = current (a,b)
                r ; radius
	  	ba bb ; last point of arc
		at bt
		bat bbt)))

(defun draw-arc (stream-name circle-list angle-list)
" computes otogonal circle
  using give-arc-list
  and draws arc using svg path command, from the current point to (x, y)
 M = Move to current point ( here (a,b))
 A = elliptical Arc : rx ry   x-axis-rotation large-arc-flag sweep-flag x y"

(let* ((arc-list (give-arc-list circle-list angle-list)))
 (format stream-name "<path d=\"M~,0f ~,0f A~,0f ~,0f 0 0 0 ~,0f ~,0f\"  />~%" 
	(first arc-list)
	(second arc-list)
	(third arc-list)
	(third arc-list)
	(fourth arc-list)  ; 
 	(fifth arc-list))
  ; write labels ( angles) of arcs to image
  (format stream-name "<text x=\"~,0f\"  y=\"~,0f\">~a</text>" (sixth arc-list) (seventh arc-list) (write-to-string (first angle-list)))
  (format stream-name "<text x=\"~,0f\"  y=\"~,0f\">~a</text>" (eighth arc-list) (ninth arc-list) (write-to-string (second angle-list)))))

(defun draw-arcs (stream-name circle-list angles-list)
"draws arc from angles-list
using draw-arc procedure"
(loop for angles in angles-list do (draw-arc stream-name circle-list angles)))

; example of use : (draw-lavaurs "a.svg" 800 2)

(defun draw-lavaurs (file-name side period)
"computes 
 and draws Lavaurs model of Mandelbrot set.  "

  (let* ((x0 (/ side 2))
         (y0 x0)  ; 
 	 (r0 (- x0 50)) ; leave place for titles
	 (main-circle-list (list x0 y0 r0)))
            

            (with-open-file 
			(st file-name 
			:direction :output
			:if-exists :SUPERSEDE
			:if-does-not-exist :create )
       		; write  file header to the file
		(format st "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>~%")
		(format st "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"~%")
                (format st "\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">~%")
		(format st "<svg width=\"~dcm\" height=\"~dcm\" viewBox=\"0 0 ~d ~d\" ~%" 20 20 side side)
                (format st "xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">~%" )

		; draw main circle 	
		(format st "<circle cx=\"~f\" cy=\"~f\" r=\"~f\" fill=\"none\" stroke=\"black\" />~%" x0 y0 r0)		

		; compute and draw arcs ( chords)		
		(format st "<g  fill=\"none\" stroke=\"black\" stroke-width=\"1\">~%") ; open group
		(draw-arcs st main-circle-list (give-pairs-up-to (- period 1))) ; draw	arcs
		(format st "</g>~%") ; close group
		(format st "<g  fill=\"none\" stroke=\"red\" stroke-width=\"1\">~%") ; open group
		(draw-arcs st main-circle-list (give-pairs period)) ; draw	arcs
		(format st "</g>~%") ; close group

	(format st "</svg>~%") ; close svg
(format t "file ~S is saved ~%" file-name) ; info
        

)))

;----------global var ----------------------
 
(defparameter *period* 5 " maximal period. It is an integer >= 2 ")

(defparameter *size* 1000 " size of image in pixels. It is an integer >= 0 ") 

(defparameter *file-name*
  (make-pathname 
   :name (concatenate 'string "lavaurs-t-" (write-to-string *period*))
   :type "svg")
  "name (or pathname) of png file ")
 

;======================= run =====================================================================

(draw-lavaurs *file-name* *size* *period*)