; 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 Lavaures : ; 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 "ls.lisp") ; (draw-lavaurs "lavaurs-5.png" 2000 5) ; ;look for lavaurs-5.png file in your home directory ; ; 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/>.  ; first run ;(require :asdf) ;(require :asdf-install) ;(asdf-install:install :zpng) ;(asdf-install:install :Vecto) ; http://www.xach.com/lisp/vecto/ ; you must press 2 and 0 when the program asks    ; next times load packages from disk (asdf:operate 'asdf:load-op 'vecto)     (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) ))  ; circle-list angle-list (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 vecto:arc procedure   thru draw-arc procedure    http://classes.yale.edu/fractals/Labs/NonLinTessLab/BasicConstr3.html   "    (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) 	 ; 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)))) 	 ; angle of point (a,b) measured in new circle units 	 (phi  (atan r0 r))) 	 ; result 	 (list 	x y r  	  	(+ pi gamma phi)  ; new balpha            	(- (+ pi gamma) phi) ; new alpha 	  	a b))) ; point (a,b)  (defun draw-arc (circle-list angle-list) " computes otogonal circle   using give-arc-list   and draws arc using vecto:arcn procedure  vecto:arcn x y radius angle1 angle2 "  (let* ((arc-list (give-arc-list circle-list angle-list))) (vecto:move-to ( sixth arc-list) (seventh arc-list)) ; beginning of arc is point (a,b) (vecto:arcn 	( first arc-list)  ; x 	(second arc-list)  ; y 	(third arc-list)   ; radius 	(fourth arc-list)  ; angle1  	(fifth arc-list))) ; angle2   (vecto:stroke))  (defun draw-arcs (circle-list angles-list) "draws arc from angles-list using draw-arc procedure" (loop for angles in angles-list do (draw-arc circle-list angles)))  ; example of use : (draw-lavaurs "a.png" 800 2)  (defun draw-lavaurs (file side period) "computes   and draws Lavaurs model of Mandelbrot set.  "    (vecto:with-canvas (:width side :height side) ; vecto 	       (vecto:set-rgb-stroke 0 0 0) ; vecto 	       (vecto:set-line-width 1) ; vecto     	      (let* (	(x0 (/ side 2))                         (y0 x0)  ;   			(r0 (- x0 50)) ; leave place for titles 			(main-circle-list (list x0 y0 r0)) 			arc-list)  		; main circle 	 		 (vecto:centered-circle-path x0 y0 r0 ) 		  		; arcs ( chords)		 		(setq arc-list (give-pairs-up-to period)) ; compute 		(draw-arcs main-circle-list arc-list) ; draw	 		 		(vecto:stroke)) ; before save ( vecto procedure) 		 	       (print (vecto:save-png file)))) ; save image  ;----------global var ----------------------   (defparameter *period* 4 " maximal period. It is an integer >= 2 ")  (defparameter *size* 800 " size of image in pixels. It is an integer >= 0 ")   (defparameter *file-name*   (make-pathname     :name (concatenate 'string "lavaurs-" (write-to-string *period*))    :type "png")   "name (or pathname) of png file ")    ;======================= run =====================================================================  (draw-lavaurs *file-name* *size* *period*)