'============================================================================
' The Mandelbrot Dazibao - October 2003 - http://mandelbrot.dazibao.free.fr
'       Von Koch, Sierpinsky and Apollonia fractals in Screen 12 mode
'============================================================================
' Copyright 2003 - The Mandelbrot Dazibao - http://mandelbrot.dazibao.free.fr
'
'  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 2 of the License, or
'  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, write to the Free Software
'  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
'============================================================================
DECLARE SUB Init.Circle (x0!, y0!, Radius!)
DECLARE SUB Init.Line (x0!, y0!, xf!, yf!, Rank%)
DECLARE SUB Init.Triangle (x0!, y0!, Span!, Rank%)

DECLARE SUB Plot.Circle (Id%)
DECLARE SUB Plot.Line (Id%)
DECLARE SUB Plot.Triangle (Id%)

DECLARE SUB Apollo.Create ()
DECLARE SUB Gasket.Create ()
DECLARE SUB Koch.Create ()
DECLARE SUB Sierpinsky.Create ()

DECLARE SUB Apollo.Plot ()
DECLARE SUB Gasket.Plot ()
DECLARE SUB Koch.Plot ()
DECLARE SUB Sierpinsky.Plot ()

DECLARE SUB Soddy.Circle (Id1%, Id2%, Id3%)
DECLARE SUB Split.Line (Id%)
DECLARE SUB Split.Triangle (Id%)

' Object for Apollonia fractals
TYPE Disk
x0 AS SINGLE
y0 AS SINGLE
Radius AS SINGLE
Id1 AS INTEGER
Id2 AS INTEGER
Id3 AS INTEGER
Rank AS INTEGER
END TYPE

' Object for Von Koch fractal
TYPE KochLine
x0 AS SINGLE
y0 AS SINGLE
xf AS SINGLE
yf AS SINGLE
Rank AS INTEGER
END TYPE

' Object for Sierpinsky's triangle
TYPE Triangle
x0 AS SINGLE
y0 AS SINGLE
Span AS SINGLE
Rank AS INTEGER
END TYPE

DIM SHARED NbObjects%, GraphType%, Disk(3275) AS Disk
DIM SHARED NbLines%, KochLine(3639) AS KochLine
DIM SHARED NbTriangles%, Triangle(4680) AS Triangle

' Program constants
CONST ScrWidth = 640, ScrHeight = 480, Pi = 3.1415927#

' Activate graphic mode
CLS
SCREEN 12

' Plot the fractals...
Sierpinsky.Plot
Koch.Plot
Apollo.Plot
Gasket.Plot

SYSTEM

SUB Apollo.Create

GraphType% = 1: NbObjects% = 0

lRef = 1.5: RefRad = lRef * SQR(3) / 2
Init.Circle lRef / 2, RefRad, RefRad
Init.Circle lRef / 2, -RefRad, RefRad
Init.Circle -lRef, 0, RefRad

FOR Id% = 1 TO 3
Disk(Id%).Id1 = 0: Disk(Id%).Id2 = 0: Disk(Id%).Id3 = 0: Disk(Id%).Rank = 0
NEXT Id%

Soddy.Circle 1, 2, 3

FOR Rank% = 1 TO 6
FOR Id% = 4 TO NbObjects%
IF Disk(Id%).Rank = Rank% THEN
Soddy.Circle Id%, Disk(Id%).Id1%, Disk(Id%).Id2%
Soddy.Circle Id%, Disk(Id%).Id1%, Disk(Id%).Id3%
Soddy.Circle Id%, Disk(Id%).Id2%, Disk(Id%).Id3%
END IF
NEXT Id%
NEXT Rank%

END SUB

SUB Apollo.Plot
Apollo.Create

FOR Rank% = 0 TO 7
CLS : Count% = 0
    FOR Id% = 1 TO NbObjects%
    IF Disk(Id%).Rank <= Rank% THEN Plot.Circle Id%: Count% = Count% + 1
    NEXT Id%
PRINT Count%; "circles"
DO WHILE INKEY$ = "": LOOP
NEXT Rank%

END SUB

SUB Gasket.Create
GraphType% = 2: NbObjects% = 0

lRef = .3
RefRad = lRef * SQR(3) / 2
Init.Circle 0, 0, lRef + RefRad
Init.Circle lRef / 2, RefRad, RefRad
Init.Circle lRef / 2, -RefRad, RefRad
Init.Circle -lRef, 0, RefRad

