I, the copyright holder of this work, hereby publish it under the following licenses:

This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license.

You are free:

• to share – to copy, distribute and transmit the work
• to remix – to adapt the work

Under the following conditions:

• attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
• share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.

https://creativecommons.org/licenses/by-sa/3.0CC BY-SA 3.0 Creative Commons Attribution-Share Alike 3.0 truetrue

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled GNU Free Documentation License.http://www.gnu.org/copyleft/fdl.htmlGFDLGNU Free Documentation Licensetruetrue

You may select the license of your choice.

/* 
computes and draws to png file apollonian gasket

initial configuration : (k1,k1,k1)
- step -1 = three equal mutually tangent circles  (curvature : k1 > 0 and ri:r0out/a1 ) packed into outer circle
- step  0 = 2 circles : one outer circles ( curvature : k0out <0) and one most inner circle  ( curvature : k0in > 0)

a1 value is from : http://www2.stetson.edu/~efriedma/cirincir/

Uses Soddy algorithm ( finds inner and outer Soddy circle using circle Descartes theorem)
http://langexplr.blogspot.com/2007/10/drawing-apollonian-gasket-with-common.html

for each step  and circle :
- make ck list
- save data to svg file

Note that initial SVG coordinate system has the origin at the top/left 
with the x-axis pointing to the right and the y-axis pointing down.
It is not usual.

When center of outer circle is at origin program is more complicated,
because one has to chose solutions in case of centers.
It is easier when the all gasket is in the first (++) quadrant of Cartesian plane

Adam Majewski  fraktal.republika.pl
2010.07.19
with help of :
 Ed Beroset ,  Jaime Villat, Richard Fateman,  Stavros Macrakis 
based on code by :
 Luis Diego Fallas,  KiHyuck Hong

*/

/* 

Lists :
ck  : [center, kurvature] Basic list 

*/

/* ---------------------------- definitions -----------------------------------------------*/

/* ---------------------------------basic svg functions -----------------------*/

load(stringproc);

BeginSVG(file_name,cm_width,cm_height,i_width,i_height,_strokeWidth):=
block(
	[destination],
	destination:openw(file_name),

	printf(destination,"<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>~%"),

	printf(destination,"<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"~%\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">~%"),

	printf(destination,"<svg width=\"~dcm\" height=\"~dcm\" viewBox=\"0 0 ~d ~d\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">~%", /* no space between number and unit */

	cm_width,cm_height,i_width,i_height),
	printf(destination,"<g  stroke-width=\"~f\">~%",_strokeWidth), /* all shapes in group have the same width and color */

	return(destination)
)$

CircleSVG(dest,_c,_k,_fill,_color):=
printf(dest,"<circle cx=\"~f\" cy=\"~f\" r=\"~f\" fill=\~s stroke=~s />~%", realpart(_c),imagpart(_c),abs(1/_k),_fill,_color)$
	

EndSVG(destination):=
(
printf(destination,"</g>~%"), /* close the group */
printf(destination,"</svg>~%"),
close(destination)
)$

/* this function : 
takes two values c and k , 
makes a list ck form them, 
saves ck to global list ck_sl and 
gives ck as an output */

make_ck(_file,c,k,fill,color):=block
(
[ck],
ck:[rectform(c),k],
CircleSVG(_file,c,k,fill,color), /* save to svg,file */
ck /* gives ck as an output */
)$

/* circle Descates theorem . It  can also be used for complex Descates theorem */
/* rectform is needed for speed */
solve_m_eq(k_1,k_2,k_3):=float(rectform(k_1+k_2+k_3 - 2*sqrt(k_1*k_2 + k_2*k_3 + k_1*k_3)))$/* outer circle */

solve_p_eq(k_1,k_2,k_3):=float(rectform(k_1+k_2+k_3 + 2*sqrt(k_1*k_2 + k_2*k_3 + k_1*k_3)))$/* inner circle */

