I wanted to make a bridge out of wooden bricks, with no glue or connectors, held up by gravity only.
The bridge is an inverted catenary:
Which follows the equation:
Now, turning this into IPython code:
L = 300.0 # 300mm wide bridge. All lengths in millimetres
A = 100.0 # Scaling factor
def y(x):
return -A*( math.cosh(x/A) - math.cosh((L/2)/A))
X = linspace(-L/2, L/2)
Y = [y(x) for x in X]
fig, ax = subplots(1,1)
fig.set_size_inches(14,7)
ax.plot(X,Y,'r')
xlabel('x')
ylabel('height')
title('Bridge curve')
show()
Which produces:
This is looking good.
Now, I want to build this out of wood, so we won’t actually have this curve, but a ‘blocky’ version of this.
class Block: '''A wooden block''' def __init__(self, start_x, end_x): self.start_x = start_x self.start_y = y(start_x) self.end_x = end_x self.end_y = y(end_x) self.length_x = end_x - start_x self.midpoint_x = start_x + self.length_x/2.0 self.tangent = (self.end_y - self.start_y) / self.length_x self.midpoint_y = self.start_y + (self.end_y - self.start_y)/2.0 self.coords = [] # filled in later def angle(self): return math.atan(self.tangent) step = 60.0 def plotTangent(block): xstart = block.midpoint_x-block.length_x/2 xend = block.midpoint_x+block.length_x/2 ystart = block.midpoint_y + block.tangent * (xstart - block.midpoint_x) yend = block.midpoint_y + block.tangent * (xend - block.midpoint_x) ax.plot([block.start_x, block.end_x], [block.start_y, block.end_y], 'k-', lw=1) blocks = [] for start_x in linspace(-L/2, L/2, L/step, endpoint=False): block = Block(start_x, start_x+step) blocks.append(block) plotTangent(block) fig
Which produces:
Now we have to decide how to cut the wood. To produce the arch, the angle of each piece of wood (the tangent lines), relative to the horizontal, is just: 
The wood will be parallel with those tangent lines, so we can draw on the normal lines for those tangent lines, of the thickness of the wood:
def dy(x): '''return an approximation of dy/dx''' return (y(x+0.05) - y(x-0.05)) / 0.1 def normalAngle(x): if x == -L/2: return 0 elif x == L/2: return math.pi return math.atan2(-1.0, dy(x)); WoodThickness = 50.0 # in millimeters def getNormalCoords(x, ymidpoint, lineLength): tangent = dy(x) if tangent != 0: if x == -L/2 or x == L/2: normalSlope = 0 else: normalSlope = -1 / dy(x) xlength = lineLength * math.cos(normalAngle(x)) xstart = x-xlength/2 xend = x+xlength/2 ystart = ymidpoint + normalSlope * (xstart - x) yend = ymidpoint + normalSlope * (xend - x) else: xstart = x xend = x ystart = ymidpoint + lineLength/2 yend = ymidpoint - lineLength/2 return (xstart, ystart, xend, yend) for block in blocks: (xstart, ystart, xend, yend) = getNormalCoords(block.midpoint_x, block.midpoint_y, WoodThickness) ax.plot([xstart, xend], [ystart, yend], 'b-', lw=1) fig
Now, the red curve is going to represent the line of force. My reasoning is that if there was any force normal to the red curve, then the string/bridge/weight would thus accelerate and change position.
So to minimize the lateral forces where the wood meets, we want the wood cuts to be normal to the red line. We add that, then draw the whole block:
def getWoodCorners(block, x): ''' Return the top and bottom of the piece of wood whose cut passes through x,y(x), and whose middle is at xmid ''' adjusted_thickness = WoodThickness / math.cos(normalAngle(x) - normalAngle(block.midpoint_x)) return getNormalCoords(x, y(x), adjusted_thickness) def drawPolygon(coords, linewidth): xcoords,ycoords = zip(*coords) xcoords = xcoords + (xcoords[0],) ycoords = ycoords + (ycoords[0],) ax.plot(xcoords, ycoords, 'k-', lw=linewidth) #Draw all the blocks for block in blocks: (xstart0, ystart0, xend0, yend0) = getWoodCorners(block, block.start_x) # Left side of block (xstart1, ystart1, xend1, yend1) = getWoodCorners(block, block.end_x) # Right side of block block.coords = [ (xstart0, ystart0), (xstart1, ystart1), (xend1, yend1), (xend0, yend0) ] drawPolygon(block.coords, 2) # Draw block fig
There’s a few things to note here:
- The code appears to do some pretty redundant calculations when calculating the angles – but much of this is to make the signs of the angles correct. It’s easier and nicer to let atan2 handle the signs for the quadrants, than to try to simplify the code and handle this ourselves.
- The blocks aren’t actually exactly the same size at the points where they meet. The differences are on the order of 1%, and not noticeable in these drawings however.
Now we need to draw and label these blocks in a way that makes it easiest to cut:
fig, ax = subplots(1,1)
fig.set_size_inches(17,7)
axis('off')
def rotatePolygon(polygon,theta):
"""Rotates the given polygon which consists of corners represented as (x,y), around the ORIGIN, clock-wise, theta radians"""
return [(x*math.cos(theta)-y*math.sin(theta) , x*math.sin(theta)+y*math.cos(theta)) for (x,y) in polygon]
def translatePolygon(polygon, xshift,yshift):
"""Shifts the polygon by the given amount"""
return [ (x+xshift, y+yshift) for (x,y) in polygon]
sawThickness = 3.0 # add 3 mm gap between blocks for saw thickess
#Draw all the blocks
xshift = 10.0 # Start at 10mm to make it easier to cut by hand
topCutCoords = []
bottomCutCoords = []
for block in blocks:
coords = translatePolygon(block.coords, -block.midpoint_x, -block.midpoint_y)
coords = rotatePolygon(coords, -block.angle())
coords = translatePolygon(coords, xshift - coords[0][0], -coords[3][1])
xshift = coords[1][0] + sawThickness
drawPolygon(coords,1)
(topLeft, topRight, bottomRight, bottomLeft) = coords
itopLeft = int(round(topLeft[0]))
itopRight = int(round(topRight[0]))
ibottomLeft = int(round(bottomLeft[0]))
ibottomRight = int(round(bottomRight[0]))
topCutCoords.append(itopLeft)
topCutCoords.append(itopRight)
bottomCutCoords.append(ibottomLeft)
bottomCutCoords.append(ibottomRight)
ax.text(topLeft[0], topLeft[1], itopLeft)
ax.text(topRight[0], topRight[1], itopRight, horizontalalignment='right')
ax.text(bottomLeft[0], bottomLeft[1], ibottomLeft, verticalalignment='top', horizontalalignment='center')
ax.text(bottomRight[0], bottomRight[1], ibottomRight, verticalalignment='top', horizontalalignment='center')
print "Top coordinates:", topCutCoords
print "Bottom Coordinates:", bottomCutCoords
Which produces:
Top coordinates: [10, 141, 144, 228, 231, 307, 310, 394, 397, 528] Bottom Coordinates: [43, 131, 156, 214, 247, 291, 324, 382, 407, 495]
Now to cut this!
John Tapsell 