DECLARE SUB master.circles ()
DECLARE SUB invert ()
DECLARE SUB HOW.IT.WORKS ()
DECLARE SUB HELP.SCREEN ()
DECLARE SUB DRAW.CIRCLE ()
COMMON SHARED U0, V0, RO, P, Q, DX, DY, SF, XC, YC, RC
COMMON SHARED N, RF, AZ, FC$, CX, CY, RAD, CC
XC = 170: YC = 100: RC = 95
CONST PI = 3.141593
BEGIN:
CALL HELP.SCREEN
U0 = 0: V0 = 0: r0 = (1 - SIN(PI / N)) / (1 + SIN(PI / N))
P = RF * r0: Q = 0
START.DRAW:
SCREEN 1: CLS
CALL master.circles
RO = SIN(PI / N) / (1 + SIN(PI / N)): rR = 1 / (1 + SIN(PI / N))
FOR I = 0 TO N - 1
U0 = rR * COS(AZ + 2 * I * PI / N)
V0 = -rR * SIN(AZ + 2 * I * PI / N)
CC = 3: CALL invert
IF FC$ = "Y" THEN PAINT (CX, CY), 1, 3
NEXT I

NEXT.MENU:
LOCATE 21, 1: PRINT "PRESS ANY"
LOCATE 22, 1: PRINT "KEY TO"
LOCATE 23, 1: PRINT "CONTINUE"
U$ = INPUT$(1)
LINE (0, 0)-(75, 200), 0, BF
LOCATE 1, 1: PRINT "OPTIONS"
LOCATE 3, 1: PRINT "1 - NEW"
LOCATE 4, 1: PRINT "AZIMUTH"
LOCATE 5, 1: PRINT "FOR CHAIN"
LOCATE 7, 1: PRINT "2 - NEW"
LOCATE 8, 1: PRINT "DIAGRAM"
LOCATE 10, 1: PRINT "3 - HOW"
LOCATE 11, 1: PRINT "IT WORKS"
LOCATE 13, 1: PRINT "PRESS ANY"
LOCATE 14, 1: PRINT "OTHER KEY"
LOCATE 15, 1: PRINT "TO QUIT"
LOCATE 17, 1: PRINT "SELECT"
LOCATE 18, 1: PRINT "OPTION"
U$ = INPUT$(1)
IF U$ = "2" THEN GOTO BEGIN
IF U$ = "3" THEN
   CALL HOW.IT.WORKS
   GOTO BEGIN
   END IF
IF U$ <> "1" THEN GOTO END.OPNS

NEW.AZIMUTH:
LINE (0, 0)-(75, 200), 0, BF
PAINT (CX, CY), 0, 2
CALL master.circles
LOCATE 1, 1: PRINT "INPUT NEW"
LOCATE 2, 1: PRINT "AZIMUTH"
LOCATE 3, 1: PRINT "FOR FIRST"
LOCATE 4, 1: PRINT "CIRCLE OF"
LOCATE 5, 1: PRINT "CHAIN"
LOCATE 7, 1: INPUT AZ$: AZ = VAL(AZ$) * PI / 180
LINE (0, 0)-(75, 200), 0, BF
GOTO START.DRAW

END.OPNS:
SCREEN 0: WIDTH 80
END

SUB DRAW.CIRCLE
X0 = RAD: Y0 = 0
IF RAD > .5 THEN
  FOR T = 0 TO .52 * PI STEP .02 * PI
   X1 = RAD * COS(T): Y1 = RAD * SIN(T)
   LINE (CX + X0, CY + Y0)-(CX + X1, CY + Y1), CC
   LINE (CX - X0, CY + Y0)-(CX - X1, CY + Y1), CC
   LINE (CX + X0, CY - Y0)-(CX + X1, CY - Y1), CC
   LINE (CX - X0, CY - Y0)-(CX - X1, CY - Y1), CC
   X0 = X1: Y0 = Y1
  NEXT T
ELSE PSET (CX, CY), CC
END IF
END SUB

SUB HELP.SCREEN
SCREEN 0: WIDTH 80: CLS
COLOR 4: PRINT "                       STEINER CHAIN": PRINT
COLOR 14: PRINT "CONSTRUCT TWO CIRCLES, ONE WITHIN THE OTHER, AND THEN DRAW A CHAIN"
PRINT "OF CIRCLES IN THE SPACE BETWEEN. IN ALMOST EVERY CASE, THE CHAIN WILL NOT"
PRINT "CLOSE EXACTLY. ONLY IN SPECIAL CASES WILL THE CHAIN CLOSE EXACTLY."
PRINT : PRINT "IF THE CHAIN DOES CLOSE EXACTLY, IT IS A REMARKABLE FACT THAT IT WILL CLOSE"
COLOR 2: PRINT "REGARDLESS OF WHERE THE FIRST CIRCLE IN THE CHAIN IS DRAWN"
COLOR 14: PRINT : PRINT "SUCH A CHAIN IS CALLED A STEINER CHAIN"
PRINT "THIS PROGRAM DEMONSTRATES STEINER CHAINS"
COLOR 3: PRINT : PRINT "INPUT THE NUMBER OF CIRCLES (AT LEAST 3) IN THE CHAIN  "
COLOR 9: PRINT "    (THIS VALUE DETERMINES THE RELATIVE DIAMETERS OF THE TWO CIRCLES:"
PRINT "    THE LARGER THE VALUE OF N, THE LARGER THE INNER CIRCLE)"
COLOR 3: INPUT "INPUT THE NUMBER OF CIRCLES IN THE CHAIN  "; N$
N = VAL(N$): IF N < 3 THEN N = 3
COLOR 2: PRINT : PRINT "INPUT THE OFFSET OF THE TWO CIRCLES (0 TO 100)"
COLOR 10: PRINT "    (A VALUE OF ZERO MEANS CONCENTRIC CIRCLES,"
PRINT "    A VALUE OF 100 MEANS THE CIRCLES TOUCH)"
COLOR 2: INPUT "INPUT THE OFFSET VALUE    "; RF$
RF = VAL(RF$) / 100
IF RF < 0 THEN RF = 0
IF RF >= 1 THEN RF = .99
COLOR 5: PRINT : INPUT "INPUT THE AZIMUTH OF THE FIRST CIRCLE IN THE CHAIN "; AZ$
AZ = VAL(AZ$) * PI / 180
COLOR 12: PRINT : INPUT "DO YOU WANT TO COLOR IN THE CIRCLES - Y/N "; FC$
IF FC$ = "y" THEN FC$ = "Y"
COLOR 7: SCREEN 0: WIDTH 80: CLS
END SUB

