#! /usr/bin/python
# Relaxation simulation of the temperature at the center of the pentagon.
# Four edges are held at 10 degrees, and the fifth at 80 degrees.
from scipy import *
# rounds to nearest integer, and returns an integer
def intround(f):
return int(round(f))
# Bresenham's algorithm: returns list of lattice points on the line connecting
# r1 and r2
def line(r1, r2):
x0, y0 = intround(r1[0]), intround(r1[1])
x1, y1 = intround(r2[0]), intround(r2[1])
points = []
steep = abs(y1 - y0) > abs(x1 - x0)
if steep:
x0,y0 = y0,x0
x1,y1 = y1,x1
if x0 > x1:
x0,x1=x1,x0
y0,y1=y1,y0
deltax = x1 - x0
deltay = abs(y1 - y0)
error = 0
deltaerr = float(deltay) / deltax
y = y0
if y0 < y1:
ystep = 1
else:
ystep = -1
for x in range(x0,x1+1):
if steep:
points.append((y,x))
else:
points.append((x,y))
error += deltaerr
if error >= 0.5:
y += ystep
error -= 1.0
return points
def complex2pair(c):
return (real(c),imag(c))
def set_edge_temps(grid):
for e in lo_temp_edges:
grid[e[0]][e[1]] = 10
for e in hi_temp_edges:
if e in lo_temp_edges:
grid[e[0]][e[1]] = 45 # corner joining high- and lo-temp edges
else:
grid[e[0]][e[1]] = 80
angle = 72.0/180*pi # 72 degrees
# use complex plane to find the vertices
r = exp(angle*1j)
# pentagon vertices in the complex plane, with first vertex duplicated
# at the end of the list
corners = array([r**(i+0.25) for i in range(6)])
# translate pentagon into first quadrant
corners -= complex(min(real(corners)), min(imag(corners)))
corners *= 50 # grid spacing (bigger means finer spacing)
center = sum(corners[0:5])/5 # center of pentagon
# use Bresenham's algorithm to find the lattice points on the edges
lo_temp_edges = []
hi_temp_edges = []
for i in range(4):
lo_temp_edges += line(complex2pair(corners[i]), complex2pair(corners[i+1]))
hi_temp_edges = line(complex2pair(corners[4]), complex2pair(corners[5]))
# figure out the grid dimensions
max_x = max([r[0] for r in lo_temp_edges+hi_temp_edges])
max_y = max([r[1] for r in lo_temp_edges+hi_temp_edges])
grid = zeros((max_x+1,max_y+1))
dirs = [(-1,0), (1,0), (0,1), (0,-1)]
while True:
newgrid = zeros((max_x+1,max_y+1))
set_edge_temps(grid) # impose constraint
for x in range(max_x+1): # relax each location to avg of its neighbors
for y in range(max_y+1):
total = 0.0
n = 0
for d in dirs: # use each neighbor that's within the grid
try:
total += grid[x+d[0]][y+d[1]]
n += 1
except: # that neighbor was not inside the grid
pass
newgrid[x][y] = total/n # but save new value in a new grid
grid, newgrid = newgrid, grid # swap new and old grid
# print temperature at the center of the pentagon
print grid[intround(real(center))][intround(imag(center))]