파이썬 예제 (저항력, 중력)
Drag force and Gravity
이번에는 Python, PyQt5를 이용해 구현한 저항력(Drag force), 중력(Gravity) 예제를 소개하려 합니다.자세한 소스코드는 맨 뒤에 살펴보고, 일단 수식 유도과정과 결과만 살펴보겠습니다.
앞 게시물인 벡터(1편), 힘과 가속도(2편)에 이어지는 내용입니다.
1. 유체 저항
어떤 물체가 공기나 물속에서 받는 저항을 어떻게 프로그래밍으로 구현할 수 있을까요?
아래는 유체 저항(벡터, 힘)을 구하는 수식입니다.
각각, F 는 저항력, C는 저항 계수, A는 유체와 닿는 앞면적, p는 액체밀도, V는 속도의 방향을 뜻하는 단위벡터입니다.
저항력 F 벡터를 구해 속도의 방향과 반대로 적용해 주면 힘이 감소됩니다.
|
| [유체저항 공식] |
앞 게시물에서 가속도는 속도에 영향을 미치고, 속도는 위치에 영향을 미친다는 사실을 살펴보았습니다.
또한 가속도, 속도, 위치는 벡터라는 것도 살펴보았습니다.
실제 우리가 사는 세계나 정확한 계산이 필요한 시뮬레이션 소프트웨어라면 위의 수식에 들어가는 계수(Coefficient)들을 정확하게 측정해서 적용해야 하겠지만 우리가 만드는 코드는 꼭 정확해야 할 필요는 없습니다.
왜냐하면, 프로그래머인 내가 창조하는 소프트웨어 속 세상에서 유체저항을 내가 원하는 대로 만든다고 해도 누가 뭐라 할 사람은 없습니다.
이해를 돕기 위해 위의 유체저항 공식을 재 구성해 보겠습니다.
위 수식에서 저항계수(C), 물체 앞면적(A), 액체 밀도(p)는 프로그래머가 편리하게 상수를 정해 쓰면 됩니다.
예) 저항계수(C) : 0.1, 물체 앞면적(A) : 무시, 액체 밀도(p) : 1 ...
이렇게 생각하면 위 유체저항 수식은 아래와 같이 간략히 만들어 볼 수 있습니다.
- 저항력(힘, 벡터) = 속도(벡터) * 속도(벡터) * 저항계수
이렇게 얻어진 저항력은 뉴턴의 F = MA (A = F/M) 라는 수식에 따라 질량으로 나눈
후, 얻어진 힘을 가속도 벡터에 더하기(벡터합연산)만 하면 됩니다.
가속도는 속도에 영향을 미치므로,
-
속도 = 속도 + 가속도 (점점 빨라지게 되므로 속도 제한 필요)
- 위치 = 위치 + 속도
따라서 물체는 이동하게 됩니다.
질량이 낮은 물체는 유체의 표면장력(저항력)을 뚫지 못하고 튕겨난 후, 속도가
감소되면 유체로 진입하는 것을 볼 수 있습니다.
2. 중력 모방
왜 모방이냐면, 내가 만든 프로그램에서는 개발자가 신입니다.
아래 중력수식은 적용하되, 중력 상수는 적절한 계수(Coefficient)를 아무거나 입력해도 상관없다는 의미입니다.
참고로 중력 상수는 아래와 같습니다. (저는 그냥 1로)
| = (6.673 84 ± 0.000 0080) ×10−11 N m2 kg−2 |
중력은 질량이 있는 물체가 서로를 당기는 힘입니다.
아래는 중력을 구하는 공식입니다.
|
| [중력 공식] |
F는 중력이며, G는 중력상수, r은 두 물체 사이의 거리, m1, m2는 각 물체의 질량입니다.
예제에서 원에 그려진 숫자가 각 물체의 질량을 의미합니다.
r, 즉 두 벡터간 거리(크기)를 구하기 위해 먼저, 해당 벡터간 차연산을 통해 방향 벡터(Fd)를 구합니다. 이제 이 벡터의 크기는 피타고라스 정리를 이용해 구할 수 있습니다.
만약 아래 그램에서 A(2,2) B(6,5) 두 벡터가 있다면 직각 삼각형을 만들어 밑변과 높이를 구할 수 있습니다.
예) 벡터 A->B 간 거리(r)는 빗변제곱 = (4*4) + (3*3) 이므로, 5 입니다.
|
| [두 벡터간 거리 구하기] |
이제 중력은 위 수식에 따라 아래와 같이 구합니다.
만약 중력상수:1 , 물체1 질량:100, 물체2 질량:10, 거리:5 로 가정하면,
- 중력 크기(스칼라) = (1 * 100 * 10) / (5*5) = 40 입니다.
아직 끝이 아닙니다.
위에서 구한 방향 벡터(Fd)를 정규화(normalize) 시킨 후, 중력 크기(스칼라)와 곱하기 연산(벡터와 스칼라)를 수행하면 최종적으로 중력의 힘벡터를 구할 수 있습니다.
이후 과정은 저항력을 적용하는 것과 동일합니다.
이렇게 얻어진 중력(벡터)은 뉴턴의 F = MA (A = F/M) 라는 수식에 따라 질량으로 나눈 후, 얻어진 힘을 가속도 벡터에 더하기(벡터합연산)만 하면 됩니다.
위 예제에서 각 물체는 서로가 서로를 당기는 중력은 구현하지 않았으며, 중심 물체와 자신을 대상으로만 중력을 구해 적용하였습니다.
사실 이 부분(질량이 있는 물체끼리의 중력 구현)도 별로 어렵지는 않습니다.
반복문을 이용해 나 자신을 제외한 모든 물체 사이의 거리를 구하고 해당 수식을 적용해 중력을 구한 후, 가속도벡터에 중력을 더해주면 끝입니다.
어떤가요? 현실의 물리법칙이 벡터 연산을 통해 동작하는 것을 확인 할 수 있습니다.
직접 코드로 만들어 보면 훨씬 더 재미있을 것입니다.
중력의 방향을 조절해 위로 물체가 떠다니게도 만들 수 있고, 유체의 저항계수를 바꿔 진흙탕을 구현하는것도 어렵지 않겠지요.
마지막으로, 이 말씀을 드리고 싶습니다.
개발 속도가 중요한 세상에서 게임엔진 or 물리엔진을 이용해 쉽고 편리하게 물리 현상을 바로 구현할 수도 있습니다.
다양한 물리 연산 라이브러리가 존재하며, 심지어 무료인 경우도 많습니다.
또한, 우리가 공부해가며 만드는 것보다 버그도 적을 것입니다.
"그런데 왜 우리가 벡터를 공부해가며 물리 알고리즘을 직접 구현해야 하는 하는거지?" 라는 의문이 머리속을 맴돌 것입니다.
주된 이유는 벡터를 이용한 이동, 속도, 회전, 힘, 등 기본적인 물리, 수학에 대한 관심과, 이해가 없다면 물리엔진을 제대로 사용할 수 없습니다.
또한 스스로 코드를 만들어가는 과정에서 효율적인 코드를 고민하게 되고, 객체 지향적인 코드를 만들려는 시도와 오류처리 과정에서 배우는 지식은 좋은 프로그래머가 되는 밑거름이 될 것입니다.
내가 직접 Low-Level 에서 힘이 무엇인지, 벡터를 어떻게 활용하는지 등을 코드로 만들어 보고, 이해하는 시도는 스스로를 진정한 Code Guru의 길로 안내할 것임을 확신합니다.
감사합니다.
[참고 자료]
Nature of Code , 다니엘 쉐프만
예제에 사용된 전체 소스코드
-
유체 저항은 아래 1번 코드 + 3번 코드(mover.py, vector.py), 총 3개의 파일로 구성.
-
중력 모방은 아래 2번 코드 + 3번 코드(mover.py, vector.py), 총 3개의 파일로 구성.
각 프로젝트별 3개의 파일을 같은 경로에 두고, 시작 파일은 1 or 2번 코드로
잡고 실행하면 됩니다.
1. 유체 저항 소스코드
import sys
from PyQt5.QtWidgets import QApplication, QWidget
from PyQt5.QtCore import Qt, QRectF, QPointF
from PyQt5.QtGui import QPainter, QBrush, QColor
from Physics.vector import vector
from Physics.mover import Mover
from threading import Thread
import time
from random import randint
QApplication.setAttribute(Qt.AA_EnableHighDpiScaling, True)
class CDragForce:
def __init__(self, rect, c=0.1):
# 저항력 영역(QRectF)
self.rect = rect
# 저항계수(Coefficient)
self.c = c
def isInRect(self, pt):
return self.rect.contains(pt)
def calcDragForce(self, mover):
#저항력 = 속력제곱 * 저항계수 * 속도반대방향
speed = mover.velocity.magnitude()
dragMag = self.c * speed * speed
#속도 반대
dragforce = mover.velocity
dragforce *= vector(-1,-1)
#벡터 정규화
dragforce.normalize()
dragforce *= vector(dragMag, dragMag)
return dragforce
class CWidget(QWidget):
def __init__(self):
super().__init__()
#이동 물체 다수
self.mover = []
gap = 50
for i in range(12):
mover = Mover(mass=(i+1)*2, limit = 10)
mover.location.x = (i*gap)+gap
mover.location.y = 0
self.mover.append(mover)
#저항영역
self.drag = CDragForce(QRectF(0, self.height()/2, self.width(), self.height()/2), 1)
self.thread = Thread(target=self.threadFunc)
self.bThread = False
self.initUI()
def initUI(self):
self.setWindowTitle('force')
self.bThread = True
self.thread.start()
self.show()
def mousePressEvent(self, e):
gap = 50
i=0
for m in self.mover:
m.location.x = (i*gap)+gap
m.location.y = 0
m.velocity = vector(0,0)
m.acceleration = vector(0,0)
i+=1
def paintEvent(self, e):
qp = QPainter()
qp.begin(self)
for i in self.mover:
d = i.mass
r = d/2
rect = QRectF(i.location.x-r, i.location.y-r, d, d)
qp.setBrush(i.color)
qp.drawEllipse(rect)
qp.setBrush(QColor(0,162,232,120))
qp.drawRect(self.drag.rect)
qp.drawText(self.rect(), Qt.AlignTop|Qt.AlignLeft, '바람 x방향: +0.1\n중력 y방향: +0.1\n마찰력(속도반대방향):마찰계수 0.1, 수직력 1\n저항력: 저항계수 1')
qp.end()
def threadFunc(self):
while self.bThread:
for i in self.mover:
d = i.mass
r = d/2
# 중력
gravity = vector(0.0, 0.1*i.mass)
i.applyForce(gravity)
# 마찰력 = -1 * 마찰계수(μ) * 수직력(N) * 단위벡터
u = 0.1
N = 1.0
# 마찰계수 * 수직력 = 마찰력 크기
frictionMag = u*N
# 마찰력은 속도의 반대방향으로 생성
friction = i.velocity
friction *= vector(-1, -1)
# 마찰력 방향구하기
friction.normalize()
# 마찰력의 크기 * 단위벡터 = 마찰력
friction *= vector(frictionMag, frictionMag)
i.applyForce(friction)
#저항력
if self.drag.isInRect(QPointF(i.location.x, i.location.y)):
dragforce = self.drag.calcDragForce(i)
i.applyForce(dragforce)
i.update()
if i.location.x+r > self.width():
i.location.x = self.width()-r
i.velocity.x *= -1
elif i.location.x-r < 0:
i.velocity.x *= -1
i.location.x = 0+r
if i.location.y+r > self.height():
i.location.y = self.height()-r
i.velocity.y *= -1
elif i.location.y-r < 0:
i.velocity.y *= -1
i.location.y = 0+r
self.update()
time.sleep(0.01)
def closeEvent(self, e):
self.bThread = False
if __name__ == '__main__':
app = QApplication(sys.argv)
w = CWidget()
sys.exit(app.exec_())
2. 중력 모방 소스코드
import sys
from PyQt5.QtWidgets import QApplication, QWidget
from PyQt5.QtCore import Qt, QRectF
from PyQt5.QtGui import QPainter, QBrush, QColor
from Physics.vector import vector
from Physics.mover import Mover
from threading import Thread
from random import randint
import time
QApplication.setAttribute(Qt.AA_EnableHighDpiScaling, True)
# 끌어 당기는 자(움직임 X)
class Attractor:
def __init__(self, x, y, mass):
#움직이지 않으므로, 질량과 위치만
self.mass = mass
self.location = vector(x, y)
#중력 상수
self.G = 1.0
# 최대 최소 거리 제한
self.min = self.mass/5
self.max = self.mass*3
def attract(self, mover):
# 두 물체의 방향 계산
force = self.location-mover.location
# 두 벡터간 거리
distance = force.magnitude()
# 최대 최소거리 제한(0나누기 문제 등)
distance = self.constrain(distance, self.min, self.max)
force.normalize()
# 중력공식
F = (self.G * self.mass) / (distance*distance)
force *= vector(F, F)
return force
def constrain(self, val, min, max):
if val < min:
return min
elif val > max:
return max
else:
return val
class CWidget(QWidget):
def __init__(self):
super().__init__()
# 끌어 당기는 자
self.attractor = Attractor(self.width()/2, self.height()/2, 100)
#이동 물체
self.mover = []
for i in range(10):
mover = Mover(randint(0, self.width()) , randint(0,self.height()), randint(10,30))
self.mover.append(mover)
self.thread = Thread(target=self.threadFunc)
self.bThread = False
self.initUI()
def initUI(self):
self.setWindowTitle('force')
self.bThread = True
self.thread.start()
self.show()
def paintEvent(self, e):
qp = QPainter()
qp.begin(self)
# attractor의 지름을 질량으로 설정
d = self.attractor.mass
r = d/2
# attractor의 위치를 항상 중심으로
self.attractor.location.x = self.rect().width()/2
self.attractor.location.y = self.rect().height()/2
# attractor 그리기
rect = QRectF(self.attractor.location.x-r, self.attractor.location.y-r, d, d)
qp.setBrush(QColor(0,0,0,128))
qp.drawEllipse(rect)
qp.drawText(rect, Qt.AlignCenter, 'Mass:{}'.format(self.attractor.mass))
# mover 그리기
for m in self.mover:
d = m.mass
r = d/2
rect = QRectF(m.location.x-r, m.location.y-r, d, d)
qp.setBrush(m.color)
qp.drawEllipse(rect)
qp.drawText(rect, Qt.AlignCenter, '{}'.format(m.mass))
# 중력 범위
qp.setPen(Qt.NoPen)
qp.setBrush(QColor(0,0,0,30))
x = self.attractor.location.x
y = self.attractor.location.y
range = self.attractor.max
rect = QRectF(x-range/2, y-range/2, range, range)
qp.drawEllipse(rect)
qp.end()
def threadFunc(self):
while self.bThread:
for m in self.mover:
v = self.attractor.attract(m)
m.applyForce(v)
m.velocity += m.acceleration
m.velocity.setLimit(1)
m.location += m.velocity
d = m.mass
r = d/2
if m.location.x+r > self.width():
m.location.x = self.width()-r
m.velocity.x *= -1
elif m.location.x-r < 0:
m.velocity.x *= -1
m.location.x = 0+r
if m.location.y+r > self.height():
m.location.y = self.height()-r
m.velocity.y *= -1
elif m.location.y-r < 0:
m.velocity.y *= -1
m.location.y = 0+r
# 가속도를 0 설정
# 뉴턴의 제1법칙(관성)에 따라 속도는 유지
m.acceleration*=vector(0,0)
self.update()
time.sleep(0.01)
def closeEvent(self, e):
self.bThread = False
if __name__ == '__main__':
app = QApplication(sys.argv)
w = CWidget()
sys.exit(app.exec_())
3. 공통 소스코드 (vector.py, mover.py)
1, 2번 코드에 공통적으로 필요한 모듈
3-1 vector.py
import math
class vector:
def __init__(self, x=0.0, y=0.0, limit=9999):
self.x = x
self.y = y
self.limit = limit
def __add__(self,other):
return vector(self.x+other.x, self.y+other.y)
def __sub__(self, other):
return vector(self.x-other.x, self.y-other.y)
def __mul__(self, other):
return vector(self.x*other.x, self.y*other.y)
def __truediv__(self, other):
return vector(self.x/other.x, self.y/other.y)
#벡터 크기
def magnitude(self):
# 피타고라스 정의 (빗변 구하기)
return math.sqrt(self.x*self.x + self.y*self.y)
# 벡터 내적
def dotVector(self, other):
return (self.x*other.x) + (self.y*other.y)
# 벡터 외적 (3D에서 사용)
def corssVector(self, other1, other2):
pass
def angleVector(self):
# -y 윈도우 좌표계는 y가 꺼꾸로
theta = math.atan2(0-self.y, self.x)
deg = theta * 180.0 / math.pi;
if deg < 0:
deg += 360;
return deg
def angleBetweenVector(self, v2):
v = vector(self.x, self.y, self.limit)
v.normalize()
v2.normalize()
theta = v.dotVector(v2)
theta = math.acos(theta)
deg = theta * 180.0 / math.pi
return deg
# 벡터 정규화(방향유지, 크기 1로)
def normalize(self):
# 피타고라스 정의 (빗변 구하기)
mag = self.magnitude()
if mag > 0:
self.x /= mag
self.y /= mag
def setLimit(self, limit):
self.limit = limit
#copysign (x, y) y의 부호만 취해 x에 적용
if abs(self.x) > self.limit:
self.x = math.copysign(limit, self.x);
if abs(self.y) > self.limit:
self.y = math.copysign(limit, self.y)
3.2 mover.py
from Physics.vector import vector
from PyQt5.QtGui import QColor
from random import randint
class Mover:
def __init__(self, x=0.0, y=0.0, mass = 1.0, limit = 1.0):
self.location = vector(x,y)
self.velocity = vector()
self.acceleration = vector()
# 질량
self.mass = mass
self.G = 1.0
self.limit = limit
self.color = QColor(randint(0,255), randint(0,255), randint(0,255), 128)
def applyForce(self, force):
# 뉴턴 운동 2법칙 (F=MA or A = F/M)
force/=vector(self.mass, self.mass)
self.acceleration += force
def update(self):
self.velocity += self.acceleration
self.velocity.setLimit(self.limit)
self.location += self.velocity
# 가속도를 0 설정
# 뉴턴의 제1법칙(관성)에 따라 속도는 유지
self.acceleration*=vector(0,0)
댓글
댓글 쓰기