Welcome, guest | Sign In | My Account | Store | Cart
```# Dynamical Billiards Simulation Map (Fractal?)
# FB - 201011077
# http://en.wikipedia.org/wiki/Dynamical_billiards
# http://www.scholarpedia.org/article/Dynamical_billiards
import math
import random
import time
from PIL import Image, ImageDraw
imgx = 300
imgy = 200
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)

coloring = random.randint(0, 7) # choose a coloring method
print 'Using the coloring method: ' + str(coloring)

maxSteps = 200 # of steps of ball motion (in constant speed)

n = random.randint(1, 7) # of circular obstacles
crMax = int(min(imgx - 1, imgy - 1) / 4) # max circle radius
crMin = 10 # min circle radius

# create circular obstacle(s)
cxList = []
cyList = []
crList = []
for i in range(n):
while(True): # circle(s) must not overlap
cr = random.randint(crMin, crMax) # circle radius
cx = random.randint(cr, imgx - 1 - cr) # circle center x
cy = random.randint(cr, imgy - 1 - cr) # circle center y
flag = True
if i > 0:
for j in range(i):
if math.hypot(cx - cxList[j], cy - cyList[j]) < cr + crList[j]:
flag = False
break
if flag == True:
break
draw.ellipse((cx - cr, cy - cr, cx + cr, cy + cr))
cxList.append(cx)
cyList.append(cy)
crList.append(cr)

# initial direction of the ball
a = 2.0 * math.pi * random.random()

t0 = time.time()

for y0 in range(imgy):
for x0 in range(imgx):
x = float(x0)
y = float(y0)
s = math.sin(a)
c = math.cos(a)

# print '%completed' every 10 seconds
t = time.time()
if t - t0 >= 10:
print '%' + str(int(100 * (imgx * y0 + x0) / (imgx * imgy)))
t0 = t

# initial location of the ball must be outside of the circle(s)
flag = True
for i in range(n):
if math.hypot(x - cxList[i], y - cyList[i]) <= crList[i]:
flag = False
break

if flag:
for i in range(maxSteps):
xnew = x + c
ynew = y + s

# reflection from the walls
if xnew < 0 or xnew > imgx - 1:
c = -c
xnew = x
if ynew < 0 or ynew > imgy - 1:
s = -s
ynew = y

# reflection from the circle(s)
for i in range(n):
if math.hypot(xnew - cxList[i], ynew - cyList[i]) <= crList[i]:
# angle of the circle point
ca = math.atan2(ynew - cyList[i], xnew - cxList[i])
# reversed collision angle of the ball
rca = math.atan2(-s, -c)
# reflection angle of the ball
rab = rca + (ca - rca) * 2
s = math.sin(rab)
c = math.cos(rab)
xnew = x
ynew = y

x = xnew
y = ynew

# Color the starting point according to the final point!
# The color can be decided in many different ways.
# Only 8 methods implemented here.
if coloring == 0:
# absolute distance method
d = math.hypot(x, y) / math.hypot(imgx - 1, imgy - 1)
elif coloring == 1:
# relative distance method
d = math.hypot(x - x0, y - y0) / math.hypot(imgx - 1, imgy - 1)
elif coloring == 2:
# x+y method
d = (x + y) / (imgx + imgy - 2)
elif coloring == 3:
# x*y method
d = x * y / ((imgx - 1) * (imgy - 1))
elif coloring == 4:
# x-coordinate method
d = x / (imgx - 1)
elif coloring == 5:
# y-coordinate method
d = y / (imgy - 1)
elif coloring == 6:
# absolute angle method by taking the image center as origin
ang = math.atan2(y - (imgy - 1) / 2, x - (imgx - 1) / 2)
elif coloring == 7:
# relative angle method by taking the starting point as origin
ang = math.atan2(y - y0, x - x0)
if coloring >= 6:
# convert the angle from -pi..pi to 0..2pi
if ang < 0: ang = 2 * math.pi - math.fabs(ang)
d = ang / 2 * math.pi # convert the angle to 0..1

k = int(d * 255)
rd = k % 8 * 32
gr = k % 16 * 16
bl = k % 32 * 16
image.putpixel((x0, y0), (rd, gr, bl))

print 'Calculations completed.'
image.save('Dynamical_Billiards_Map_color' + str(coloring) + '.png', 'PNG')
```

#### Diff to Previous Revision

```--- revision 3 2010-11-06 20:35:51
+++ revision 4 2010-11-07 18:21:06
@@ -1,5 +1,5 @@
# Dynamical Billiards Simulation Map (Fractal?)
-# FB - 201011066
+# FB - 201011077
# http://en.wikipedia.org/wiki/Dynamical_billiards
# http://www.scholarpedia.org/article/Dynamical_billiards
@@ -12,7 +12,7 @@
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)

-coloring = random.randint(0, 5) # choose a coloring method
+coloring = random.randint(0, 7) # choose a coloring method
print 'Using the coloring method: ' + str(coloring)

maxSteps = 200 # of steps of ball motion (in constant speed)
@@ -100,26 +100,32 @@

# Color the starting point according to the final point!
# The color can be decided in many different ways.
-            # Only 6 methods implemented here.
+            # Only 8 methods implemented here.
if coloring == 0:
+                # absolute distance method
+                d = math.hypot(x, y) / math.hypot(imgx - 1, imgy - 1)
+            elif coloring == 1:
+                # relative distance method
+                d = math.hypot(x - x0, y - y0) / math.hypot(imgx - 1, imgy - 1)
+            elif coloring == 2:
+                # x+y method
+                d = (x + y) / (imgx + imgy - 2)
+            elif coloring == 3:
+                # x*y method
+                d = x * y / ((imgx - 1) * (imgy - 1))
+            elif coloring == 4:
# x-coordinate method
d = x / (imgx - 1)
-            elif coloring == 1:
+            elif coloring == 5:
# y-coordinate method
d = y / (imgy - 1)
-            elif coloring == 2:
-                # absolute distance method
-                d = math.hypot(x, y) / math.hypot(imgx - 1, imgy - 1)
-            elif coloring == 3:
-                # relative distance method
-                d = math.hypot(x - x0, y - y0) / math.hypot(imgx - 1, imgy - 1)
-            elif coloring == 4:
+            elif coloring == 6:
# absolute angle method by taking the image center as origin
-                ang = math.atan2(x - (imgx - 1) / 2, y - (imgy - 1) / 2)
-            elif coloring == 5:
+                ang = math.atan2(y - (imgy - 1) / 2, x - (imgx - 1) / 2)
+            elif coloring == 7:
# relative angle method by taking the starting point as origin
-                ang = math.atan2(x - x0, y - y0)
-            if coloring >= 4:
+                ang = math.atan2(y - y0, x - x0)
+            if coloring >= 6:
# convert the angle from -pi..pi to 0..2pi
if ang < 0: ang = 2 * math.pi - math.fabs(ang)
d = ang / 2 * math.pi # convert the angle to 0..1
```