#! /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))]