Статья Довольно странный, но действенный пример интеграции физики в Python.

[neopaket]

Предисловие от @neopaket.

Приветствую, недавно начал углублёно изучать модули для Python, связанные с физическими и инженерными расчётами (Например SciPy, NumPy, SumPy). Представьте, вы решили построить для каких то целей арку-мост из дерева (ну или же любого другого материала). Что вам для этого потребуется?: лист бумаги (Разлинованный (Конечно если не будете писать поперёк ) ), пишущие принадлежности, линейка, циркуль, транспортир и всё, что только душе угодно. Нет! Не обязательно. Конечно есть много программ для 3D-моделирования, но это можно сделать при помощи простой среди разработки и сейчас мы узнаем как!

Джон Тапселл.

Я хотел сделать мост из деревянных кирпичей, без клея и соединителей, удерживаемых только гравитацией.
Мост представляет собой перевернутую канатную сеть, которую вы можете увидеть в любом парке.

Что мы можем узнать из простейшего физического уравнения:

y = a * cosh(x / a)

Где: cosh - гиперболическая функция. В России cosh трактуется в виде ch.

Как же перевести это уравнение в IPython? А, ну да, вот:

L = 300.0 #300 миллиметров ширина моста. Все размеры будут произведены в миллиметрах.
A = 100.0 #Масштабный фактор  
 
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()

В итоге код производит:

Выглядит круто, не правда ли?

Теперь, я хочу превратить его из зародышевого состояния в деревянный мост. На самом деле он не будет выглядеть как кривая. Представьте, что вы строите мост из ступеней в майнкрафт

[TABLE]
[TR]
[TD]class 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[/TD]
[/TR]
[/TABLE]

И вывод:

Чертёж почти закончен, стоит заняться и подготовкой материалов для мостика.
Теперь нам нужно решить, как же правильно рубить дрова. Для того, чтобы создать арку, нужно рассчитать угол каждого куска дерева (Линии прикосновения), по отношению к горизонтали:

. Таким образом, угол между двумя смежными кусками дерева (Вспомним квадрат. Смежные линии это две стороны, которые в точке соприкосновения с внутренней части угла образуют угол в 90*) - это разница между их углами.

Древесина будет параллельна этим касательным линиям, из за чего мы можем нарисовать нормальные линии для этих касательных линий толщины дерева (Слишком много линий):

def dy(x):
    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 #В миллиметрах
 
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

Выводик:

Красная кривая представляет силовую линию (Пока не понятно). Пакетик расскажет вам. В центрах силовых линий находятся очень опасные места. Если на них наступит тяжёлый человек, то она деформируется или мост вообще сломается.

Таким образом, чтобы минимизировать силу, которая давит на плиты, где встречается древесина, нам нужно, чтобы срезы дерева были нормальны к красной линии. Мы кое что добавим.

def getWoodCorners(block, x):
    '''Возвращаем вершину и низ блока дерева, разрез которого проходит через x, y (x), а середина которого находится в 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)
 
#Нарисуем все линии
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

И как обычно:

Здесь следует кое что отметить:

• Код, видимо, выполняет некоторые довольно избыточные вычисления при расчётах, но большая часть этих излишек состоит в том, чтобы сделать правильными знаки углов. Проще позволить atan2 обрабатывать знаки для квадрантов, чем пытаться упростить код и справиться с этим самостоятельно.
• Блоки на самом деле не одинакового размера в местах их встречи. Различия составляют порядка 1% и, однако, на этих чертежах не заметны.

Теперь нам нужно нарисовать и пометить эти блоки так, чтобы их было легче вырезать:

fig, ax = subplots(1,1)
fig.set_size_inches(17,7)
axis('off')
 
def rotatePolygon(polygon,theta):    
    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):
    """смещает многоугольник на заданные координаты"""
    return [ (x+xshift, y+yshift) for (x,y) in polygon]
 
sawThickness = 3.0 #Добавить толщину 3мм между блоками 
#Draw all the blocks
xshift = 10.0 #Начните с 10 мм для того, чтобы легче было резать вручную
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

Конечный вывод:

Top coordinates: [10, 141, 144, 228, 231, 307, 310, 394, 397, 528]
Bottom Coordinates: [43, 131, 156, 214, 247, 291, 324, 382, 407, 495

Послесловие от меня.

Я впервые перевожу материал как то связанный с инженерией и физикой. Если вам понравилось, то был бы рад вашему комментарию. И помните - "погрешность = необходимость".

Спасибо за прочтение данной статьи и хорошего дня.
Оригинальная статья: https://johnflux.com/category/bridge/
neopaket