FOR Id% = 1 TO 4
Disk(Id%).Id1 = 0: Disk(Id%).Id2 = 0: Disk(Id%).Id3 = 0: Disk(Id%).Rank = 0
NEXT Id%

Soddy.Circle 2, 3, 4
Soddy.Circle 1, 2, 3
Soddy.Circle 1, 2, 4
Soddy.Circle 1, 3, 4

FOR Rank% = 1 TO 5
FOR Id% = 5 TO NbObjects%
IF Disk(Id%).Rank = Rank% THEN
Soddy.Circle Id%, Disk(Id%).Id1%, Disk(Id%).Id2%
Soddy.Circle Id%, Disk(Id%).Id1%, Disk(Id%).Id3%
Soddy.Circle Id%, Disk(Id%).Id2%, Disk(Id%).Id3%
END IF
NEXT Id%
NEXT Rank%

END SUB

SUB Gasket.Plot
Gasket.Create
FOR Rank% = 0 TO 6
CLS : Count% = 0
    FOR Id% = 1 TO NbObjects%
    IF Disk(Id%).Rank <= Rank% THEN Plot.Circle Id%: Count% = Count% + 1
    NEXT Id%
PRINT Count%; "circles"
DO WHILE INKEY$ = "": LOOP
NEXT Rank%
END SUB

SUB Init.Circle (x0!, y0!, Radius)
NbObjects% = NbObjects% + 1
Disk(NbObjects%).Radius = Radius
Disk(NbObjects%).x0 = x0!
Disk(NbObjects%).y0 = y0!
END SUB

SUB Init.Line (x0!, y0!, xf!, yf!, Rank%)
NbLines% = NbLines% + 1
KochLine(NbLines%).x0 = x0!
KochLine(NbLines%).y0 = y0!
KochLine(NbLines%).xf = xf!
KochLine(NbLines%).yf = yf!
KochLine(NbLines%).Rank = Rank%
END SUB

SUB Init.Triangle (x0!, y0!, Span!, Rank%)
NbTriangles% = NbTriangles% + 1
Id% = NbTriangles%
Triangle(Id%).x0 = x0!
Triangle(Id%).y0 = y0!
Triangle(Id%).Span = Span!
Triangle(Id%).Rank = Rank%
END SUB

SUB Koch.Create

NbLines% = 0
lRef = .4: yRef = lRef * TAN(Pi / 6)
Init.Line -lRef, yRef, lRef, yRef, 0

FOR Rank% = 1 TO 5
FOR Id% = 1 TO NbLines%
Split.Line Id%
NEXT Id%
NEXT Rank%

END SUB

SUB Koch.Plot
Koch.Create

FOR Rank% = 0 TO 5
CLS : Count% = 0
lRef = .4: yRef = lRef * TAN(Pi / 6)
Radius% = INT(SQR(yRef ^ 2 + lRef ^ 2) * 480 + .5)
CIRCLE (320, 240), Radius%, 14
    FOR Id% = 1 TO NbLines%
    IF KochLine(Id%).Rank = Rank% THEN Count% = Count% + 1: Plot.Line Id%
    NEXT Id%
PRINT Count% * 3; "lines"
DO WHILE INKEY$ = "": LOOP
NEXT Rank%

END SUB

SUB Plot.Circle (Id%)
xp0% = INT(Disk(Id%).x0 * ScrHeight / 1.2 + .5) + ScrWidth / 2
yp0% = INT(Disk(Id%).y0 * ScrHeight / 1.2 + .5) + ScrHeight / 2
Rad0% = INT(Disk(Id%).Radius * ScrHeight / 1.2 + .5)
CIRCLE (xp0%, yp0%), Rad0%, Disk(Id%).Rank + 1
END SUB

SUB Plot.Line (Id%)
x0! = KochLine(Id%).x0
y0! = KochLine(Id%).y0
xf! = KochLine(Id%).xf
yf! = KochLine(Id%).yf

x10! = -.5 * x0! + SQR(3) / 2 * y0!
y10! = -SQR(3) / 2 * x0! - .5 * y0!
x1f! = -.5 * xf! + SQR(3) / 2 * yf!
y1f! = -SQR(3) / 2 * xf! - .5 * yf!