SUB HOW.IT.WORKS
SCREEN 0: WIDTH 80: CLS
PRINT "HOW IT WORKS": PRINT
PRINT "The trick here is finding those special pairs of circles that allow"
PRINT "Steiner Chains. Devising a formula directly for the geometry of Steiner"
PRINT "Chains is brutal. Using a technique called Inversion Geometry, the problem"
PRINT "becomes rather simple": PRINT
PRINT "Inversion geometry works like this. Pick a point (the inversion point)."
PRINT "Call it I. Any other point at distance R from I is mapped to a new point"
PRINT "with the same azimuth but distance 1/R. In inversion geometry:"
PRINT "   - If I is on a line, the line inverts onto itself"
PRINT "   - If I is not on a line, the line inverts into a circle"
PRINT "   - If I is on the circumference of a circle, the circle inverts into a line"
PRINT "       (We can see the last two relations in polar coordinates. The equation"
PRINT "        of line x=1 is r=sec a. The equation of a circle with diameter "
PRINT "        (0,0)-(1,0) is r=cos a.)"
PRINT "   - If I is not on the circumference of a circle, the circle inverts"
PRINT "     into another circle.": PRINT
PRINT "Press any key to see next page": x$ = INPUT$(1)
CLS
PRINT "The fact that circles invert into circles is the key to the construction"
PRINT : PRINT "We proceed like this:"
PRINT "   1.  Decide on the number of circles in the chain"
PRINT "   2.  Draw the outer bounding circle"
PRINT "   3.  Working IN INVERSE SPACE, calculate the radii for a chain of"
PRINT "       equal circles"
PRINT "   4.  Calculate the radius of the inner circle. It will be concentric"
PRINT "       with the outer circle. It takes nothing more than simple trig"
PRINT "       to do operations 1-4."
PRINT "   5.  Decide on the offset of the bounding circles in real space"
PRINT "   6.  Find the inversion point that maps the concentric circles in"
PRINT "       inverse space into the desired offset circles in real space."
PRINT "   7.  Calculate the locations of the bounding circles and the chain"
PRINT "       circles in inverse space"
PRINT "   8.  Invert the locations through the inversion point to find the"
PRINT "       coordinates in real space.": PRINT
PRINT "Press any key to continue": x$ = INPUT$(1)
END SUB

SUB invert
REM INVERT A CIRCLE
REM ORIGINAL CIRCLE AT U0,V0, RADIUS RO
REM INVERSION POINT AT P,Q
REM INVERTED CIRCLE AT I0,J0, RADIUS RI
REM COMPUTE DISTANCE U0,V0 - P,Q
L = SQR((U0 - P) ^ 2 + (V0 - Q) ^ 2)
IF L = 0 THEN
     REM INVERT CIRCLE CENTERED ON INVERSION POINT
     I0 = P: J0 = Q: RI = 1 / RO
     CX = I0 * SF + DX: CY = J1 * SF + DY: RAD = RI * SF
     CALL DRAW.CIRCLE
ELSE
     REM INVERT CIRCLE NOT CENTERED ON INVERSION POINT
     REM INVERT NEAR SIDE OF CIRCLE
     U1 = U0 + RO * (P - U0) / L: V1 = V0 + RO * (Q - V0) / L
     R = SQR((U1 - P) ^ 2 + (V1 - Q) ^ 2)
     I1 = P + ((U1 - P) / R) / R: J1 = Q + ((V1 - Q) / R) / R
     REM INVERT FAR SIDE OF CIRCLE
     U2 = U0 - RO * (P - U0) / L: V2 = V0 - RO * (Q - V0) / L
     R = SQR((U2 - P) ^ 2 + (V2 - Q) ^ 2)
     I2 = P + ((U2 - P) / R) / R: J2 = Q + ((V2 - Q) / R) / R
     REM PLOT INVERSE CIRCLE
     I0 = (I1 + I2) / 2: J0 = (J1 + J2) / 2
     RI = SQR((I1 - I2) ^ 2 + (J1 - J2) ^ 2) / 2
     CX = I0 * SF + DX: CY = J0 * SF + DY: RAD = RI * SF
     CALL DRAW.CIRCLE
END IF
END SUB

SUB master.circles
U0 = 0: V0 = 0: RO = (1 - SIN(PI / N)) / (1 + SIN(PI / N))
R = SQR(P * P + Q * Q): SF = RC * (RO * RO - R * R) / RO
DX = XC - P * SF / (RO * RO - R * R)
DY = YC - Q * SF / (RO * RO - R * R)
CC = 2: CALL invert
U0 = 0: V0 = 0: RO = 1: CC = 2: CALL invert
IF FC$ = "Y" THEN PAINT (CX, CY), 3, 2
END SUB