Welcome, guest | Sign In | My Account | Store | Cart
# Dynamical Billiards Simulation Map (Fractal?)
# FB - 201011066
# More info:
# 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)

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

            # Plot the starting point according to the final point!
            # The color can be decided in many different ways:
            # absolute distance method
            d = math.hypot(x, y) / math.hypot(imgx - 1, imgy - 1)
            # relative distance method
            # d = math.hypot(x - x0, y - y0) / math.hypot(imgx - 1, imgy - 1)
            # absolute angle method by taking the image center as origin
            # ang = math.atan2(x - (imgx - 1) / 2, y - (imgy - 1) / 2)
            # relative angle method by taking the starting point as origin
            # ang = math.atan2(x - x0, y - y0)
            # 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.png", "PNG")

Diff to Previous Revision

--- revision 1 2010-11-06 07:06:23
+++ revision 2 2010-11-06 07:17:11
@@ -42,8 +42,6 @@
     
 # initial direction of the ball
 a = 2.0 * math.pi * random.random()
-s = math.sin(a)
-c = math.cos(a)
 
 t0 = time.time()
 
@@ -51,6 +49,8 @@
     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()
@@ -102,12 +102,12 @@
             # relative distance method
             # d = math.hypot(x - x0, y - y0) / math.hypot(imgx - 1, imgy - 1)
             # absolute angle method by taking the image center as origin
-            # a = math.atan2(x - (imgx - 1) / 2, y - (imgy - 1) / 2)
+            # ang = math.atan2(x - (imgx - 1) / 2, y - (imgy - 1) / 2)
             # relative angle method by taking the starting point as origin
-            # a = math.atan2(x - x0, y - y0)
+            # ang = math.atan2(x - x0, y - y0)
             # convert the angle from -pi..pi to 0..2pi
-            # if a < 0: a = 2 * math.pi - math.fabs(a)
-            # d = a / 2 * math.pi # convert the angle to 0..1
+            # 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

History