x20! = -.5 * x10! + SQR(3) / 2 * y10!
y20! = -SQR(3) / 2 * x10! - .5 * y10!
x2f! = -.5 * x1f! + SQR(3) / 2 * y1f!
y2f! = -SQR(3) / 2 * x1f! - .5 * y1f!

xp0% = INT(x0! * ScrHeight + .5) + ScrWidth / 2
yp0% = INT(y0! * ScrHeight + .5) + ScrHeight / 2
xpf% = INT(xf! * ScrHeight + .5) + ScrWidth / 2
ypf% = INT(yf! * ScrHeight + .5) + ScrHeight / 2

xp10% = INT(x10! * ScrHeight + .5) + ScrWidth / 2
yp10% = INT(y10! * ScrHeight + .5) + ScrHeight / 2
xp1f% = INT(x1f! * ScrHeight + .5) + ScrWidth / 2
yp1f% = INT(y1f! * ScrHeight + .5) + ScrHeight / 2

xp20% = INT(x20! * ScrHeight + .5) + ScrWidth / 2
yp20% = INT(y20! * ScrHeight + .5) + ScrHeight / 2
xp2f% = INT(x2f! * ScrHeight + .5) + ScrWidth / 2
yp2f% = INT(y2f! * ScrHeight + .5) + ScrHeight / 2

LINE (xp0%, yp0%)-(xpf%, ypf%), 2
LINE (xp10%, yp10%)-(xp1f%, yp1f%), 2
LINE (xp20%, yp20%)-(xp2f%, yp2f%), 2

END SUB

SUB Plot.Triangle (Id%)
' Based on lines. Could be reprogged with DRAW

x0! = Triangle(Id%).x0
y0! = Triangle(Id%).y0
Span! = Triangle(Id%).Span
Rank% = Triangle(Id%).Rank

x1! = x0! + Span! / 2
y1! = y0! + Span! / 2 * TAN(Pi / 6)
x2! = x0! - Span! / 2
y2! = y1!
x3! = x0!
y3! = y0! - Span! * SQR(.5 ^ 2 + (.5 * TAN(Pi / 6)) ^ 2)

yp1% = INT(y1! * ScrHeight + .5) + ScrHeight / 2
yp3% = INT(y3! * ScrHeight + .5) + ScrHeight / 2

Slope1! = (x3! - x1!) / (y3! - y1!)
Slope2! = (x3! - x2!) / (y3! - y2!)
FOR y% = yp1% TO yp3% STEP SGN(yp3% - yp1%)
y! = (y% - ScrHeight / 2) / ScrHeight
xl1! = x1! + (y! - y1!) * Slope1!
xl2! = x2! + (y! - y2!) * Slope2!
xp1% = INT(xl1! * ScrHeight + .5) + ScrWidth / 2
yp2% = INT(yl1! * ScrHeight + .5) + ScrHeight / 2
xp2% = INT(xl2! * ScrHeight + .5) + ScrWidth / 2
yp2% = INT(yl2! * ScrHeight + .5) + ScrHeight / 2
LINE (xp1%, y%)-(xp2%, y%), Rank% + 1
NEXT y%

END SUB

SUB Sierpinsky.Create
Init.Triangle 0, 0, .7, 0

FOR Rank% = 0 TO 6
    FOR Id% = 1 TO NbTriangles%
    IF Triangle(Id%).Rank = Rank% THEN Split.Triangle Id%
    NEXT Id%
NEXT Rank%

END SUB

SUB Sierpinsky.Plot
Sierpinsky.Create
FOR Rank% = 0 TO 6
CLS
Count% = 0
    FOR Id% = 1 TO NbTriangles%
    IF Triangle(Id%).Rank = Rank% THEN Plot.Triangle Id%: Count% = Count% + 1
    NEXT Id%
PRINT Count%; "Triangles"
DO WHILE INKEY$ = "": LOOP
NEXT Rank%
END SUB

SUB Soddy.Circle (Id1%, Id2%, Id3%)

' Get the three ranks
Rank1% = Disk(Id1%).Rank
Rank2% = Disk(Id2%).Rank
Rank3% = Disk(Id3%).Rank
IF Rank1% < Rank2% THEN SWAP Rank1%, Rank2%
IF Rank1% < Rank3% THEN SWAP Rank1%, Rank3%
Rank% = Rank1%

