function [] = ap_gasket()
% apollonian_gasket.m:
% draw a Apollonian gasket.
% by David Cary <d.cary@ieee.org>
% I first saw this in
% _The Fractal Geometry of Nature_
% by Benoit Mandelbrot (1983) (pp. 169-172 ?)
% and a quick summary is at
% http://www.astro.virginia.edu/~eww6n/math/ApollonianGasket.html
% which reminds us that
% "The points which are never inside a circle form a set of measure 0
% having fractal dimension approximately 1.3058 (Mandelbrot 1983, p. 172)."
% There is a closed-form solution to the sub-problem
% of finding the circle tangent to 3 circles at
% http://www.astro.virginia.edu/~eww6n/math/ApolloniusProblem.html
% .
% which has a link to
% http://www.astro.virginia.edu/~eww6n/math/notebooks/PlaneGeometry.m
% .

% ap_gasket: draw a
% appolonian gasket

% Default starting condition: equilateral triangle.
% columns x,y. 3 rows for the 3 points.
corners = [...
  0 0;...
  sqrt(3) 1;...
  sqrt(3) -1;...
  ];


% rb = bx^2 + by^2 - axbx - ayby = d2(b) - a dot b
% ra = ax^2 + ay^2 - axbx - ayby = d2(a) - a dot b
% rc = ( d(b-a) - d(a) - d(a) )/2

% For any 2 points a and b,
% |a-b|^2 =
% = (ax-bx)^2 + (ay-by)^2 =
% = ax^2 + bx^2 - 2*ax*bx + ay^2 + by^2 - 2*ay*by =
% = |a|^2 + |b|^2 - 2(a dot b).

% given a triangle with points a, b, c,
% and with line segments opposite them A, B, C,
% we can find the radius of the 3 circles
% centered at a,b,c that are just big enough to touch each other.
% (If you think about making one of them really big,
% and slowly deflating it while inflating the others
% just enough to keep them touching the big one,
% you can see that there is only one place
% where all 3 circles just touch each other).
% Since we know
% A = rb + rc = |b-c|
% B = ra + rc = |a-c|
% C = ra + rb = |a-b|
% therefore
% ra = (B + C - A)/2
% rb = (A + C - B)/2
% rc = (A + B - C)/2.
%
%
% With a new point N, the center of the new circle,
% since we inflate rn just big enough to touch the other circle,
% Na = rn + ra = |a-n|
% Nb = rn + rb = |b-n|
% Nc = rn + rc = |c-n|
% therefore
% rn = (Na + Nb - C)/2
% rn = (Na + Nc - B)/2
% rn = (Nb + Nc - A)/2.
%
% If you consider all possible lines
% between each pair of the 4 points a,b,c,n,
% the sum of the lengths of any 2 line segments that
% do not share a endpoint is special:
% Na + A = Nb + B = Nc + C = rn + ra + rb + rc =
% (A + B + C)/2 + rn = A + B + rn - rc

% if cy,cx = 0,0,
% since we know
% |n| = rn + rc
% |a-n| = rn + ra
% |b-n| = rn + rb
% then we can rearrange to get
% rn = |n| - rc
% and once we eliminate rn, that leaves
% |a-n| = |n| - rc + ra
% |b-n| = |n| - rc + rb.
% which is equivalent to
% |a-n| - |n| = ra - rc
% |b-n| - |n| = rb - rc.
% Really the only 2 unknowns in these 2 equations
% are nx and ny.
% sqrt( (ax-nx)^2 + (ay-ny)^2 ) - sqrt( nx^2 + ny^2) = ra - rc.
% ... umm, that didn't help ...
% square both sides to get
% |a-n|^2 = |n|^2 + (rc + ra)^2 + 2*|n|*(ra-rc)
% |b-n|^2 = |n|^2 + (rc + rb)^2 + 2*|n|*(rb-rc).
% Then we expand to get
% (ax-nx)^2 + (ay-ny)^2 = |n|^2 + (rc+ra)^2 + 2*|n|*(ra-rc) =
% = ax^2 + nx^2 - 2*ax*nx + ay^2 + ny^2 - 2*ay*ny = |n|^2 + 2*|n|*(ra-rc)
% = |a|^2 + |n|^2 - 2*(a dot n) = |n|^2 + 2*|n|*(ra-rc) = 
% = |a|^2 - 2*(a dot n) = 2*|n|*(ra-rc).
% and then square both sides to get
% (|a|^2 - 2*(a dot n))^2 = 4*|n|^2*(ra-rc)^2 = 
% = |a|^4 + 4*(a dot n)^2 - 4*|a|^2*(a dot n) = 4*|n|^2*(ra-rc)^2.
% which we can expand into
% |a|^4 + 4*(ax*nx + ay*ny)^2 - 4*|a|^2*(ax*nx + ay*ny) = 4*(nx^2 + ny^2)*(ra-rc)^2.
% |b|^4 + 4*(bx*nx + by*ny)^2 - 4*|b|^2*(bx*nx + by*ny) = 4*(nx^2 + ny^2)*(rb-rc)^2.
% ... umm, that didn't seem to help either.




% end ap_gasket.m