Welcome, guest | Sign In | My Account | Store | Cart

It creates fractal-like map plots from the simulation.

(See my other post titled "Dynamical Billiards Simulation" first!)

I had to keep image size and maxSteps small otherwise the calculation takes too long!

(It shows what percentage of calculations completed every 10 seconds also.)

Python, 140 lines
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
# Dynamical Billiards Simulation Map (Fractal?)
# FB - 201011077
# 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)

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')

2 comments

FB36 (author) 11 years, 1 month ago  # | flag

Mapping algorithm is this:

Select a random initial direction for the ball.

Start the ball moving from the first pixel of the image (0,0) using that initial direction.

Color the first pixel of the image according to the final location coordinates of the ball after moving maxSteps.

Next time start the ball moving from the next image pixel using the same initial direction and so on.

FB36 (author) 10 years, 11 months ago  # | flag

This version uses n by m grid of circles:

import math
import random
from PIL import Image,ImageDraw
imgx=300
imgy=200
image=Image.new("RGB",(imgx,imgy))
draw=ImageDraw.Draw(image)
coloring=random.randint(0,7)
maxSteps=200
a=2.0*math.pi*random.random()
n=random.randint(1,9)
m=random.randint(1,9)
crMax=int(min(imgx/(n+1)/2,imgy/(m+1)/2))-1
crMin=10
cr=random.randint(crMin,crMax)
cxList=[]
cyList=[]
crList=[]
for j in range(m):
    cy=int((j+1)*imgy/(m+1))
    for i in range(n):
        cx=int((i+1)*imgx/(n+1))
        cxList.append(cx)
        cyList.append(cy)
        crList.append(cr)

for y0 in range(imgy):
    for x0 in range(imgx):
        x=float(x0)
        y=float(y0)
        s=math.sin(a)
        c=math.cos(a)
        flag=True
        for i in range(n*m):
            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
                if xnew<0 or xnew>imgx-1:
                    c=-c
                    xnew=x
                if ynew<0 or ynew>imgy-1:
                    s=-s
                    ynew=y
                for i in range(n*m):
                    if math.hypot(xnew-cxList[i],ynew-cyList[i])<=crList[i]:
                        ca=math.atan2(ynew-cyList[i],xnew-cxList[i])
                        rca=math.atan2(-s,-c)
                        rab=rca+(ca-rca)*2
                        s=math.sin(rab)
                        c=math.cos(rab)
                        xnew=x
                        ynew=y
                x=xnew
                y=ynew
            if coloring==0:
                d=math.hypot(x,y)/math.hypot(imgx-1,imgy-1)
            elif coloring==1:
                d=math.hypot(x-x0,y-y0)/math.hypot(imgx-1,imgy-1)
            elif coloring==2:
                d=(x+y)/(imgx+imgy-2)
            elif coloring==3:
                d=x*y/((imgx-1)*(imgy-1))
            elif coloring==4:
                d=x/(imgx-1)
            elif coloring==5:
                d=y/(imgy-1)
            elif coloring==6:
                ang=math.atan2(y-(imgy-1)/2,x-(imgx-1)/2)
            elif coloring==7:
                ang=math.atan2(y-y0,x-x0)
            if coloring>=6:
                if ang < 0: ang=2*math.pi-math.fabs(ang)
                d=ang/2*math.pi
            k=int(d*255)
            rd=k%8*32
            gr=k%16*16
            bl=k%32*16
            image.putpixel((x0,y0),(rd,gr,bl))
image.save('Dynamical_Billiards_Map_color'+str(coloring)+'.png','PNG')