' Get the three radius
R1! = Disk(Id1%).Radius
R2! = Disk(Id2%).Radius
R3! = Disk(Id3%).Radius

' Outer Soddy Circle
IF GraphType% = 2 THEN
IF Id1% = 1 THEN R1! = -R1!
IF Id2% = 1 THEN R2! = -R2!
IF Id3% = 1 THEN R3! = -R3!
END IF


' Get the three centers coordinates
x1! = Disk(Id1%).x0
x2! = Disk(Id2%).x0
x3! = Disk(Id3%).x0
y1! = Disk(Id1%).y0
y2! = Disk(Id2%).y0
y3! = Disk(Id3%).y0

' The three equations to be solved are:
' (x-x1)^2 + (y-y1)^2 = (R + R1)^2      (1)
' (x-x2)^2 + (y-y2)^2 = (R + R2)^2      (2)
' (x-x3)^2 + (y-y3)^2 = (R + R3)^2      (3)

' Substract the equations:
' (1) - (2):
' (x-x1)^2 + (y-y1)^2 - (x-x2)^2 - (y-y2)^2 = (R + R1)^2 - (R + R2)^2
' -2*x*x1 + x1^2 - 2*y*y1 + y1^2 + 2*x*x2 - x2^2 + 2*y*y2 - y2^2 = 2*(R1-R2)*R +(R1^2-R2^2)
' x*2*(x2-x1) + y*2*(y2-y1) + (x1^2-x2^2) + (y1^2-y2^2) =  2*(R1-R2)*R +(R1^2-R2^2)

' (A) x*2*(x2-x1) + y*2*(y2-y1) + R*2*(R2-R1) = (x2^2-x1^2) + (y2^2-y1^2) + (R1^2-R2^2)
' (B) x*2*(x3-x2) + y*2*(y3-y2) + R*2*(R3-R2) = (x3^2-x2^2) + (y3^2-y2^2) + (R2^2-R3^2)
' (C) x*2*(x1-x3) + y*2*(y1-y3) + R*2*(R1-R3) = (x1^2-x3^2) + (y1^2-y3^2) + (R3^2-R1^2)

' x * Ax + y * Ay + r * Ar = Aa
' x * Bx + y * By + r * Br = Bb
' x * Cx + y * Cy + r * Cr = Cc

ax! = 2 * (x2! - x1!): Ay! = 2 * (y2! - y1!): Ar! = 2 * (R2! - R1!)
bx! = 2 * (x3! - x2!): By! = 2 * (y3! - y2!): Br! = 2 * (R3! - R2!)
cx! = 2 * (x1! - x3!): Cy! = 2 * (y1! - y3!): Cr! = 2 * (R1! - R3!)
Aa! = (x2! ^ 2 - x1! ^ 2) + (y2! ^ 2 - y1! ^ 2) + (R1! ^ 2 - R2! ^ 2)
Bb! = (x3! ^ 2 - x2! ^ 2) + (y3! ^ 2 - y2! ^ 2) + (R2! ^ 2 - R3! ^ 2)
Cc! = (x1! ^ 2 - x3! ^ 2) + (y1! ^ 2 - y3! ^ 2) + (R3! ^ 2 - R1! ^ 2)

DAB! = ax! * By! - Ay! * bx!
DAC! = ax! * Cy! - Ay! * cx!
DBC! = bx! * Cy! - By! * cx!

' 1st case : R1 = R2 = R3
IF Ar! = 0 AND Br! = 0 THEN
Det! = ax! * By! - bx! * Ay!
x! = (Aa! * By! - Bb! * Ay!) / Det!
y! = (ax! * Bb! - Ay! * Aa!) / Det!
p! = 1
q! = 2 * R1!
r! = R1 ^ 2 - (x! - x1!) ^ 2 - (y! - y1!) ^ 2
d! = q! ^ 2 - 4 * p! * r!
Radius! = (-q! + SQR(d!)) / 2
GOTO NewDisk:
END IF

' 2nd case : R1, R2 and R3 are <>
IF DAB! <> 0 THEN
' x * Ax + y * Ay = Aa - r * Ar
' x * Bx + y * By = Bb -  r * Br
' x = ((Aa - r * Ar) * By - (Bb - r * Br) * Ay) / DAB!
' y = (Ax * (Bb - r * Br) - Bx * (Aa - r * Ar)) / DAB!
' x = (r*(Br*Ay-Ar*By)+(Aa*By-Bb*Ay))/DAB!
' y = (r*(Ar*Bx-Br*Ax)+(Bb*Ax-Aa*Bx))/DAB!
' x = xr*r+xc
' y = yr*r+yc