/* 
gives a ck list for 4-th circle 
using plus and plus solve
*/
give_pp_ck4(ck1,ck2,ck3):=block
(
[c4,k4,ck4],
k4:solve_p_eq(ck1[2],ck2[2],ck3[2]),
c4:solve_p_eq(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
ck4:make_ck(f,c4,k4,none,black)

)$

/* 
gives a ck list for 4-th circle 
using plus and minus solve
*/
/* 
gives a ck list for 4-th circle 
using positive and negative solve
*/
give_pm_ck4(ck1,ck2,ck3):=block
(
[c4,k4,ck4],
k4:solve_p_eq(ck1[2],ck2[2],ck3[2]),
c4:solve_m_eq(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
ck4:make_ck(f,c4,k4,red,black)

);

/* computes  gap and saves it to global var ck_sl
thru give_ ... procedures  */
fill_gap_pp(ckla,cklb,cklc,max_step):=block
(
[ckm],
if max_step>0 then (
		/* circle in the middle */
		ckm:give_pp_ck4(ckla,cklb,cklc),
		/* 3 circles around */
		fill_gap_pp(ckla,cklb,ckm,max_step-1),
		fill_gap_pp(cklb,cklc,ckm,max_step-1),
		fill_gap_pp(ckla,cklc,ckm,max_step-1))
)$

/* it is very specyfic procedure = only for 1c gap */
/* check the order of input values */
fill_hard_gap(ck_ic,ck_ia,ck_0out,max_step):=block
(
[ck_1c,ck_2cc,ck_3ccc,ck_4ccca,ck_5cccac], /* hard circles */
if 	max_step>0 then (
	/* step 1 = circle in the middle */
	ck_1c:give_pm_ck4(ck_ia,ck_ic,ck_0out),      /* 1c */
	/* step 2 = 3 circles around */
	fill_gap_pp(ck_ia,ck_0out,ck_1c,max_step-1), /* 2ca */
	fill_gap_pp(ck_ia,ck_ic,ck_1c,max_step-1), /* 2cb */
	/* hard subgap 2cc */
	if max_step>1 then ck_2cc:give_pm_ck4(ck_ic,ck_1c,ck_0out), /* 2cc */
	/* step 3 for subgap 2cc */
	fill_gap_pp(ck_0out,ck_1c,ck_2cc,max_step-2), /* 3cca */
	fill_gap_pp(ck_ic,ck_1c,ck_2cc,max_step-2),  /* 3ccb */
	if max_step>2 then ck_3ccc:give_pm_ck4(ck_ic,ck_0out,ck_2cc),	/* 3ccc */
	/* step 4 for subgap 3ccc */
	if max_step>3 then ck_4ccca:give_pm_ck4(ck_0out,ck_2cc,ck_3ccc),
	fill_gap_pp(ck_ic,ck_2cc,ck_3ccc,max_step-3),  /* 4cccb */
	fill_gap_pp(ck_ic,ck_0out,ck_3ccc,max_step-3),  /* 4cccc */
	/* step 5 for subgap 4ccca */
	fill_gap_pp(ck_0out,ck_2cc,ck_4ccca,max_step-4),  /*  5cccaa */
	fill_gap_pp(ck_2cc,ck_3ccc,ck_4ccca,max_step-4),  /* 5cccab */
	if max_step>4 then ck_5cccac:give_pm_ck4(ck_0out,ck_3ccc,ck_4ccca)
	)
)$

/*------*/

Apollonian(steps,radius):=
block(

[x0,y0,r0out,a1,ck0out,ck0in, /* circles of stage 0 */
/* 3 equall circles of stage -1 */
ri,ki,
xia,yia,ckia,
xib,yib,ckib,
xic,yic,ckic,
e2,e3,sib,
/* svg */
file_title,
f],

/* ------------- drawing code -----------------------*/
file_title:concat("apollonian_rc_",string(steps),".svg"),
f:BeginSVG(file_title,20,20,2*radius,2*radius,0.8),

/* initial values */
r0out:radius,
x0:r0out,   /* realpart of center */
y0:r0out,   /* imagpart of center */

a1:float(1 + 2 / sqrt(3)), /* http://www2.stetson.edu/~efriedma/cirincir/ */ /* float is needed for solve */
/* step 0 = outer circle  0out with  negative curvature*/
ck0out:make_ck(f,x0+y0*%i,-1/r0out,none,black), /* ck list is a list : center, kurvature */
/* step -1 =  3 equall circles packed in ck0out */
ri:r0out/a1, /* I can compute its common radius */
ki:1/ri,

/* 
I will configure circles in this way :
* one circle in upper row ( realpart of its center will be equal to realpart of outer circle )
* two circles in lower row
I can compute first circle :
*/

xia:x0,
yia:ri,/* one in upper row , remember that in svg y axis is inverted */
ckia:make_ck(f,xia +yia*%i,ki,none,black),
/* now compute second inner circle ib , tangent to circles 0out and ia */

e2:(x0-xib)^2  +(y0-yib)^2  -(r0out-ri)^2=0, /* 0 and ia */
e3:(xia-xib)^2 +(yia-yib)^2 -(ri+ri)^2=0, /* ia and ib */
sib:solve([e2,e3],[xib,yib]),
xib:float(rectform(rhs(sib[2][1]))),
yib:float(rectform(rhs(sib[2][2]))),
ckib:make_ck(f,xib+yib*%i,ki,none,black),
/*  third inner circle ic , tangent to circles 0out , ia  and ib */
xic:float(rectform(rhs(sib[1][1]))),
yic:float(rectform(rhs(sib[1][2]))),
cic:xic+yic*%i,
ckic:make_ck(f,cic,ki,none,black),
/* most inner circle, next from step 0 */
ck0in:give_pp_ck4(ckia,ckib,ckic),
/* fill 6 gaps (curvilinear triangles) */
fill_gap_pp(ckia,ckib,ck0in,steps),
fill_gap_pp(ckib,ckic,ck0in,steps),
fill_gap_pp(ckic,ckia,ck0in,steps),
fill_gap_pp(ckia,ckib,ck0out,steps),
fill_gap_pp(ckib,ckic,ck0out,steps),
fill_hard_gap(ckic,ckia,ck0out,steps),/* !!!!! */

EndSVG(f),
/* ------------------------info --------- */
disp(concat("file  ",file_title, " is made "))
)$

/* --------------------------- compilation ---------------------------------------------------*/

compile(all)$

/* ---------------------------------- computing ----------------------------------------*/

Apollonian(6,500)$

/*
file  apollonian_r_1.svg is made  Evaluation took    	  3.4700 seconds 
file  apollonian_r_2.svg is made  Evaluation took 	  3.3600 seconds 
file  apollonian_r_3.svg is made  Evaluation took 	  3.4400 seconds   
file  apollonian_r_4.svg is made  Evaluation took 	  3.6200 seconds 
file  apollonian_r_5.svg is made  Evaluation took 	  4.1100 seconds 
file  apollonian_r_6.svg is made  Evaluation took 	  5.4200 seconds 
file  apollonian_r_7.svg is made  Evaluation took 	  9.2400 seconds 
file  apollonian_r_8.svg is made  Evaluation took 	 21.5800 seconds 
file  apollonian_r_9.svg is made  Evaluation took 	 55.9000 seconds 
file  apollonian_r_10.svg is made Evaluation took       161.5100 seconds  
file  apollonian_r_11.svg is made Evaluation took  
file  apollonian_r_12.svg is made Evaluation took 

*/