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