xr! = (Br! * Ay! - Ar! * By!) / DAB!
yr! = (Ar! * bx! - Br! * ax!) / DAB!
xc! = (Aa! * By! - Bb! * Ay!) / DAB!
yc! = (Bb! * ax! - Aa! * bx!) / DAB!

' (x-x3)^2 + (y-y3)^2 = (R + R3)^2
' x^2-2*x*x3*x3^2 + y^2-2*y*y3*y3^2 = r^2 + 2*r*r3 + r3^2
' (xr*r+xc)^2 - 2*(xr*r+xc)*x3 + x3^2 + (yr*r+yc)^2 - 2*(yr*r+yc)*y3 + y3^2 = r^2 + 2*r*r3 + r3^2

Zp! = xr! ^ 2 + yr! ^ 2 - 1
zq! = 2 * xc! * xr! - 2 * xr! * x3! + 2 * yc! * yr! - 2 * yr! * y3! - 2 * R3!
zr! = xc! ^ 2 - 2 * xc! * x3! + x3! ^ 2 + yc! ^ 2 - 2 * yc! * y3! + y3! ^ 2 - R3! ^ 2

Det! = zq! ^ 2 - 4 * Zp! * zr!

Radius! = (-zq! - SQR(Det!)) / 2 / Zp!

x! = xr! * Radius! + xc!
y! = yr! * Radius! + yc!


GOTO NewDisk:

END IF

NewDisk:
Init.Circle x!, y!, Radius!
Disk(NbObjects%).Id1 = Id1%
Disk(NbObjects%).Id2 = Id2%
Disk(NbObjects%).Id3 = Id3%
Disk(NbObjects%).Rank = Rank% + 1
EXIT SUB

END SUB

SUB Split.Line (Id%)
x0! = KochLine(Id%).x0
y0! = KochLine(Id%).y0
xf! = KochLine(Id%).xf
yf! = KochLine(Id%).yf
Rank% = KochLine(Id%).Rank

x1! = x0! + (xf! - x0!) / 3
y1! = y0! + (yf! - y0!) / 3
x2! = x0! + (xf! - x0!) * 2 / 3
y2! = y0! + (yf! - y0!) * 2 / 3

x4! = (x0! + xf!) / 2 + (y0! - yf!) / 2 / SQR(3)
y4! = (y0! + yf!) / 2 + (xf! - x0!) / 2 / SQR(3)

Init.Line x0!, y0!, x1!, y1!, Rank% + 1
Init.Line x1!, y1!, x4!, y4!, Rank% + 1
Init.Line x4!, y4!, x2!, y2!, Rank% + 1
Init.Line x2!, y2!, xf!, yf!, Rank% + 1
END SUB

SUB Split.Triangle (Id%)

x0! = Triangle(Id%).x0
y0! = Triangle(Id%).y0
Span! = Triangle(Id%).Span
Rank% = Triangle(Id%).Rank

x1! = x0! + Span! / 2: x2! = x0! - Span! / 2
y1! = y0! + Span! / 2 * TAN(Pi / 6): y2! = y1!
x3! = x0!: y3! = y0! - Span! * SQR(.5 ^ 2 + (.5 * TAN(Pi / 6)) ^ 2)

xa! = (x2! + x3!) / 2: ya! = (y2! + y3!) / 2
xb! = (x3! + x1!) / 2: yb! = (y3! + y1!) / 2
xc! = (x1! + x2!) / 2: yc! = (y1! + y2!) / 2

' 1st Triangle 1, b, c
xn! = (x1! + xc!) / 2: yn! = y1! - Span! / 4 / SQR(3)
Init.Triangle xn!, yn!, Span! / 2, Rank% + 1
' 2nd Triangle 2, a, c
xn! = (x2! + xc!) / 2: yn! = y2! - Span! / 4 / SQR(3)
Init.Triangle xn!, yn!, Span! / 2, Rank% + 1
' 3rd Triangle 3, a, b
xn! = (xa! + xb!) / 2: yn! = ya! - Span! / 4 / SQR(3)
Init.Triangle xn!, yn!, Span! / 2, Rank% + 1

END SUB