; console program
; Lisp / Common Lisp / SBCL 
; based on : 
; http://www.mostlymaths.net/2009/08/lavaurs-algorithm.html
; lisp code by R Berenguel

; draws orbit portrait of orbit 
; orbit is a list of pairs of ratios
; for example (list (list 1/7 2/7 ))
;
; arcs are a part of 
; 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 "o.lisp")
; ; look for svg file in your home directory
; ; after loading file you can :
; (draw-orbit "a.svg" 800 (list (list 1/7 2/7 )))
;
; Adam Majewski
; fraktal.republika.pl
; 2010.11.22
;
;
;;  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-orbit (ratio-angle length)
" gives orbit of angle under doubling map.
result is a list of pairs of angles.
but pairs are in ascending order !!!!!
(give-orbit 3/7 3)"
(let* ((old-angle ratio-angle)
       (new-angle (doubling-map old-angle))
       (orbit (list (list old-angle new-angle)))
	(i 1))
(setq old-angle new-angle)
(loop while (< i length) do

	(setq new-angle (doubling-map old-angle))
	(if (< old-angle new-angle) ; for drawing it is better
		(setq orbit (append orbit (list (list old-angle new-angle))))
		(setq orbit (append orbit (list (list new-angle old-angle)))))
	(setq old-angle new-angle)
	(setq i (+ i 1)))
orbit))

; --------------------------------------  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))
	 (yMax (fourth 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 with reversed y axis 
	 (list 	a (- yMax b) ; first point of arc = current (a,b)
                r ; radius
	  	ba (- yMax bb) ; last point of arc
		at (- yMax bt)
		bat (- yMax 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 1 ~,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-orbit "a.svg" 800 (list (list 1/7 2/7 )))

(defun draw-orbit (file-name side orbit)
"draws orbit portrait  "

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

            (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 orbit) ; draw	arcs
		(format st "</g>~%") ; close group

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

)))

;----------global var ----------------------
(defvar *period*  " period of angle ( orbit) under doubling map ")

(defparameter *angle* 1/31 
" external angle in turns. 
 It is a ratio.
  = proper rational rational fraction with odd denominator ") 

(defvar *orbit*  
" Orbit of angle  under doubling map. 
  Orbit is a list of pairs of ratios
  For example (list (list 1/7 2/7 ))")

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

(defparameter *file-name* 
  (make-pathname 
   :name (concatenate 'string "orbit-" (write-to-string (numerator *angle*))"-"(write-to-string (denominator *angle*)))
   :type "svg")
  "name (or pathname) of svg file ")
 

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

(setq *period* (give-period *angle*))
(setq *orbit* (give-orbit *angle* *period*))

(draw-orbit *file-name* *size* *orbit*)