Given A=(a,0), B=(0,0), and C=(0,a)
Let f(a)=the total number of unit equilateral triangles XYZ that can be formed such that the lengths AX, BY, and CZ are all 1 unit.
Give a piecewise definition by intervals for f(a)
0.000 0
0.000 4
0.000 4
0.001 0
0.491 0
0.492 4 -- intermediate with 2 ??
0.517 4
0.518 6
0.635 6
0.636 8
0.703 8
0.704 12
0.729 12
0.730 11
0.738 11
0.739 12 -- one or another of two points added to the 11 ?
0.739 12
0.740 11
0.740 11
0.741 12
1.065 12
1.066 8
1.152 8
1.153 6
1.414 6
1.415 2
1.931 2
1.932 4
1.955 4
1.956 3
1.956 3 -- tangency ?
1.957 4
1.957 4
1.958 0
The vacillation around a=0 can be attributed to the degenerate nature of that case, so the summary will void the zero case.
The comments above, ask questions about whether certain conditions hold in transition areas as to whether one or two points are added at a time somewhere between the chosen points at .001 intervals.
Summarized:
.000 - .491 0
.492 - .517 4
.518 - .635 6
.636 - .703 8
.704 - .729 12
.730 - .738 11
.739 12
.740 11
.741 - 1.065 12
1.066 - 1.152 8
1.153 - 1.414 6
1.415 - 1.931 2
1.932 - 1.955 4
1.956 3
1.957 4
1.958 - 0
again, with the possibility that transitional numbers might exist between the end of one range and the beginning of the next, as there's a .001 gap.
The program was based upon my experimentation with Geometer's Sketchpad, and that product was also used to spot check the answers given by the program.
For any given multiple of .001 for a, points at 1/50 of a degree on a unit circle about C were used for Z. Then the intersections of a unit circle about Z with the unit circle about A were used as two values for the location of X. Sixty degrees in either direction about the unit circle around Z were then checked to see if they were inside, outside or directly on the unit circle about B; any time that changed inside/outside or hit exactly on the circumference, it was counted. One tricky part was keeping the identities of the two X points and potential Y points preserved upon iteration of the location of Z, as a transition from in to out or vice versa had to refer to the same moving point, rather than an identity shift dependent on the solution to the simultaneous equations (quadratic) for the circle intersection. Another tricky spot is the assurance that the circle that Z made about C was completely closed exactly, as the increment of its angular rotation, .02 degrees is not exact when represented in binary, as it is, internally.
The code below includes commented out debugging code (with leading apostrophes).
DECLARE SUB quadSolve (a#, b#, c#, sol1#, sol2#, flag#)
DECLARE SUB pol2rect (rho#, theta#, x#, y#)
DECLARE SUB rect2pol (x#, y#, rho#, theta#)
DEFDBL A-Z
DIM SHARED pi
pi = 4 * ATN(1)
dr = pi / 180
a = 1.9576
zx = .862: zy = 1.4507
' for linear eq. y=ax + b for intersections of circles about A and Z:
linA = (a - zx) / zy
linB = (zy * zy + zx * zx - a * a) / (2 * zy)
quadA = 1 + linA * linA
quadB = -2 * a + 2 * linA * linB
quadC = a * a + linB * linB - 1
quadSolve quadA, quadB, quadC, x1, x2, solns
PRINT solns
y1 = linA * x1 + linB: y2 = linA * x2 + linB
PRINT USING "####.#######"; x1; y1; x2; y2
CLS
FOR a = 0# TO 2.5# STEP .001#
FOR Ypt = 1 TO 4: hitCt(Ypt) = 0: NEXT
stepper = .02#
FOR zAngle = 0 TO 360 STEP stepper
IF zAngle > 360 - stepper / 2 THEN
pol2rect 1, 0, zx, zy
ELSE
pol2rect 1, zAngle * dr, zx, zy
END IF
zy = zy + a ' translate ctr of coords
' for linear eq. y=ax + b for intersections of circles about A and Z:
IF zy <> 0 THEN
linA = (a - zx) / zy
linB = (zy * zy + zx * zx - a * a) / (2 * zy)
quadA = 1 + linA * linA
quadB = -2 * a + 2 * linA * linB
quadC = a * a + linB * linB - 1
quadSolve quadA, quadB, quadC, x1, x2, solns
ELSE
x1 = (zx + 1) / 2: x2 = x1
solns = -1
END IF
IF solns THEN
IF x2 <> x1 AND ABS(linA) < 10000 AND ABS(linB) < 10000 THEN
y1 = linA * x1 + linB: y2 = linA * x2 + linB
ELSE
y1 = SQR(1 - x1 * x1)
y2 = -SQR(1 - x2 * x2)
END IF
rect2pol x1 - zx, y1 - zy, r1, th1
rect2pol x2 - zx, y2 - zy, r2, th2
diffnce = (th2 - th1) / (2 * pi)
diffnce = diffnce - INT(diffnce): IF diffnce > .5 THEN diffnce = diffnce - 1
IF diffnce > 0 THEN
SWAP x1, x2: SWAP y1, y2
END IF
rect2pol x1 - zx, y1 - zy, r, th
pol2rect r, th - 60 * dr, yx(1), yy(1)
yx(1) = yx(1) + zx: yy(1) = yy(1) + zy
pol2rect r, th + 60 * dr, yx(2), yy(2)
yx(2) = yx(2) + zx: yy(2) = yy(2) + zy
rect2pol x2 - zx, y2 - zy, r, th
pol2rect r, th - 60 * dr, yx(3), yy(3)
yx(3) = yx(3) + zx: yy(3) = yy(3) + zy
pol2rect r, th + 60 * dr, yx(4), yy(4)
yx(4) = yx(4) + zx: yy(4) = yy(4) + zy
FOR Ypt = 1 TO 4
cDist = 1 - SQR(yx(Ypt) * yx(Ypt) + yy(Ypt) * yy(Ypt))
distPrev(Ypt) = distnce(Ypt)
distnce(Ypt) = cDist
IF zAngle > 0 AND zAngle - AngPrev < 2 * stepper THEN
IF (distPrev(Ypt) * distnce(Ypt) < 0 AND prevSolns) OR distnce(Ypt) = 0 THEN
hitCt(Ypt) = hitCt(Ypt) + 1
' PRINT USING "###.###"; zAngle;
IF Ypt <= 2 THEN
' PRINT USING "###.####"; x1; y1;
ELSE
' PRINT USING "###.####"; x2; y2;
END IF
' PRINT SPACE$(10);
' PRINT USING "###.####"; yx(Ypt); yy(Ypt);
' PRINT SPACE$(10);
' PRINT USING "###.####"; zx; zy
' PRINT USING "###.###"; x1prev(Ypt); y1prev(Ypt); x2prev(Ypt); y2prev(Ypt)
' PRINT USING "###.###"; x1; y1; x2; y2
END IF
END IF
x1prev(Ypt) = x1: x2prev(Ypt) = x2
y1prev(Ypt) = y1: y2prev(Ypt) = y2
NEXT
AngPrev = zAngle
END IF
prevSolns = solns
NEXT zAngle
totHit = 0
FOR i = 1 TO 4
' PRINT hitCt(i)
totHit = totHit + hitCt(i)
NEXT
IF prevTot <> totHit THEN
PRINT USING "##.### ###"; prevA; prevTot
PRINT USING "##.### ###"; a; totHit
END IF
prevTot = totHit: prevA = a
NEXT a
SUB pol2rect (rho, theta, x, y)
x = rho * COS(theta)
y = rho * SIN(theta)
END SUB
SUB quadSolve (a, b, c, sol1, sol2, flag)
discrim = b * b - 4 * a * c
IF discrim < 0 THEN flag = 0: EXIT SUB
diff = SQR(discrim)
sol1 = (-b + diff) / (2 * a)
sol2 = (-b - diff) / (2 * a)
flag = -1
END SUB
SUB rect2pol (x, y, rho, theta)
IF x > 0 THEN
theta = ATN(y / x)
rho = SQR(x * x + y * y)
ELSEIF x < 0 THEN
theta = ATN(y / x) + pi
rho = SQR(x * x + y * y)
ELSE
theta = pi / 2 * SGN(y)
rho = ABS(y)
END IF
END SUB
|
|
Posted by Charlie
on 2009-09-14 23:17:34 |
computer exploration
