깜지로 쓰려는 블로그

2차원 공간 발사체의 최적 발사 각도 구하기

(1) 발사체 궤적 플로팅

발사체의 최대거리 함수 r(theta)가 있을시, 거리를 최대로 하는 theta를 찾는 문제를 다룬다고 하자

일단 발사체 궤적과 플로팅 함수는 다음과 같다

vx, vz는 속도를 xy축으로 분해한 성분으로

궤적 루프는 z가 0이 될때 즉 발사체가 땅으로 떨어지면 종료한다.

import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from math import sin, cos, pi

def trajectory(theta, speed=20, height=0, dt=0.01, g=-9.81):
    vx = speed * cos(pi * theta / 180)
    vz = speed * sin(pi * theta / 180)
    t, x, z = 0, 0, height
    ts, xs, zs = [t], [x], [z]
    while z >= 0:
        t += dt
        vz += g * dt
        x += vx * dt
        z += vz * dt
        ts.append(t)
        xs.append(x)
        zs.append(z)
    return ts, xs, zs

def plot_trajectories(*trajs, show_seconds=False):
    for traj in trajs:
        xs, zs = traj[1], traj[2]
        plt.plot(xs, zs)
        if show_seconds:
            second_indices = []
            second = 0
            for i, t in enumerate(traj[0]):
                if t >= second:
                    second_indices.append(i)
                    second +=1
            plt.scatter([xs[i] for i in second_indices], [zs[i] for i in second_indices])
    xl = plt.xlim()
    plt.plot(plt.xlim(), [0, 0], c='k')
    plt.xlim(*xl)

    width = 7
    coords_height = (plt.ylim()[1] - plt.ylim()[0])
    coords_width = (plt.xlim()[1] - plt.xlim()[0])
    plt.gcf().set_size_inches(width, width * coords_height / coords_width)

일단 발사각 45, 60인 경우 발사 궤적은 이런식인데

45도일때 가장 멀리 나간다.

plot_trajectories(
    trajectory(45),
    trajectory(60)
)

(2) 발사체 궤적 성질보기

최대 발사각 구하려면 지금 처럼 모든 각도에 대한 궤적을 플로팅해서 가장 멀리가는걸 찾거나

r(theta)를 최대로 하는 theta를 최적화 기법을 찾아내야한다.

일단 궤적 성질이라하면 도착거리, 최대 높이, 체공 시간 정도가 있는데

궤적 함수로부터 다음과같이 뽑아낼 수 있음

def landong_position(traj):
    return traj[1][-1]

def hang_time(traj):
    return traj[0][-1]

def max_height(traj):
    return max(traj[2])

traj_45 = trajectory(theta=45)

print(landong_position(traj_45))
print(hang_time(traj_45))
print(max_height(traj_45))

(3) 여러 각도에서 도착 거리 보기

위에서 도착거리, 체공시간, 최대 높이에 대한 함수를 구현했다.

r(theta)는 발사체의 도착 거리이므로

여러 각도 값을 주어 산점도로 플로팅해보자

아까 구현한 궤적 함수의 경우 시간에 따른 도착 거리지만

지금 플로팅한것은 발사각 세타에 대한 높이를 나타낸다.

대충 봐도 45도쯤에서 발사 거리가 최대로 나온다.

angles = range(0, 90, 5)
landong_positions = [landong_position(trajectory(theta)) for theta in angles]
plt.scatter(angles, landong_positions)

(4) 여러 발사 각도에서 채공시간 보기

1초간 발사체 위치는 아까 플로팅 궤적 함수에서 옵션이 있어 true를 주면 이런식으로 나온다.

아깐 theta별 거리에 대한 함수를 만들어서 플로팅했다면

이번에는 theta별 채공 시간에 대해 만들어서 플로팅해보자

45도 즘에서 최대 거리가 됫던것과는 다르게

90도쯤에서 가장 채공시간이 높게 나온다.

test_angles = range(0, 181, 5)
hang_times = [hang_time(trajectory(theta)) for theta in test_angles]
plt.scatter(test_angles, hang_times)

 (5) 지표와 각도가 주어질때 산점도 플로팅 함수

아까 발사체 지표가 도착거리, 채공시간, 최대 높이 3가지가 있었는데

방금 구현한것처럼 지표를 주면 바로 산점도 띄우는 함수를 만들자

이런식으로 채공 시간을 플로팅했고,

최대 높이나 거리도 플로팅 가능

def plot_trajectories(metric, thetas, **settings):
    plt.scatter(thetas, [metric(trajectory(theta, **settings)) for theta in thetas])

plot_trajectories(hang_time, range(0, 181, 5))

최대 거리의 경우 발사 시작점 높이가 10인 가정시 다음과 같음

plot_trajectories(landong_position, range(0,90, 5), height=10)

(6) 함수로 최적 거리구하기

위에서 한내용은 대충 오일러 방법으로 구하고, 그래프 플로팅해서 눈대중으로 찾아보았다.

이번에는 함수로 구해봄

일단 45도때 궤적을 플롯해보자

지금까지 계속 보던 포물선 그래프가 나온다.

발사 거리는 z(t) = 0인 시점이 되는 t로 정해짐

z 방향 속도 z'(t) = 초기 속도 + 중력 속도(가속도 아님)로 정해지고 이를 적분하면

z 방향 거리 z(t)가 아래와 같이 나온다.

z(t)가 t에 대한 이차식이 되었으므로 근의 공식으로 해를 구할수 있음

t = -2vz/g를 

발사거리 r = vx * dt = |v| cos(theta) dt에 대입하여

각도별 거리

r(theta) = - 2|v|**2/g * sin(theta) * cos(theta)로 정리가 가능

다시 정리하자

1. z축 거리가 0이 되는 시간 t를 구함, z(t) = 0이되는 t를 구함

2. x축 거리 함수에 대입하여 theta별 거리 r(t)를 구함

이를 플로팅하면 아까 궤적 함수를 플로팅했을때랑 동일한 형태 결과가 나온다.

최대 거리를 구하려면? r(theta) = 발사 거리 함수가 주어졌을때

r'(theta) = 0, 극점이 되는 theta를 구하면 된다.

trj = trajectory(45)
ts, zs = trj[0], trj[2]
plt.plot(ts, zs)


def z(t):
    return 20 * sin(45*pi/180) * t + (-9.81/2)*t**2

def plot_function(f, xmin, xmax, **kwargs):
    ts = np.linspace(xmin, xmax, 1000)
    plt.plot(ts, [f(t) for t in ts], **kwargs)

plt.plot(ts,zs)
plot_function(z,0,2.9)



def r(theta):
    return (-2 * 20 * 20/-9.81) * sin(theta * pi / 180) * cos(theta * pi / 180)
plot_function(r, 0, 90)

3치원 궤적 다루기

(1) 3차원 궤적 플로팅

위에서 한 경우는 발사각 theta가 주어졌을때 xz축에서의 거리와 높이를 다루었음

이번엔 y축 phi를 추가하여 xyz 3d 공간에서의 궤적을 다룸

r(t) = 거리 했던거랑은 달리

이런 식으로 표현하면 되나 싶음

r(t, theta, phi) = (x, y, z)

def trajectory3d(theta, phi, speed= 20, height=0, dt=0.01, g=-9.81):
    vx = speed * cos(pi * theta / 180) * cos(pi * phi / 180)
    vy = speed * cos(pi * theta / 180) * sin(pi * phi / 180)
    vz = speed * sin(pi * theta / 180)
    t, x, y, z = 0, 0, 0, height
    ts, xs, ys, zs = [t], [x], [y], [z]
    while z >= 0:
        t += dt
        vz += g * dt
        x += vx * dt
        y += vy * dt
        z += vz * dt
        ts.append(t)
        xs.append(x)
        ys.append(y)
        zs.append(z)
    return ts, xs, ys, zs

def plot_trajectory3d(traj):
    fig = plt.gcf()
    fig.set_size_inches(7, 7)
    ax = fig.gca(projection='3d')
    ax.plot(traj[1], traj[2], traj[3])

plot_trajectory3d(trajectory3d(45, 45))

(2) 지형을 추가한 3차원 궤적 플로팅

위 3차원 궤적 플로팅에서는 지형이 존재하지 않음

3d 궤적함수와 플롯팅 함수에 지형 요소를 추가

지형 함수로는 릿지 함수 사용

(x^2 - 5 y ^2)/2500

궤적 계산시 z >= 지형(x, y)인 경우에만 루프하여 궤적 추적

3차원 플로팅시 지형 서피스도 띄움

traj(20, 0)

traj(20, 270)인 두 경우에 대한 3차원 궤적

def flat_ground(x, y):
    return 0

def ridge(x, y):
    return (x ** 2 - 5 * y **2) / 2500

def trajectory3d(theta, phi, speed=20, height=0, dt=0.01, g=-9.81, elevation=flat_ground):
    vx = speed * cos(pi * theta / 180) * cos(pi * phi / 180)
    vy = speed * cos(pi * theta / 180) * sin(pi * phi / 180)
    vz = speed * sin(pi * theta / 180)
    t, x, y, z = 0, 0, 0, height
    ts, xs, ys, zs = [t], [x], [y], [z]
    while z >= elevation(x, y):
        t += dt
        vz += g * dt
        x += vx * dt
        y += vy * dt
        z += vz * dt
        ts.append(t)
        xs.append(x)
        ys.append(y)
        zs.append(z)
    return ts, xs, ys, zs

def plot_trajectories_3d(*trajs, elevation=flat_ground, bounds=None, zbounds=None, shadows=False):
    fig, ax = plt.gcf(), plt.gca(projection='3d')
    fig.set_size_inches(7, 7)
    if not bounds:
        xmin = min([x for traj in trajs for x in traj[1]])
        xmax = max([x for traj in trajs for x in traj[1]])
        ymin = min([x for traj in trajs for x in traj[2]])
        ymax = max([x for traj in trajs for x in traj[2]])

        padding_x = 0.1 * (xmax-xmin)
        padding_y = 0.1 * (ymax-ymin)
        xmin -= padding_x
        xmax += padding_x
        ymin -= padding_y
        ymax += padding_y
    else:
        xmin, xmax, ymin, ymax = bounds
    plt.plot([xmin, xmax], [0, 0], [0, 0], c='k')
    plt.plot([0, 0], [ymin, ymax], [0, 0], c='k')
    g = np.vectorize(elevation)
    ground_x = np.linspace(xmin, xmax, 20)
    ground_y = np.linspace(ymin, ymax, 20)
    ground_x, ground_y = np.meshgrid(ground_x, ground_y)
    ground_z = g(ground_x, ground_y)
    ax.plot_surface(ground_x, ground_y, ground_z, cmap=cm.coolwarm, 
                    alpha=0.5,linewidth=0, antialiased=True)
    
    for traj in trajs:
        ax.plot(traj[1], traj[2], traj[3])

    if zbounds:
        ax.set_zlim(*zbounds)

plot_trajectories_3d(
    trajectory3d(20, 0, elevation=ridge),
    trajectory3d(20, 270, elevation=ridge),
    bounds=[0,40,-40,0],
    elevation=ridge
)

(3) 3차원 공간에서 발사체 궤적 거리 구하기

2d 차원에서 theta를 x축에 놓고 거리를 y축에 놓아 그래프로 표현했었음

이번에는 3d에서 다룸

고도 함수 h(x, y) = (x^2 - 5y^2) / 2500이므로

h(x, y) = Bx^2 - Cy^2로 놓으면

B = 1/2500

C= 5/2500이 된다.

발사체의 고도는 시간 t에 따라 변하므로 h(t)라 일단 하고

x, y는 vx * t, vy * t로 구할수 있다.

h(x, y)에 이를 대입하면

h(vx * t, vy * t) = B vx^2  t^2 - C vy^2 t^2 가 된다.

vx, vy는 나두고 t에 대한 이차방정식으로 정리후 근의 공식으로 t에 대한 해를 구하면

t = -vz / (g/s - B vx^2 + cvy^2)가 온다.

다시 정리해보자

지금까지 릿지 경사에서 채공시간을 구하는 과정을 다루었다.

1. 고도 함수 h(x, y)가 있다

2. 발사체 x, y는 vx * t, vy * t로 구한다.

3. 발사체 xy 좌표를 고도함수에 대입한다.

4. x, y에 대한 고도함수를 t에 대한 이차식으로 정리한다.

5. h(vx * t, vy * t) = 0가 되는 t를 근의 공식으로 구한다. <- 채공시간

채공 시간을 구하였으면 xy 평면에서의 이동 거리를 구하면 된다.

vxy * 채공시간을 하면 되는데

vxy = 피타고라스정리 vx, vy로 구하면 됨

2차원에서는 theta 한 변수만 하면됬지만

지금은 theta와 phi 두 변수를 다루므로 따로 플로팅 x

from math import sqrt
B = 0.001
C = 0.005
v = 20
g = -9.81

def velocity_component(v, theta, phi):
    vx = v * cos(theta * pi / 180) * cos(phi * pi / 180)
    vy = v * cos(theta * pi / 180) * sin(phi * pi / 180)
    vz = v * sin(theta * pi / 180)
    return vx, vy, vz

def landing_distance(theta, phi):
    vx, vy, vz = velocity_component(v, theta, phi)
    v_xy = sqrt(vx**2 + vy**2)
    a = (g/2) - B * vx**2 + C * vy **2
    b = vz
    landing_time = -b/a
    landing_distance = v_xy * landing_time
    return landing_distance

(4) 항력 요소를 발사체 궤적에 가미하기

실제 발사체는 공기저항을 받으므로 위 식그대로 날라가지는 않음

vx, vy, vz요소에 항력 영향을 가중치 0.1로하여 가미하면 다음과 같음

def trajectory3d(theta, phi, speed=20, height=0, dt=0.01, g=-9.81, elevation=flat_ground,drag=0):
    vx = speed * cos(pi * theta / 180) * cos(pi * phi / 180)
    vy = speed * cos(pi * theta / 180) * sin(pi * phi / 180)
    vz = speed * sin(pi * theta / 180)
    t, x, y, z = 0, 0, 0, height
    ts, xs, ys, zs = [t], [x], [y], [z]
    while z >= elevation(x, y):
        t += dt
        vx -= (drag * vx) * dt
        vy -= (drag * vy) * dt
        vz += (g - (drag * vz)) * dt
        x += vx * dt
        y += vy * dt
        z += vz * dt
        ts.append(t)
        xs.append(x)
        ys.append(y)
        zs.append(z)
    return ts, xs, ys, zs

plot_trajectories_3d(
    trajectory3d(20, -20, elevation=ridge),
    trajectory3d(20, -20, elevation=ridge, drag=0.1),
    bounds=[0,40,-40,0],
    elevation=ridge
)

(5) xy 방향각에 따른 발사체 거리 스칼라 필드를 히트맵으로 시각화

아까 landing_distance 함수 구하긴 했는데 theta, phi에 따라 거리를 어떻게 시각화하지는 못했는데

거리 스칼라장을 히트맵으로 보여줄수 있음

히트맵으로 플로팅한 결과는 다음과 같다

발사각(y or x 축 회전) theta는 0 ~ 90도 사이

발사 방향(z축 회전) phi는 0 ~ 360도

릿지 지형에서 theta는 40, phi는 100, 270즘에서 가장 멀리 날라간다.

def scalar_field_heatmap(f, xmin, xmax, ymin, ymax, xsteps=100, ysteps=100):
    fig = plt.figure()
    fig.set_size_inches(7, 7)
    fv = np.vectorize(f)
    X = np.linspace(xmin, xmax, xsteps)
    Y = np.linspace(ymin, ymax, ysteps)
    X, Y = np.meshgrid(X, Y)
    z = fv(X, Y)
    
    fig, ax = plt.subplots()
    c = ax.pcolormesh(X, Y, z, cmap='plasma')
    ax.axis([X.min(), X.max(), Y.min(), Y.max()])
    fig.colorbar(c, ax=ax)

scalar_field_heatmap(landing_distance, 0, 90, 0, 360)

(6) 스칼라장을 히트맵 대신 3d에서 플로팅하기

스칼라장을 히트맵대신 서피스로 플로팅한 결과

봉우리 2개가 나오는데

이 두 지점에선 기울기가 0이 됨

def plot_scalar_field(f, xmin, xmax, ymin, ymax, xsteps=100, ysteps=100, c=None,cmap=cm.coolwarm, alpha=1,antialiased=False):
    fig = plt.gcf()
    fig.set_size_inches(7, 7)
    ax = fig.gca(projection='3d')
    fv = np.vectorize(f)
    X = np.linspace(xmin, xmax, xsteps)
    Y = np.linspace(ymin, ymax, ysteps)
    X, Y = np.meshgrid(X, Y)
    Z = fv(X, Y)
    surf = ax.plot_surface(X, Y, Z, cmap=cmap, color=c, alpha=alpha,
                          linewidth =0, antialiased=antialiased)

plot_scalar_field(landing_distance, 0, 90, 0, 360)

경사 상승을 이용한 최적화

(1) 거리 히트맵과 경사 벡터장

경사 상승이라니 딥러닝 공부할때 생각나네

딥러닝에서는 비용 함수가 최소가 되는 방향으로 학습을 시켰는데

여기서는 비행체 거리가 최대가 되는 방향을 찾아서 경사 상승이라 하나보다

일단 경사와 관련된 함수부터 구현하자

도함수 공식으로 시컨트 슬로프 함수 만들고

시컨트 슬로프로 근사 미분 함수를

근사 미분 함수로 근사 그라디언트 xy 방향 편미분 반환하는 함수 만든다..

그리고 이 함수와 도착거리함수로 도착 거리 그라디언트 함수 구현

def secant_slope(f, xmin, xmax):
    return (f(xmax) - f(xmin)) / (xmax - xmin)

def approx_derivative(f, x, dx=1e-6):
    return secant_slope(f, x-dx, x+dx)

def approx_gradient(f, x0, y0, dx=1e-6):
    partial_x = approx_derivative(lambda x: f(x, y0), x0, dx=dx)
    partial_y = approx_derivative(lambda y: f(x0, y), y0, dx=dx)
    return (partial_x, partial_y)

def landing_distance_gradient(theta, phi):
    return approx_gradient(landing_distance, theta, phi)

벡터장 플로팅 함수로 히트맵에다가 화살표들을 올림 

from vectors import to_polar, to_cartesian

def plot_vector_field(f, xmin, xmax, ymin, ymax, xsteps=10, ysteps=10, color='k'):
    X, Y = np.meshgrid(np.linspace(xmin,xmax,xsteps), np.linspace(ymin, ymax, ysteps))
    U = np.vectorize(lambda x, y : f(x, y)[0])(X, Y)
    V = np.vectorize(lambda x, y : f(x, y)[1])(X, Y)
    plt.quiver(X, Y, U, V, color=color)
    fig = plt.gcf()

scalar_field_heatmap(landing_distance, 0, 90, 0, 360)
plot_vector_field(landing_distance_gradient, 0, 90, 0, 360, xsteps=10, ysteps=10, color='k')
plt.xlabel('theta')
plt.ylabel('phi')
plt.gcf().set_size_inches(9, 7)

(2) 경사 상승법 구현하기

위에서 근사치 그라디언트 함수를 작성하여 벡터장 플로팅까지 했다

이번에는 이 근사 그라디언트로 경사 방향을 계속 올라 최대 지점까지 가는 함수 구현

36, 83을 줄시

끝점 37.58, 89.9에 도달함

from vectors import length

def gradient_ascent(f, xstart, ystart, tolerance=1e-6):
    x = xstart
    y = ystart
    grad = approx_gradient(f, x, y)
    while length(grad) > tolerance:
        x += grad[0]
        y += grad[1]
        grad = approx_gradient(f, x, y)
    return x, y

gradient_ascent(landing_distance, 36, 83)

(3) 경사 상승 궤적 그리기

그라디언트 어센트 함수로 시작점을 주면 최적의 지점까지 찾는건 잘 됬는데

궤적을 드로잉해보면 다음과 같음

def gradient_ascent_points(f, xstart, ystart, tolerance=1e-6):
    x = xstart
    y = ystart
    xs, ys = [x], [y]
    grad = approx_gradient(f, x, y)
    while length(grad) > tolerance:
        x += grad[0]
        y += grad[1]
        grad = approx_gradient(f, x, y)
        xs.append(x)
        ys.append(y)
    return xs, ys

scalar_field_heatmap(landing_distance, 35, 40, 80, 100)

plt.plot(*gradient_ascent_points(landing_distance, 36, 83), c='k')
plt.scatter([36,37.58114751557887],[83,89.99992616039857],c='k',s=75)
plt.xlabel('theta')
plt.ylabel('phi')
plt.gcf().set_size_inches(9, 7)

드디어 장을 다뤄보게 됬다

벡터장 플로팅 하기

좌표 x,y를 입력받아 반환하는 다음 함수를 주면 벡터를 반환

f(x, y) = (-2y, x)

import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return (-2*y, x)

def plot_vector_field(f, xmin, xmax, ymin, ymax, xstep=1, ystep=1):
    X, Y = np.meshgrid(np.arange(xmin, xmax, xstep), np.arange(ymin, ymax, ystep))
    U = np.vectorize(lambda x, y : f(x, y)[0])(X, Y)
    V = np.vectorize(lambda x, y : f(x, y)[1])(X, Y)
    plt.quiver(X, Y, U, V, color='red')
    fig = plt.gcf()
    fig.set_size_inches(7,7)


plot_vector_field(f, -5, 5, -5, 5)

x, y = 1, 2 인경우

f(1, 2) = (-5, 1)인 벡터를 반환

meshgrid xy 둘다 -5 ~ 5 범위로 step은 1로 해서

이런 식으로 반환

-5, -5가 f에 들어갔다고 하자

f(-5, -5) = (10, -5)가 된다.

U, V 행렬 첫 값을 보면 생각한대로 나온다.

그리고 plt.quiver()는 2d 화살표 장을 플로팅 해주는 함수라고 한다.

위 사진에서 잘렸는데

x,y = -5, -5에서 벡터는 (10, -5)인데 사진 밖에 나가서 보이질 않는다

대신 x, y = -3, -3이라고 하면 벡터는 (6, -3)

-3, -3위치의 화살표를 보면 대충 느낌은 맞는거같다

이번엔 함수를 f(x,y) = (-x, -y)로 바꿔보면 이런식으로 중앙을 향하는 벡터장이 나온다.

벡터장의 중심을 -2, 4로 이동시키기

함수 -2-x, 4-y로 하면 벡터장의 중심이 -2, 4로 이동한 결과가 나온다.

스칼라장 플로팅하기

이번에는 스칼라장을 띄워봄

포텐셜 함수 u(x, y) = 0.5 * (x^2 + y^2)를 플로팅하면

from matplotlib.pyplot import cm
def plot_scalar_field(
            f, xmin, xmax, ymin, ymax, xstep=0.25, ystep=0.25, 
            c=None, cmap=cm.coolwarm, alpha=1, antialiased=False
        ):
    fig = plt.figure()
    fig.set_size_inches(7, 7)
    ax = fig.gca(projection='3d')
    fv = np.vectorize(f)
    X = np.arange(xmin, xmax, xstep)
    Y = np.arange(ymin, ymax, ystep)
    X, Y = np.meshgrid(X, Y)
    Z = fv(X, Y)
    surf = ax.plot_surface(X, Y, Z, cmap=cmap, color=c, alpha=alpha,
                             linewidth=0, antialiased=antialiased)
def u(x, y):
    return 0.5 * (x**2 + y**2)
plot_scalar_field(u, -5, 5, -5, 5)

히트맵으로 플로팅하면 이런식

def scalar_field_heatmap(f, xmin, xmax, ymin, ymax, xstep=0.1, ystep=0.1):
    fig = plt.figure()
    fig.set_size_inches(7,7)
    fv = np.vectorize(f)
    X = np.arange(xmin, xmax, xstep)
    Y = np.arange(ymin, ymax, ystep)
    X, Y = np.meshgrid(X, Y)
    z = fv(X, Y)
    fig, ax = plt.subplots()
    c = ax.pcolormesh(X, Y, z, cmap='plasma')
    ax.axis([X.min(), X.max(), Y.min(), Y.max()])
    fig.colorbar(c, ax=ax)
    plt.xlabel('x')
    plt.ylabel('y')
scalar_field_heatmap(u, -5, 5, -5, 5)

이번엔 스칼라장을 등고선으로 플로팅

def scalar_field_contour(f, xmin, xmax, ymin, ymax, levels=None):
    fv = np.vectorize(f)
    X = np.arange(xmin, xmax, 0.1)
    Y = np.arange(ymin, ymax, 0.1)
    X, Y = np.meshgrid(X, Y)
    Z = fv(X, Y)
    fig, ax = plt.subplots()
    CS = ax.contour(X, Y, Z, levels=levels)
    ax.clabel(CS, inline=1, fontsize=10, fmt='%1.1f')
    plt.xlabel('x')
    plt.ylabel('y')
    fig.set_size_inches(7, 7)

scalar_field_contour(u, -10, 10, -10, 10, levels=[10, 20, 30, 40, 50, 60])

지금까지 벡터장과 스칼라장을 플로팅 시켰고

이번에는 우주선 게임에 중력장을 추가

우주선 게임에 중력장 추가하기

진행할때마다 내용이 달라져서 처음부터 구현함

1. 폴리곤 모델 구현하기

이번 구현 게임은 화면 벽에 닿으면 반대로 튀어나갈수 있게 bounce 여부를 추가

맨앞에 bounce 있으며 bounce 켜짐 여부에 따라 반대편에 튀어나오게 할지, 벽에 부딪히면 튕겨지도록 move에 구현됨

그외 차이는 변수에 gravity와 mass가 추가되고

move 함수 내에 가속도를 나타내는 추진력 벡터와 중력장 파트가 추가됨

최종 가속도 = 추력 가속도 + 중력 가속도

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform
from linear_solver import do_segments_intersect
import numpy as np

boundce = True

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.gravity = 0
        self.mass = 1
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
        self.draw_center = False
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds, thrust_vector=(0, 0), gravity_sources=[]):
        ax, ay = thrust_vector
        gx, gy = gravitational_field(gravity_source, self.x, self.y)
        ax += gx
        ay += gy
        self.vx += ax * milliseconds / 1000
        self.vy += ay * milliseconds / 1000
        
        dx, dy = self.vx * milliseconds / 1000, self.vy * milliseconds / 1000
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        
        if boundce:
            if self.x < -10 or self.x > 10:
                self.vx = - self.vx
            if self.y < -10 or self.y > 10:
                self.vy = - self.vy
        else:
            if self.x < -10:
                self.x += 20
            if self.y < -10:
                self.y += 20
            if self.x > 10:
                self.x -= 20
            if self.y > 10:
                self.y -= 20
        self.rotation_angle += self.angular_velocity * milliseconds / 1000
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]
    
    def does_collide(self, other_poly):
        for other_segment in other_poly.segments():
            if self.does_intersect(other_segment):
                return True
        return False

    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False

2. 중력장 구현하기

폴리곤 모델로 블랙홀 클래스 구현

블랙홀은 폴리곤 모델 move 함수의 gravity sources로 들어감

블랙홀 1개만이 아닌 여러개 영향을 줄수 있음

gravitational_field에서는 (x, y) 벡터 반환하나

블랙홀들의 x축 중력, y축 중력들을 합하여 반환해줌

class BlackHole(PolygonModel):
    def __init__(self, gravity):
        vs = [vectors.to_cartesian((0.5, 2 * pi * i / 20)) for i in range(0, 20)]
        super().__init__(vs)
        self.gravity = gravity
    

def gravitational_field(sources, x, y):
    fields = [vectors.scale(-source.gravity, (x-source.x, y-source.y))
              for source in sources]
    return vectors.add(*fields)

black_hole = BlackHole(0.1)
black_hole.x, black_hole.y = 0, 0
black_holes = [
    black_hole
]

3. 메인함수 제외한 나머지 구현

먼저 우주선, 소행성 구현

변경사항은 우주선 시작점 빼곤 없음

class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0,), (-0.25, 0.25), (-0.25, -0.25)])
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x + dist * cos(self.rotation_angle), y + dist*sin(self.rotation_angle))

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i/sides))
              for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)


ship = Ship()
ship.x = 7
ship.y = 3
asteroid_count = 10
default_asteroids = [Asteroid() for _ in range(0, asteroid_count)]
for ast in default_asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)

나머지는 컬러, 화면크기, 드로잉 관련 함수

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE =  (0, 0, 255)
GREEN = (0, 255, 0)
RED =   (255, 0, 0)
LIGHT_GRAY =  (240, 240, 240)
DARK_GRAY = (128, 128, 128)
width, height = 400, 400

def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=BLACK, fill=False):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    if fill:
        pygame.draw.polygon(screen, color, pixel_points, 0)
    else:
        pygame.draw.lines(screen, color, True, pixel_points, 2)


def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.line(screen, color, to_pixels(*v1), to_pixels(*v2), 2)

def draw_grid(screen):
    for x in range(-9, 10):
        draw_segment(screen, (x, -10), (x, 10), color=LIGHT_GRAY)
    for y in range(-9, 10):
        draw_segment(screen, (-10, y), (10, y), color=LIGHT_GRAY)
    draw_segment(screen, (-10, 0), (10, 0), color=DARK_GRAY)
    draw_segment(screen, (0, -10), (0, 10), color=DARK_GRAY)

thrust = 3

메인함수도 크게 바뀐건 없지만

우주선, 소행성 무브에 추력 벡터와 중력 소스가 추가되었으며

추력 벡터를 계산하는 과정

블랙홀을 드로잉하는 과정등이 추가되었음

전후진 입력없으면 추력 벡터는 (0, 0)이지만 

입력있을시 우주선 방향에 크기 3인 벡터를 x, y축 방향으로 분해해서 추력 벡터로 무브에 전달

바운스도 잘되고 블랙홀을 중심으로 움직인다.

def main(asteroids=default_asteroids):
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    pygame.display.set_caption('asteroids')
    done = False
    clock = pygame.time.Clock()

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        milliseconds = clock.get_time()
        keys = pygame.key.get_pressed()

        if keys[pygame.K_LEFT]:
            ship.rotation_angle += milliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= milliseconds * (2 * pi / 1000)
        for ast in asteroids:
            ast.move(milliseconds, (0, 0), gravity_sources=black_holes)
        
        ship_thrust_vector = (0, 0)
        if keys[pygame.K_UP]:
            ship_thrust_vector = vectors.to_cartesian((thrust, ship.rotation_angle))
        elif keys[pygame.K_DOWN]:
            ship_thrust_vector = vectors.to_cartesian((-thrust, ship.rotation_angle))
        ship.move(milliseconds, ship_thrust_vector, black_holes)
        laser = ship.laser_segment()


        screen.fill(WHITE)
        draw_grid(screen)
        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)
        draw_poly(screen, ship)
        for bh in black_holes:
            draw_poly(screen, bh, fill=True)
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast, color=GREEN)
        pygame.display.flip()
    pygame.quit()

if __name__ == "__main__":
    main()

이번엔 블랙홀을 2개로 늘린경우

소행성들이 두 블랙홀의 영향을 받아 움직인다.

black_hole = BlackHole(0.1)
black_hole.x, black_hole.y = 0, 0


black_hole2 = BlackHole(0.2)
black_hole.x, black_hole.y = -4, -4

black_holes = [
    black_hole,
    black_hole2
]

전체 코드

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform
from linear_solver import do_segments_intersect
import numpy as np

boundce = True

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.gravity = 0
        self.mass = 1
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
        self.draw_center = False
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds, thrust_vector=(0, 0), gravity_sources=[]):
        ax, ay = thrust_vector
        gx, gy = gravitational_field(gravity_sources, self.x, self.y)
        ax += gx
        ay += gy
        self.vx += ax * milliseconds / 1000
        self.vy += ay * milliseconds / 1000
        
        dx, dy = self.vx * milliseconds / 1000, self.vy * milliseconds / 1000
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        
        if boundce:
            if self.x < -10 or self.x > 10:
                self.vx = - self.vx
            if self.y < -10 or self.y > 10:
                self.vy = - self.vy
        else:
            if self.x < -10:
                self.x += 20
            if self.y < -10:
                self.y += 20
            if self.x > 10:
                self.x -= 20
            if self.y > 10:
                self.y -= 20
        self.rotation_angle += self.angular_velocity * milliseconds / 1000
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]
    
    def does_collide(self, other_poly):
        for other_segment in other_poly.segments():
            if self.does_intersect(other_segment):
                return True
        return False

    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False





class BlackHole(PolygonModel):
    def __init__(self, gravity):
        vs = [vectors.to_cartesian((0.5, 2 * pi * i / 20)) for i in range(0, 20)]
        super().__init__(vs)
        self.gravity = gravity
    

def gravitational_field(sources, x, y):
    fields = [vectors.scale(-source.gravity, (x-source.x, y-source.y))
              for source in sources]
    return vectors.add(*fields)

black_hole = BlackHole(0.1)
black_hole.x, black_hole.y = 0, 0


black_hole2 = BlackHole(0.2)
black_hole.x, black_hole.y = -4, -4

black_holes = [
    black_hole,
    black_hole2
]




class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0,), (-0.25, 0.25), (-0.25, -0.25)])
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x + dist * cos(self.rotation_angle), y + dist*sin(self.rotation_angle))

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i/sides))
              for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)


ship = Ship()
ship.x = 7
ship.y = 3
asteroid_count = 10
default_asteroids = [Asteroid() for _ in range(0, asteroid_count)]
for ast in default_asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)


BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE =  (0, 0, 255)
GREEN = (0, 255, 0)
RED =   (255, 0, 0)
LIGHT_GRAY =  (240, 240, 240)
DARK_GRAY = (128, 128, 128)
width, height = 400, 400

def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=BLACK, fill=False):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    if fill:
        pygame.draw.polygon(screen, color, pixel_points, 0)
    else:
        pygame.draw.lines(screen, color, True, pixel_points, 2)

def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.line(screen, color, to_pixels(*v1), to_pixels(*v2), 2)

def draw_grid(screen):
    for x in range(-9, 10):
        draw_segment(screen, (x, -10), (x, 10), color=LIGHT_GRAY)
    for y in range(-9, 10):
        draw_segment(screen, (-10, y), (10, y), color=LIGHT_GRAY)
    draw_segment(screen, (-10, 0), (10, 0), color=DARK_GRAY)
    draw_segment(screen, (0, -10), (0, 10), color=DARK_GRAY)

thrust = 3


def main(asteroids=default_asteroids):
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    pygame.display.set_caption('asteroids')
    done = False
    clock = pygame.time.Clock()

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        milliseconds = clock.get_time()
        keys = pygame.key.get_pressed()

        if keys[pygame.K_LEFT]:
            ship.rotation_angle += milliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= milliseconds * (2 * pi / 1000)
        for ast in asteroids:
            ast.move(milliseconds, (0, 0), gravity_sources=black_holes)
        
        ship_thrust_vector = (0, 0)
        if keys[pygame.K_UP]:
            ship_thrust_vector = vectors.to_cartesian((thrust, ship.rotation_angle))
        elif keys[pygame.K_DOWN]:
            ship_thrust_vector = vectors.to_cartesian((-thrust, ship.rotation_angle))
        ship.move(milliseconds, ship_thrust_vector, black_holes)
        laser = ship.laser_segment()


        screen.fill(WHITE)
        draw_grid(screen)
        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)
        draw_poly(screen, ship)
        for bh in black_holes:
            draw_poly(screen, bh, fill=True)
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast, color=GREEN)
        pygame.display.flip()
    pygame.quit()

if __name__ == "__main__":
    main()

대수식을 심볼릭 프로그래밍으로 풀기 위해

자료 구조로 만들기 진행

제곱, 수, 심볼 변수 정의하기

설명할 내용없이 단순함

class Power():
    def __init__(self, base, exponent):
        self.base = base
        self.exponent = exponent
    
class Number():
    def __init__(self, number):
        self.number = number

class Variable():
    def __init__(self, symbol):
        self.symbol = symbol

Power(Variable("x"), Number(2))

곱 클래스 만들기

실제 연산하는 내용은 없지만 단순하게 이것도 구현

이건 3 x^2

class Product():
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2

Product(Number(3), Power(Variable("x"), Number(2)))

약간 복잡한 함수 정의하기

Function은 함수

Apply는 함수와 해당 함수 인자

다음 함수식은

(3x^2 + x) * sin(x)

class Sum():
    def __init__(self, *exps):
        self.exps = exps
    
class Function():
    def __init__(self, name):
        self.name = name

class Apply():
    def __init__(self, function, argument):
        self.function = function
        self.argument = argument

f_expression = Product(
    Sum(
        Product(
            Number(3),
            Power(
                Variable("x"),
                Number(2)
            )
        ),
        Variable("x")
    ),
    Apply(
        Function("sin"),
        Variable('x')
    )
)

cos(x^3 - 5)는 다음과 같이 표현

Apply(
    Function("cos"),
    Sum(
        Power(
            Variable("x"),
            Number("3")
        ),
        Number(-5)
    )
)

로그 연산 표현해보기

ln(y^z)는 원래 쉽게 표현할수 있으나

지금 구현한 함수로는 다음과 같이 표현가능

from math import log
def f(y, z):
    return log(y ** z)

Apply(
    Function("ln"),
    Power(
        Variable("y"),
        Variable("z")
    )
)

나눗셈 구현하기

(a + b) / 2같은 경우는 다음과 같이 표현

class Quotient():
    def __init__(self, numerator, denominator):
        self.numerator = numerator
        self.denominator = denominator

Quotient(
    Sum(
        Variable("a"),
        Variable("b")
    ),
    Number(2)
)

차 연산 구현해서 이차 방정식 판별식 구현하기

이차 방정식 y = ax^2 + bx + c의 판별식 D = b  ^ 2 - 4 ac를 다음과 같이 표현 가능

class Difference():
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2

Difference(
    Power(
        Variable("b"),
        Number(2)
    ),
    Product(
        Number(4),
        Product(
            Variable("a"),
            Variable("c")
        )
    )
)

차연산 표현하기

차연산도 간단하게 이런식으로

class Negative():
    def __init__(self, exp):
        self.exp = exp

Negative(
    Sum(
        Power(
            Variable("x"),
            Number(2)
        ),
        Variable("y")
    )
)

2차 방정식 근의 공식 표현하기

근의 공식은 길어지지만 이런식으로

A = Variable('a')
B = Variable('b')
C = Variable('c')
Sqrt = Function('sqrt')

Quotient(
    Sum(
        Negative(B),
        Apply(
            Sqrt,
            Difference(
                Power(
                    B, Number(2)
                ),
                Product(
                    Number(4),
                    Product(A, C)
                )
            )
        )
    ),
    Product(
        Number(2),
        A
    )
)

식에 있는 모든 변수 찾아내기

def f(x):

   return x^3

기존의 파이썬 함수는 결과만 반환하지 안에 뭐가 들어있는지 모르므로

앞서 동작은 안하지만 대수 방정식 구조를 만듬

일단 기존에 만든 것에 어떤 변수가 있는지 반환하는 함수 만듬

일단 변수만 찾아내는 함수니까 넣은 변수는 잘 찾는거같다.

def distinct_variables(exp):
    if isinstance(exp, Variable):
        return (set(exp.symbol))
    elif isinstance(exp, Number):
        return set()
    elif isinstance(exp, Sum):
        return set().union(*[distinct_variables(exp) for exp in exp.exps])
    elif isinstance(exp, Product):
        return distinct_variables(exp.exp1).union(distinct_variables(exp.exp2))
    elif isinstance(exp, Power):
        return distinct_variables(exp.base).union(distinct_variables(exp.exponent))
    elif isinstance(exp, Apply):
        return distinct_variables(exp.argument)
    else:
        raise TypeError("Not a valid expression")

print(distinct_variables(Variable("z")))
print(distinct_variables(Number(3)))
print(distinct_variables(Sum(Number(3), Number(4))))
print(distinct_variables(Sum(Number(3), Variable("x"))))
print(distinct_variables(Sum(Variable("x"), Variable("y"))))
print(distinct_variables(Power(Variable("x"), Number(3))))
print(distinct_variables(Product(Power(Variable("x"), Number(3)), Variable("y"))))

방정식 동작하게 하기

아까 구현한 변수, 곱같은 클래스에 evaluate 함수 추가

evaluate로 변수를 바인딩해서 계산 수행

from abc import ABC, abstractmethod

class Expression(ABC):
    @abstractmethod
    def evaluate(self, **bindings):
        pass

class Number(Expression):
    def __init__(self, number):
        self.number = number
    def evaluate(self, **bindings):
        return self.number

class Variable(Expression):
    def __init__(self, symbol):
        self.symbol = symbol
    def evaluate(self, **bindings):
        try:
            return bindings[self.symbol]
        except:
            raise KeyError("Variable '{}' is not bound.".format(self.symbol))

class Product(Expression):
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2
    def evaluate(self, **bindings):
        return self.exp1.evaluate(**bindings) * self.exp2.evaluate(**bindings)

Product(Variable('x'), Variable('y')).evaluate(x=2, y=5)

타 연산 구현하기

방금 변수와 곱을 구현했고

나머지 연산들을 구현해보자

sin, cos, ln 같은 초월함수는 따로 Apply에서 바인딩 되도록 하고

합, 제곱, 차, 나눗셈 연산은 기본 연산대로 evaluate 추가

import math
from math import sin, cos, log

_function_bindings = {
    "sin":math.sin,
    "cos":math.cos,
    "ln":math.log
}

class Apply(Expression):
    def __init__(self, function, argument):
        self.function = function
        self.argument = argument
    def evaluate(self, **bindings):
        return _function_bindings[self.function.name](self.argument.evaluate(**bindings))


class Sum(Expression):
    def __init__(self, *exps):
        self.exps = exps
    def evaluate(self, **bindings):
        return sum([exp.evaluate(**bindings) for exp in self.exps])

class Power(Expression):
    def __init__(self, base, exponent):
        self.base = base
        self.exponent = exponent
    def evaluate(self, **bindings):
        return self.base.evaluate(**bindings) ** self.exponent.evaluate(**bindings)

class Difference(Expression):
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2
    def evaluate(self, **bindings):
        return self.exp1.evaluate(**bindings) - self.exp2.evaluate(**bindings)

class Quotient(Expression):
    def __init__(self, numerator, denominator):
        self.numerator = numerator
        self.denominator = denominator
    def evaluate(self, **bindings):
        return self.numerate.evaluate(**bindings) / self.denominator.evaluate(**bindings)

(3 x ^2 + x) * sin(x)

이 식을 심볼릭 방법으로 한것과

그냥 파이썬 함수에 대입한  결과를 보면 동일하게 잘 연산되었음

f_expression = Product(
                Sum(
                    Product(
                        Number(3),
                        Power(
                            Variable("x"),
                            Number(2)
                        )
                    ),
                    Variable("x")
                ),
                Apply(
                    Function("sin"),
                    Variable("x")
                )
            )

print(f_expression.evaluate(x=5))

def f(x):
    return (3 * x ** 2 + x) * sin(x)
print(f(5))

표현식 전개하기

지금까지 구한 표현식을 그대로 출력하면 클래스랑 주소를 알려주는데

내용을 전개해주면 좋으니 우선 추상 클래스 expression을 다음과 같이 구현

class Expression(ABC):
    @abstractmethod
    def evaluate(self, **bindings):
        pass
    @abstractmethod
    def expand(self):
        pass
    @abstractmethod
    def display(self):
        pass
    def __repr__(self):
        return self.display()

지금까지 구현한

숫자, 변수, 연산들을 다시 구현해보자

먼저 가장 간단한 숫자와 변수

아까와는 달리 주소가 아니라 클래스명과 값이 잘 나온다

class Number(Expression):
    def __init__(self, number):
        self. number = number
    def evaluate(self, **bindings):
        return self.number
    def expand(self):
        return self
    def display(self):
        return "Number({})".format(self.number)

class Variable(Expression):
    def __init__(self, symbol):
        self.symbol = symbol
    def evaluate(self, **bindings):
        return bindings[self.symbol]
    def expand(self):
        return self
    def display(self):
        return "Variable(\"{}\")".format(self.symbol)
print(Number(4))
print(Variable("x"))

합, 차연산 구현하고 동작시켜보면 이런식

class Sum(Expression):
    def __init__(self, *exps):
        self.exps = exps
    def evaluate(self, **bindings):
        return sum([exp.evaluate(**bindings) for exp in self.exps])
    def expand(self):
        return (Sum(*[exp.expand() for exp in self.exps]))
    def display(self):
        return "Sum({})".format(",".join([e.display() for e in self.exps]))

class Difference(Expression):
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2
    def evaluate(self, **bindings):
        return self.exp1.evaluate(**bindings) - self.exp2.evaluate(**bindings)
    def expand(self):
        return self
    def display(self):
        return "Difference({},{})".format(self.exp1.display(), self.exp2.display())

print(Sum(Number(3), Number(5)))
print(Sum(Variable("x"), Number(3)))
print(Difference(Number(3), Number(5)))
print(Difference(Variable("x"), Number(3)))
print(Difference(Variable("x"), Number(3)).evaluate(x=4))
print(Sum(Variable("x"), Difference(Number(3), Number(5))))

tmp_exp = Sum(Variable("x"), Number(3))
print(tmp_exp)
tmp_exp.expand()

나머지 연산도 구현

class Product(Expression):
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2
    def evaluate(self, **bindings):
        return self.exp1.evaluate(**bindings) * self.exp2.evaluate(**bindings)
    def expand(self):
        expanded1 = self.exp1.expand()
        expanded2 = self.exp2.expand()
        if isinstance(expanded1, Sum):
            return Sum(*[Product(e, expanded2).expand() for e in expanded1.exps])
        elif isinstance(expanded2, Sum):
            return Sum(*[Product(expanded1, e) for e in expanded2.exps])
        else:
            return Product(expanded1, expanded2)
    def display(self):
        return "Product({},{})".format(self.exp1.display(), self.exp2.display())

class Quotient(Expression):
    def __init__(self, numerator, denominator):
        self.numerator = numerator
        self.denominator = denominator
    def evaluate(self, **bindings):
        return self.numerator.evaluate(**bindings) / self.denominator.evaluate(**bindings)
    def expand(self):
        return self
    def display(self):
        return "Quotient({},{})".format(self.numerator.display(),self.denominator.display())

class Negative(Expression):
    def __init__(self, exp):
        self.exp = exp
    def evaluate(self, **bindings):
        return - self.exp.evaluate(**bindings)
    def expand(self):
        return self
    def display(self):
        return "Negative({})".format(self.exp.display())
    
class Power(Expression):
    def __init__(self, base, exponent):
        self.base = base
        self.exponent = exponent
    def evaluate(self, **bindings):
        return self.base.evaluate(**bindings) ** self.exponent.evaluate(**bindings)
    def expand(self):
        return self
    def display(self):
        return "Power({},{})".format(self.base.display(),self.exponent.display())

class Function():
    def __init__(self, name, make_latex=None):
        self.name = name
        self.make_latex = make_latex
    def latex(self, arg_latex):
        if self.make_latex:
            return self.make_latex(arg_latex)
        else:
            return " \\operatorname{{ {} }} \\left( {} \\right)".format(self.name, arg_latex)

class Apply(Expression):
    def __init__(self, function, argument):
        self.function = function
        self.argument = argument
    def evaluate(self, **bindings):
        return _function_bindings[self.function.name](self.argument.evaluate(**bindings))
    def expand(self):
        return Apply(self.function, self.argument.expand())
    def display(self):
        return "Apply(Function(\"{}\"),{})".format(self.function.name, self.argument.display())

식 간단하게 만들고 그냥 출력했을때랑

전개했을때 차이는 이런식

(a + b) * (y + z)

= ay + az + by + bz

출력 결과는 길어서 좀 그렇긴한데 결과도 잘 나오긴 함

Y = Variable('y')
Z = Variable('z')
A = Variable('a')
B = Variable('b')
Product(Sum(A, B), Sum(Y, Z))


Product(Sum(A, B), Sum(Y, Z)).expand()




f_expression = Product(
                Sum(
                    Product(
                        Number(3),
                        Power(
                            Variable("x"),
                            Number(2)
                        )
                    ),
                    Variable("x")
                ),
                Apply(
                    Function("sin"),
                    Variable("x")
                )
            )
f_expression.expand()

미분, 적분 내용까지 가긴 너무 힘드니 대충 sympy로 마무리

Sympy

기호로 수학 푸는 파이썬 lib

지금까지 구현한거랑 거의 비슷

숫자 그대로 써도되고, 심볼 쓸수있고

아래 내용은 방정식에 대입, 미분, 적분 결과

from sympy import *
from sympy.core.core import *

Mul(Symbol('y'), Add(3, Symbol('x')))



y = Symbol('y')
x = Symbol('x')
print(y * (3+x))
print(y * (3+x).subs(x, 1)) # substitude
print((x**2).diff(x)) # derivative
print((3*x**2).integrate(x)) # integrate

이번에는 오일러 방법으로 지난번에 만든 게임에 이동을 추가함

사실 지난번에 만든 게임에는 우주선 이동을 안넣엇지 운석 이동을 시키면서 이미 오일러 방법을 사용했었다.

변화 거리 = 속도 * 소요시간(밀리세컨드) / 1000

기존 위치 = 이전 위치 + 변화 거리

그래도 게임을 개선한건 맞으니 그냥 쭉 정리하면

오일러 방법으로 무브는 하는데

벽을 넘어가면 반대편에 나오도록 수정

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
        self.draw_center = True
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]
    
    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds / 1000
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        if self.x < -10:
            self.x += 20
        if self.y < -10:
            self.y += 20
        if self.x > 10:
            self.x -= 20
        if self.y > 10:
            self.y -= 20
        self.rotation_angle += self.angular_velocity * milliseconds / 1000
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]

    def does_collide(self, other_poly):
        for other_segment in other_poly.segments():
            if self.does_intersect(other_segment):
                return True
        return False

    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False

우주선, 운석도 변경사항은 없음

class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25, -0.25)])
    
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x+dist*cos(self.rotation_angle), y + dist * sin(self.rotation_angle))

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5, 1.0), 2 * pi * i / sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)

폴리 그리기 함수에서 중심점 파트랑

그리드 그리기 함수가 추가

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE =  (0, 0, 255)
GREEN = (0, 255, 0)
RED =   (255, 0, 0)
LIGHT_GRAY =  (240, 240, 240)
DARK_GRAY = (128, 128, 128)
width, height = 400, 400

def to_pixel(x, y):
    return (width / 2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=BLACK):
    pixel_points = [to_pixel(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.lines(screen, color, True, pixel_points, 2)
    if polygon_model.draw_center:
        cx, cy = to_pixel(polygon_model.x, polygon_model.y)
        pygame.draw.circle(screen, BLACK, (int(cx), int(cy)), 4, 4)
    
def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.line(screen, color, to_pixel(*v1), to_pixel(*v2), 2)

def draw_grid(screen):
    for x in range(-9, 10):
        draw_segment(screen, (x, -10), (x, 10), color=LIGHT_GRAY)
    for y in range(-9, 10):
        draw_segment(screen, (-10, y), (10, y), color=LIGHT_GRAY)
    
    draw_segment(screen, (-10, 0), (10, 0), color=DARK_GRAY)
    draw_segment(screen, (0, -10), (0, 10), color=DARK_GRAY)

메인 함수 진입 전에

우주선이랑, 운석, 가속도 변수 초기화

메인함수가 지난번이랑 다른건

우주선 좌우전후, 그리드 그리기가 추가된게 다다

ship = Ship()
asteroid_count = 10
default_asteroids = [Asteroid() for _ in range(0, asteroid_count)]
for ast in default_asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)


acceleration = 3
def main(asteroids=default_asteroids):
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    done = False
    clock = pygame.time.Clock()

    while not done:
        clock.tick()

        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        milliseconds = clock.get_time()
        keys = pygame.key.get_pressed()

        for ast in asteroids:
            ast.move(milliseconds)
        
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += milliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= milliseconds * (2 * pi / 1000)
        if keys[pygame.K_UP]:
            ax = acceleration * cos(ship.rotation_angle)
            ay = acceleration * sin(ship.rotation_angle)
            ship.vx += ax * milliseconds/1000
            ship.vy += ay * milliseconds/1000
        if keys[pygame.K_DOWN]:
            ax = acceleration * cos(ship.rotation_angle)
            ay = acceleration * sin(ship.rotation_angle)
            ship.vx -= ax * milliseconds/1000
            ship.vy -= ay * milliseconds/1000
        ship.move(milliseconds)

        laser = ship.laser_segment()
        screen.fill(WHITE)
        draw_grid(screen)
        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)
        draw_poly(screen, ship)
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast, color=GREEN)
        pygame.display.flip()
    pygame.quit()

동작 결과는 이럼

우주선도 움직이니까 이젠 게임같긴하다

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform
from linear_solver import do_segments_intersect
import numpy as np
import sys

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
        self.draw_center = True
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]
    
    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds / 1000
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        if self.x < -10:
            self.x += 20
        if self.y < -10:
            self.y += 20
        if self.x > 10:
            self.x -= 20
        if self.y > 10:
            self.y -= 20
        self.rotation_angle += self.angular_velocity * milliseconds / 1000
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]

    def does_collide(self, other_poly):
        for other_segment in other_poly.segments():
            if self.does_intersect(other_segment):
                return True
        return False

    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False

class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25, -0.25)])
    
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x+dist*cos(self.rotation_angle), y + dist * sin(self.rotation_angle))

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5, 1.0), 2 * pi * i / sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)




BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE =  (0, 0, 255)
GREEN = (0, 255, 0)
RED =   (255, 0, 0)
LIGHT_GRAY =  (240, 240, 240)
DARK_GRAY = (128, 128, 128)
width, height = 400, 400

def to_pixel(x, y):
    return (width / 2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=BLACK):
    pixel_points = [to_pixel(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.lines(screen, color, True, pixel_points, 2)
    if polygon_model.draw_center:
        cx, cy = to_pixel(polygon_model.x, polygon_model.y)
        pygame.draw.circle(screen, BLACK, (int(cx), int(cy)), 4, 4)
    
def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.line(screen, color, to_pixel(*v1), to_pixel(*v2), 2)

def draw_grid(screen):
    for x in range(-9, 10):
        draw_segment(screen, (x, -10), (x, 10), color=LIGHT_GRAY)
    for y in range(-9, 10):
        draw_segment(screen, (-10, y), (10, y), color=LIGHT_GRAY)
    
    draw_segment(screen, (-10, 0), (10, 0), color=DARK_GRAY)
    draw_segment(screen, (0, -10), (0, 10), color=DARK_GRAY)


ship = Ship()
asteroid_count = 10
default_asteroids = [Asteroid() for _ in range(0, asteroid_count)]
for ast in default_asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)


acceleration = 3
def main(asteroids=default_asteroids):
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    done = False
    clock = pygame.time.Clock()

    while not done:
        clock.tick()

        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        milliseconds = clock.get_time()
        keys = pygame.key.get_pressed()

        for ast in asteroids:
            ast.move(milliseconds)
        
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += milliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= milliseconds * (2 * pi / 1000)
        if keys[pygame.K_UP]:
            ax = acceleration * cos(ship.rotation_angle)
            ay = acceleration * sin(ship.rotation_angle)
            ship.vx += ax * milliseconds/1000
            ship.vy += ay * milliseconds/1000
        if keys[pygame.K_DOWN]:
            ax = acceleration * cos(ship.rotation_angle)
            ay = acceleration * sin(ship.rotation_angle)
            ship.vx -= ax * milliseconds/1000
            ship.vy -= ay * milliseconds/1000
        ship.move(milliseconds)

        laser = ship.laser_segment()
        screen.fill(WHITE)
        draw_grid(screen)
        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)
        draw_poly(screen, ship)
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast, color=GREEN)
        pygame.display.flip()
    pygame.quit()

if __name__ == "__main__":
    main()

ㅡㅇㅂ재ㅔㄹㅈ배렂벫재

잠깐 딴짓했다가 다시 할 기분이 잘 안들었는데 겨우 시작

유량 다루기 시작 - 먼저 그래프 드로잉 함수부터

plot function은 이름 그대로 함수와 시간 구간 주면 그대로 드로잉 해주는 함수

여기서 사용한 함수는

유체량과 유체흐름률 두개가 있음

유체량은 (t-4)^3 / 64 + 3.3

유체흐름률은 유체량 함수의 도함수 3 * (t - 4)^2 / 64

0 ~ 10초 사이 유체량 함수를 드로잉한 결과는 다음과 같음

import matplotlib.pyplot as plt
import numpy as np

def plot_function(f, tmin, tmax, tlabel=None, xlabel=None, axes=False, **kwargs):
    ts = np.linspace(tmin, tmax, 1000)
    if tlabel:
        plt.xlabel(tlabel, fontsize=18)
    if xlabel:
        plt.ylabel(xlabel, fontsize=18)
    plt.plot(ts, [f(t) for t in ts], **kwargs)
    if axes:
        total_t = tmax-tmin
        plt.plot([tmin-total_t,tmax+total_t/10], [0, 0], c='k', linewidth=1)
        plt.xlim(tmin-total_t/10,tmax+total_t/10)
        xmin, xmax = plt.ylim()
        plt.plot([0, 0], [xmin, xmax], c='k', linewidth=1)
        plt.ylim(xmin, xmax)
def plot_volume(f,tmin,tmax,axes=False,**kwargs):
    plot_function(f, tmin, tmax, tlabel='time (hr)', xlabel='volume (bbl)', axes=axes, **kwargs)
def plot_flow_rate(f, tmin, tmax, axes=False, **kwargs):
    plot_function(f, tmin, tmax, tlabel='time (hr)', xlabel='flow rate(bbl/hr)', axes=axes, **kwargs)

def volume(t):
    return (t-4)**3 / 64 + 3.3
def flow_rate(t):
    return 3 * (t - 4) ** 2 / 64

plot_volume(volume, 0, 10)

특정 시간에서 유체의 량과 유체 량의 평균 변화를 구해보자

유체량에 대한 함수 volume에다가 시간 4, 9일때 량을 구하면 3.3, 5.2 가 나옴

평균 유체 흐름량 함수를 구현해보면 0.39가 나온다.

def average_flow_rate(v, t1, t2):
    return (v(t2) - v(t1))/(t2 - t1)

print(volume(4))
print(volume(9))
print(average_flow_rate(volume, 4, 9))

음의 변화율 다루기

아까 구현한 볼륨 함수는 점점 증가하는 형태였으나 이번엔

시간에 따라 내려가도록 구현

0 ~ 4 초사이 평균 흐름율은 -0.8

플로팅 결과는 저렇다

def decreasing_volume(t):
    if t < 5:
        return 10 - (t ** 2) / 5
    else:
        return 0.2 * (10 - t) ** 2

print(average_flow_rate(decreasing_volume, 0, 4))
plot_volume(decreasing_volume, 0, 10)

시컨트 라인 그리기

12시에 오도메트리가 77600 마일인데, 4시 30분에는 77900 마일이 됬다. 평균 속도 변화율은 어떨까?

으악 삼각함수도 잘 모르는데 시컨트라니

잠깐 찾아보니 삼각함수의 역수란다 이게 왜쓰는지 모르겠는데

코사인은 밑변/빗변이니 시컨트는 빗변/밑변 

흐름은 잘 모르겠지만 대충 시컨트 라인은 x1, x2사이 주어진 함수에 대한 직선 방정식을 구하는 함순거같음

def secant_line(f, x1, x2):
    def line(x):
        return f(x1) + (x - x1) * (f(x2) - f(x1)) / (x2 - x1)
    return line

def plot_secant(f, x1, x2, color='k'):
    line = secant_line(f, x1, x2)
    plot_function(line,x1,x2,c=color)
    plt.scatter([x1, x2], [f(x1), f(x2)], c=color)

plot_volume(volume, 0, 10)
plot_secant(volume, 4, 8)

구간별 평균 변화율

이번에는 1초 간격 평균 변화율을 구해보자

코드는 여전히 단순

def interval_flow_rate(v, t1, t2, dt):
    return [(t, average_flow_rate(v, t, t+dt)) for t in np.arange(t1, t2, dt)]

interval_flow_rate(volume, 0, 10, 1)

구간별 평균 유체 변화율 산점도 그리기

방금 작성한 구간별 평균 변화율로 산점도 드로잉하는 함수

def plot_interval_flow_rates(volume, t1, t2, dt):
    series = interval_flow_rate(volume, t1, t2, dt)
    times = [t for (t, _) in series]
    rates = [q for (_, q) in series]
    plt.scatter(times, rates)

plot_interval_flow_rates(volume, 0, 10, 1)

유체 흐름률과 구간별 평균 흐름률 플로팅

유체 흐름률은 유체량 함수의 도함수로 부터 구함

구간별 흐름률은 유체량에서 구간별 평균 변화율로 구함

모양은 비슷하나 후자가 살짝 뒤처짐

plot_flow_rate(flow_rate, 0, 10)
plot_interval_flow_rates(volume, 0, 10, 1)

하지만 구간을 1 -> 1/3으로 줄이면 원 함수에 더 가까워짐

plot_flow_rate(flow_rate, 0, 10)
plot_interval_flow_rates(volume, 0, 10, 1/3)

감소 하는 함수의 구간별 평균 흐름율 플로팅하면 이럼

plot_interval_flow_rates(decreasing_volume, 0, 10, 0.5)

직선의 경우 구간별 평균 흐름률

def linear_volume_function(t):
    return 5 * t + 3

plot_interval_flow_rates(linear_volume_function, 0, 10, 0.25)

유체량의 평균 변화율은 1에 가까워질수록 어떻게 될까?

아까 구현한 평균 유체 변화율로 1에 가까워지게 해보면 특정 값에 점점 가까워진다.

print(average_flow_rate(volume, 0.5, 1.5))
print(average_flow_rate(volume, 0.9, 1.1))
print(average_flow_rate(volume, 0.99, 1.01))
print(average_flow_rate(volume, 0.999, 1.001))
print(average_flow_rate(volume, 0.9999, 1.0001))

시컨트 라인과 탄젠트

계속 시컨트 시컨트하길래 잠깐 딴데 찾아보니 이렇게 설명한다

시컨트라인은 곡선 두 점 지나는 선

접선은 곡선의 한점을 지나는 선

순간 변화율 함수 구현하기

아깐 평균 변화율 구간을 1을 중심으로 좁힘으로서 특정 수에 수렴하는 것을 봤는데

이번에는 아예 순간변화율을 구하는 함수 구현됨

먼저 톨레랑스 용인도? 얼마까지 에러 봐줄래 하는 변수를 잡아줌

디폴트로 디지트 6했으니 10 ^ -6까지 오차는 용인한다.

h = 1로 놓고 위아래 1차이와 평균 변화율 구함

2 * digit 그러니까 12회 까지 반복수행하는대 오차가 10 ^-6 보다 크지않으면 반올림하고 반환

루프마다 h는 1/10배로 줄어듬

def instantaneous_flow_rate(v, t, digits=6):
    tolerance = 10 ** (-digits)
    h = 1
    approx = average_flow_rate(v, t-h, t+h)
    for i in range(0, 2 * digits):
        h = h / 10
        next_approx = average_flow_rate(v, t-h, t+h)
        if abs(next_approx - approx) < tolerance:
            return round(next_approx, digits)
        else:
            approx = next_approx
    raise Exception("derivative did not converge")

instantaneous_flow_rate(volume, 1)

구현한 순간 변화율 함수가 옳은지 플로팅 해보기

먼저 기존에 볼륨 도함수인 변화율 함수 플로팅

이번엔 순간 변화율 구하는 함수에 볼륨 함수를 집어넣은 뒤 플로팅

둘다 동일함

def get_flow_rate_function(v):
    def flow_rate_function(t):
        return instantaneous_flow_rate(volume, t)
    return flow_rate_function

plot_function(flow_rate, 0, 10)

plot_function(get_flow_rate_function(volume),0,10)

짧은 시간 동안 유체량 변화 정도 

1초 동안 유체량 변화 정도 계산시  0.1875나 실제로는 1.09가 나온다.

대신 구간을 1초 대신 0.01로 줄여보면 0.001875 실제로는 0.00186으로

구간을 좁히면 실제 값에 근사하게 나옴

def small_volume_change(q, t, dt):
    return q(t) * dt

print(small_volume_change(flow_rate, 2, 1))
print(volume(3) - volume(2))
print("#######")
print(small_volume_change(flow_rate, 2, 0.01))
print(volume(2.01) - volume(2))

볼륨 변화량 구하기

계속 유체라고 했는데 변수 명 그대로 볼륨이라고 해야겠다.

0 ~ 10 초간 0.1초 간격으로 볼륨 변화량을 합치면 4.3289가 나온다.

하지만 실제 볼륨10, 볼륨0으로 구해보면 4.375가 뜬다.

이번엔 구간을 0.0001로 줄이면 4.3749로 실제 값에 가까워짐

def volume_change(q, t1, t2, dt):
    return (sum(small_volume_change(q, t, dt) for t in np.arange(t1, t2, dt)))

print(volume_change(flow_rate, 0, 10, 0.1))
print(volume(10) - volume(0))
print(volume_change(flow_rate, 0, 10, 0.0001))

볼륨 플로팅하기

이번엔 실제 볼륨 함수와 근사화한 볼륨 함수를 플로팅함

볼륨 함수 근사화는 아까 구현한 볼륨 변화 사용

근사시 구간을 0.5로 해서 좀 안맞아보이긴 한데 거의 비슷하게 나오긴 함

근사 구간을 0.1로하면 아까 보다 원함수에 가까워짐

구간을 0.01로 더 줄이면 더가까워짐

def approximate_volume(q, v0, dt, T):
    return v0 + volume_change(q, 0, T, dt)

def approximate_volume_function(q, v0, dt):
    def volume_function(T):
        return approximate_volume(q, v0, dt, T)
    return volume_function

plot_volume(approximate_volume_function(flow_rate, 2.3, 0.5), 0, 10)
plot_volume(volume, 0, 10)



plot_volume(approximate_volume_function(flow_rate, 2.3, 0.1), 0, 10)
plot_volume(volume, 0, 10)




plot_volume(approximate_volume_function(flow_rate, 2.3, 0.01), 0, 10)
plot_volume(volume, 0, 10)

우주선, 소행성 게임 만들기

먼저 폴리곤 모델 정의

폴리곤 모델은 우주선, 행성의미

폴리곤 구성하는 점들의 모임과

x, y, 각도값

x속도, y속도, 각속도

로 구성

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0

폴리곤 모델로 배와 소행성 만들기

폴리곤 모델을 상속받는 우주선과 소행성 클래스 정의

소행성은 5 ~ 9개의 면을 가짐

소행성의 점들은 극좌표로 만든후 카티지안 좌표계로 변환

이때 극좌표 값 r은 0.5 ~ 1.0 균일분포로 뽑고 2 * pi * i /sides 로 theta를 뽑아냄

초기 속도 vx,vy는 0이나 각속도는 균일분포로 뽑아냄

import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform

class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25,-0.25)])
    
class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i /sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = 0
        self.vy = 0
        self.angular_velocity = uniform(-pi/2, pi/2)

게임 초기화

배와 소행성들을 생성

소행성의 xy는 -9 ~ 9 사이 임의의 정수로 뽑음

ship = Ship()
asteroid_count = 10
asteroids = [Asteroid() for _ in range(0, asteroid_count)]

for ast in asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)

게임 대충 켜보기

일단 지금 드로잉하는 게임은 width, height 400, 400이고

위에서 작성한 우주선과 소행성의 x, y가 아닌 점들은 0 ~ 1 사이로 정규화된 값을 갖는다.

정규화된 x, y값을 실제 게임화면에 맞추기 위한 to_pixel 함수와 비정규화된 xy좌표 값을 가진 pixel을 드로잉하는 함수를 작성

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE = (0, 0, 255)
GREEN = (0, 255, 0)
RED = (255, 0, 0)

width, height = 400, 400
def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=GREEN):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.aalines(screen, color, True, pixel_points, 10)

폴리곤 모델은 각도에 따라 포인트 점 위치가 달라지므로 transformed() 함수가 추가됨

init 만 있던 폴리곤 모델에 변환 함수 추가하고

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

지금 구현한 사항들을 pygame으로 돌리므로

최소한 화면에 띄을수 있는 코드를 작성해서 돌리면 

def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True

        screen.fill(WHITE)

        draw_poly(screen, ship)
        for ast in asteroids:
            draw_poly(screen, ast)
        pygame.display.flip()

    pygame.quit()
if __name__== "__main__":
    main()

우주선 1개랑

소행성 10개가 파이게임 스크린에 나온다.

지금까지 전체코드

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]


class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25,-0.25)])
    
class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i /sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = 0
        self.vy = 0
        self.angular_velocity = uniform(-pi/2, pi/2)

ship = Ship()
asteroid_count = 10
asteroids = [Asteroid() for _ in range(0, asteroid_count)]

for ast in asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE = (0, 0, 255)
GREEN = (0, 255, 0)
RED = (255, 0, 0)

width, height = 400, 400
def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=GREEN):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.aalines(screen, color, True, pixel_points, 10)


def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True

        screen.fill(WHITE)

        draw_poly(screen, ship)
        for ast in asteroids:
            draw_poly(screen, ast)
        pygame.display.flip()

    pygame.quit()
if __name__== "__main__":
    main()

소행성 이동시키기

화면에 띄운건 좋은데 아무것도 움직이지 못하니 좀 그렇다.

먼저 소행성부터 시간에따라 움직이도록 해보자

아직 polygonmodel에 move는 구현안했지만

일단 메인에 이런식으로 밀리세컨드마다 무브하도록 짚어놓고

def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        
        miliseconds = clock.get_time()
        for ast in asteroids:
            ast.move(miliseconds)

폴리곤 모델에 move를 구현한다.

해당 폴리곤의 밀리세컨드 시간에 따라 x, y와 회전각을 갱신시킨다.

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds/ 1000.0
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        self.rotation_angle += self.angular_velocity * milliseconds / 1000.0

근대 생각해보니 아까 소행성 vx, vy를 0으로 해놔서

위 코드에선 소행성이 회전만 한다

이동하도록 하기위해 vx, vy를 균일분포 -1, 1로 해놓자

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i /sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)

우주선이랑 소행성 구분이 안가는것 때문에

배는 draw poly에서 색을 빨강으로 설정

아무튼 이제 소행성이 회전도 하고 이동도 한다.

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds/ 1000.0
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        self.rotation_angle += self.angular_velocity * milliseconds / 1000.0


class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25,-0.25)])
    
class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i /sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)

ship = Ship()
asteroid_count = 10
asteroids = [Asteroid() for _ in range(0, asteroid_count)]

for ast in asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE = (0, 0, 255)
GREEN = (0, 255, 0)
RED = (255, 0, 0)

width, height = 400, 400
def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=GREEN):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.aalines(screen, color, True, pixel_points, 10)


def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        
        miliseconds = clock.get_time()
        for ast in asteroids:
            ast.move(miliseconds)

        screen.fill(WHITE)

        draw_poly(screen, ship, color="red")
        for ast in asteroids:
            draw_poly(screen, ast)
        pygame.display.flip()

    pygame.quit()
if __name__== "__main__":
    main()

우주선 회전하기

원래 우주선 회전 + 이동 둘다 구현이 되어있는줄 알았는데

깃헙 자료상 회전만 있다 .... 시간만 들이면 전후진 구현할수 있긴한데 너무 복잡하니 이번엔 회전 추가

복잡한건 없고, 위 코드의 메인함수에

키입력에 따라 배의 회전각만 변경시켜주는게 전부다

def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        
        miliseconds = clock.get_time()
        for ast in asteroids:
            ast.move(miliseconds)

        keys = pygame.key.get_pressed()
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += miliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= miliseconds * (2 * pi / 1000)

        screen.fill(WHITE)

        draw_poly(screen, ship, color="red")
        for ast in asteroids:
            draw_poly(screen, ast)
        pygame.display.flip()

    pygame.quit()

레이저 구현하기

소행성도 움직이고, 우주선도 회전하겠다 이젠 레이저 쏘기 구현하자.

레이저는 직선 형태니까 draw segment 함수로 구현, 2점이 주어질때 빨간색 드로잉하는 내용

main함수에는 ship.laser_segment()로 레이저 정보 가져오고,

space 눌렀을때 드로잉 한다.(?)

스페이스 안눌러도 레이저 날라가는게 남아있어야 하지 않나 싶긴한데 일단 그냥 진행

def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.aaline(screen, color, to_pixels(*v1), to_pixels(*v2), 10)


def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        
        miliseconds = clock.get_time()
        for ast in asteroids:
            ast.move(miliseconds)

        keys = pygame.key.get_pressed()
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += miliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= miliseconds * (2 * pi / 1000)

        laser = ship.laser_segment()

        screen.fill(WHITE)

        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)

        draw_poly(screen, ship, color="red")
        for ast in asteroids:
            draw_poly(screen, ast)
        pygame.display.flip()

    pygame.quit()

ship의 laser segment를 보니 왜 눌른 동안만 드로잉하게 된지 알겠다.

소행성과 우주선은 -10 ~ 10의 좌표 범위를 갖고있지만 to_pixels로 스크린에 맞게 400, 400으로 매칭되어있다.

레이저는 우주선 중앙부터 20 * root(2)로 대충 화면 너머까지 쏘개하도록 되어있음

레이저 날라가는 동작 없이 없이 발사하면 끝까지 나가도록 되어있기 때문

레이저 세그먼트는

레이저 길이를 20 * root 2

시작점은 우주선중앙

끝점은 이 레이저 길이 x cos, sin으로 우주선 방향 끝까지 넘어가도록 되어있음

class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25,-0.25)])
    
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x + dist * cos(self.rotation_angle), y + dist*sin(self.rotation_angle))

레이저를 쏘면 화면끝까지 날라간다.

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds/ 1000.0
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        self.rotation_angle += self.angular_velocity * milliseconds / 1000.0


class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25,-0.25)])
    
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x + dist * cos(self.rotation_angle), y + dist*sin(self.rotation_angle))
    
    
class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i /sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)

ship = Ship()
asteroid_count = 10
asteroids = [Asteroid() for _ in range(0, asteroid_count)]

for ast in asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE = (0, 0, 255)
GREEN = (0, 255, 0)
RED = (255, 0, 0)

width, height = 400, 400

def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=GREEN):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.aalines(screen, color, True, pixel_points, 10)

def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.aaline(screen, color, to_pixels(*v1), to_pixels(*v2), 10)


def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        
        miliseconds = clock.get_time()
        for ast in asteroids:
            ast.move(miliseconds)

        keys = pygame.key.get_pressed()
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += miliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= miliseconds * (2 * pi / 1000)

        laser = ship.laser_segment()

        screen.fill(WHITE)

        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)

        draw_poly(screen, ship, color="blue")
        for ast in asteroids:
            draw_poly(screen, ast)
        pygame.display.flip()

    pygame.quit()
if __name__== "__main__":
    main()

레이저와 소행성 충돌 판정하기

소행성 게임도 이제 거의 다되가지만 가장 중요한

레이저와 소행성 충돌 판정이 남아있다.

일단 메인 함수에 레이저랑 충돌하면 제거하고, 충돌안하면 그대로 드로잉 하도록 고치자

여기선 space가 눌린 상태에 ast.does_intersect(laser)가 true 시 소행성 제거하게된다.

def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        
        miliseconds = clock.get_time()
        for ast in asteroids:
            ast.move(miliseconds)

        keys = pygame.key.get_pressed()
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += miliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= miliseconds * (2 * pi / 1000)

        laser = ship.laser_segment()

        screen.fill(WHITE)

        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)

        draw_poly(screen, ship, color="blue")
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast)

        pygame.display.flip()

does_intersect는 소행성이 아니라 폴리곤 모델에 구현

직선간 겹침을 판단하기위해 폴리곤을 이루는 모든 직선은 세그먼트로 가져와 사용한다.

세그먼츠는 원래 폴리곤 모델의 포인트를 회전 변환후, 현재 인덱스와 다음 인덱스 두 점의 튜플 리스트

인자인 other_segment는 레이저이므로

레이저와 세그먼츠를 이루는 직선 세금먼트간에 교점있는지 여부를 do segments intersect에서 다룬다.

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds/ 1000.0
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        self.rotation_angle += self.angular_velocity * milliseconds / 1000.0
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]
    
    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False

연립일차방정식으로 레이저와 소행성 충돌 풀기

위에서 do_segments_intersect로 두 세그먼트가 주어질떄 교차되는지 판정하는 내용을 구현해야한다.

각 세그먼트는 두 점으로 이뤄져있으므로 각각 직선의 방정식을 구해서 해를 찾으면 된다.

두 직선의 교점, 연립일차방정식의 해를 구한다.

대충 소행성의 선분 세그먼트가 (2, 4)와 (4, 2)를 지나면 직선의 방정식으로 y = -x + 6

레이저 선분 세그먼트 직선의 방정식은 대충 y = x라 하자

그러면 일차식 표준형 ax + by = c로 바꾸면

좌측 하단과 같다.

이 두 일차식 표준형은 행렬로 다루면 

여기서 역행렬을 곱해주고 계산하면 두 직선의 교점 xy를 얻을수 있음

이경우엔 3, 3이 된다.

하지만 이런 경우 처럼 두 직선의 교점이 소행성 선분 밖에 존재할 수도 있음

이런 경우를 처리하기 위해

레이저 길이 d1, 소행성 선분 길이 d2로 놓고

다음 4가지 경우가 성립할떄, 두 직선이 제대로 만난다고 판단한다.

1. 소행성 선분 한점 u1, 교점 사이 길이가 d2 보다 작음

2. 소행성 선분 다른점 u2, 교점 사이 길이가 d2보다 작음

3. 레이저 한점 v1, 교점사이 길이가 d1보다 작음

4. 레이저 한점 v2, 교점사이 길이가 d1 보다 작음

아래 그림의 경우 교점이 (3,3)에 있어 레이저 길이 d1보다 작으므로 3, 4번 조건은 만족하지만

왼쪽 점을 u1이라 할경우, u1에서 (3, 3)거리는 d2보다 크므로 1번은 만족하지 않는다.

오른쪽 점을 u2라 하는 경우, u2에서 (3, 3)까지 거리는 d2보다 작아 만족한다.

이 4가지 경우를 곱하면 1번 조건을 불만족하므로 교차하지 않는 것으로 판단한다.

does_intersect 안에있는 do_segment_intersect를 구현하면 다음과 같다.

do_segment_intersect에선 세그먼트 s1, 세그먼트 s2를 받아

선분의 길이 d1, d2를 구하고

insersection에서 두 세그먼트의 교점 xy를 구한다.

아까 위에서 말한 4가지 경우에 따라 선분 길이 d1,d2와 xy교점과 세그먼트 점 길이로 충돌 여부를 판단한다.

intersection에서는 각 세그먼트를 구성하는 두점을 일차식 표준형으로 만들고

넘파이로 해를 구해서 반환한다.

나머지는 이 과정에 사용되는 전체 함수

def subtract(v1, v2):
    return tuple(v1-v2 for (v1,v2) in zip(v1, v2))

def length(v):
    return sqrt(sum([coord ** 2 for coord in v]))

def distance(v1, v2):
    return length(subtract(v1, v2))

def standard_form(v1, v2):
    x1, y1 = v1
    x2, y2 = v2
    a = y2 - y1
    b = x1 - x2
    c = x1 * y2 - y1 * x2
    return a, b, c

def intersection(u1, u2, v1, v2):
    a1, b1, c1 = standard_form(u1, u2)
    a2, b2, c2 = standard_form(v1, v2)
    m = np.array(((a1, b1), (a2, b2)))
    c = np.array((c1, c2))
    return np.linalg.solve(m, c)

def do_segments_intersect(s1, s2):
    u1, u2 = s1
    v1, v2 = s2
    d1, d2 = distance(*s1), distance(*s2)
    try:
        x, y = intersection(u1, u2, v1, v2)
        return (distance(u1, (x,y)) <= d1 and
                distance(u2, (x, y)) <= d1 and
                distance(v1, (x, y)) <= d2 and
                distance(v2, (x, y)) <= d2)
    except np.linalg.linalg.LinAlgError:
        return False

완성

이제 우주선에서 레이저를 쏘면 동작(평행이동, 회전)하는 소행성들의 충돌 여부를 판단해서

제거하는 간단한 게임 완성

import pygame
import vectors
import numpy as np
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform

def subtract(v1, v2):
    return tuple(v1-v2 for (v1,v2) in zip(v1, v2))

def length(v):
    return sqrt(sum([coord ** 2 for coord in v]))

def distance(v1, v2):
    return length(subtract(v1, v2))

def standard_form(v1, v2):
    x1, y1 = v1
    x2, y2 = v2
    a = y2 - y1
    b = x1 - x2
    c = x1 * y2 - y1 * x2
    return a, b, c

def intersection(u1, u2, v1, v2):
    a1, b1, c1 = standard_form(u1, u2)
    a2, b2, c2 = standard_form(v1, v2)
    m = np.array(((a1, b1), (a2, b2)))
    c = np.array((c1, c2))
    return np.linalg.solve(m, c)

def do_segments_intersect(s1, s2):
    u1, u2 = s1
    v1, v2 = s2
    d1, d2 = distance(*s1), distance(*s2)
    try:
        x, y = intersection(u1, u2, v1, v2)
        return (distance(u1, (x,y)) <= d1 and
                distance(u2, (x, y)) <= d1 and
                distance(v1, (x, y)) <= d2 and
                distance(v2, (x, y)) <= d2)
    except np.linalg.linalg.LinAlgError:
        return False


class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds/ 1000.0
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        self.rotation_angle += self.angular_velocity * milliseconds / 1000.0
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]
    
    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False


class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25,-0.25)])
    
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x + dist * cos(self.rotation_angle), y + dist*sin(self.rotation_angle))
    
    
class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i /sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)

ship = Ship()
asteroid_count = 10
asteroids = [Asteroid() for _ in range(0, asteroid_count)]

for ast in asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE = (0, 0, 255)
GREEN = (0, 255, 0)
RED = (255, 0, 0)

width, height = 400, 400

def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=GREEN):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.aalines(screen, color, True, pixel_points, 10)

def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.aaline(screen, color, to_pixels(*v1), to_pixels(*v2), 10)


def main():
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    clock = pygame.time.Clock()
    done = False

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        
        miliseconds = clock.get_time()
        for ast in asteroids:
            ast.move(miliseconds)

        keys = pygame.key.get_pressed()
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += miliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= miliseconds * (2 * pi / 1000)

        laser = ship.laser_segment()

        screen.fill(WHITE)

        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)

        draw_poly(screen, ship, color="blue")
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast)

        pygame.display.flip()

    pygame.quit()
if __name__== "__main__":
    main()

단순 덧셈 연산가능한 2차원 벡터 만들기

2d 벡터니 x, y만 가짐

덧셈 연산 구현

class Vec2():
    def __init__(self, x, y):
        self.x = x
        self.y = y
    
    def add(self, v2):
        return Vec2(self.x + v2.x, self.y + v2.y)

v = Vec2(3, 4)
w = v.add(Vec2(-2, 6))
print(w.x, w.y)

2차원 벡터에 스칼라곱, 동등 비교연산 추가

스칼라 곱, 비교 연산, 덧셈 연산 추가후

벡터와 스칼라간 곱  + 덧셈 연산 결과 repr이 없어주소가 반환됨

class Vec2():
    def __init__(self, x, y):
        self.x = x
        self.y = y
    
    def add(self, v2):
        return Vec2(self.x + v2.x, self.y + v2.y)

    def scale(self, scalar):
        return Vec2(scalar * self.x, scalar * self.y)
    
    def __eq__(self, other_vec : Vec2):
        return self.x == other_vec.x and self.y == other_vec.y
    
    def __add__(self, other_vec : Vec2):
        return self.add(other_vec)
    
    def __mul__(self, scalar):
        return self.scale(scalar)
    
    def __rmul__(self, scalar):
        return self.scale(scalar)

2 * Vec2(2, 0) + 5 * Vec2(0, 2)

2d 벡터 repr 추가

2d 벡터 출력시 원하는 형태로 안나오고 주소가 나옴

repr 작성해서 잘 나오도록 수정

class Vec2():
    def __init__(self, x, y):
        self.x = x
        self.y = y
    
    def add(self, v2):
        return Vec2(self.x + v2.x, self.y + v2.y)

    def scale(self, scalar):
        return Vec2(scalar * self.x, scalar * self.y)
    
    def __eq__(self, other_vec : Vec2):
        return self.x == other_vec.x and self.y == other_vec.y
    
    def __add__(self, other_vec : Vec2):
        return self.add(other_vec)
    
    def __mul__(self, scalar):
        return self.scale(scalar)
    
    def __rmul__(self, scalar):
        return self.scale(scalar)
    
    def __repr__(self):
        return "Vec2({},{})".format(self.x, self.y)

2 * Vec2(2, 0) + 5 * Vec2(0, 2)

3차원 벡터 정의하기

2d 때랑 거의 비슷, z만 추가

class Vec3():
    def __init__(self, x, y, z):
        self.x = x
        self.y = y
        self.z = z
    
    def add(self, v2):
        return Vec3(self.x + v2.x, self.y + v2.y, self.z + v2.z)
    
    def scale(self, scalar):
        return Vec3(self.x * scalar, self.y * scalar, self.z * scalar)

    def __add__(self, v2):
        return self.add(v2)
    
    def __eq__(self, v2):
        return self.x == v2.x and self.y == v2.y and self.z == v2.z
    
    def __mul__(self, scalar):
        return self.scale(scalar)

    def __rmul__(self, scalar):
        return self.scale(scalar)
    
    def __repr__(self):
        return "Vec3({},{},{})".format(self.x, self.y, self.z)

2.0 * (Vec3(1, 0, 0) + Vec3(0, 1, 0))

메타 클래스 Vector를 만들고, 상속받는 클래스 vec2 정의하기

메타 클래스 Vector에선 vec2. vec3 둘다 __mul__, __add__ 같은 연산자 함수 포함되므로 다음과 같이 정의

실제 add, scale 수행되는 과정은 vec2와 vec3은 다르므로 추상 메서드로 비워둠

추가로 차연산은 + (우측 값 x -1)을 하면 되므로 구현 

메타 클래스 Vector를 상속받는 Vec2를 구현하면서

산술 함수는 이미 Vector에서 구현했으므로 안해도되고,

실제 산술 동작 add, scale과 기타 필요한 함수들 구현해주면 됨 

from abc import ABCMeta, abstractmethod

class Vector(metaclass=ABCMeta):
    @abstractmethod
    def scale(self, scalar):
        pass
    @abstractmethod
    def add(self, toher):
        pass
    
    def __mul__(self, scalar):
        return self.scale(scalar)
    
    def __rmul__(self, scalar):
        return self.scale(scalar)
    
    def __add__(self, other):
        return self.add(other)
    
    def subtract(self, other):
        return self.add(-1 * other)
    
    def __sub__(self, other):
        return self.subtract(other)

class Vec2(Vector):
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def add(self, other):
        return Vec2(self.x + other.x, self.y + other.y)
    
    def scale(self, scalar):
        return Vec2(self.x * scalar, self.y * scalar)
    
    def __eq__(self, other):
        return self.x == other.x and self.y == other.y and self.z == other.z
    
    def __repr__(self):
        return "Vec2({},{})".format(self.x, self.y)
    
Vec2(5, 5) - Vec2(1, 3)

메타클래스로 vec3 만들기

구현 과정은 방금 한 vec2랑 같음, z 내용만 추가

class Vec3(Vector):
    def __init__(self, x, y, z):
        self.x = x
        self.y = y
        self.z = z
    
    def add(self, other):
        return Vec3(self.x + other.x, self.y + other.y, self.z + other.z)
    
    def scale(self, scalar):
        return Vec3(self.x * scalar, self.y * scalar, self.z * scalar)
    
    def __eq__(self, other):
        return self.x == other.x and self.y == other.y and self.z == other.z

    def __repr__(self):
        return "Vec3({},{},{})".format(self.x, self.y, self.z)

vec3 = Vec3(3, 2, 1) * 3
print(vec3)

좌표 벡터 만들기

지금까지 2, 3차원 벡터 다루는 클래스를 만들었지만 상위 차원을 다룰수 있도록 좌표 벡터를 만들면

__init__ coordinates를 가변인자로 받아 저장

add와 scale도 구현. 실제 내용은 vectors.py에 이런식으로 구현되있음

def add(*vectors):
    return tuple(map(sum,zip(*vectors)))

def scale(scalar,v):
    return tuple(scalar * coord for coord in v)

repr도 출력 형태도 구현

추상 프로퍼티와 추상 메서드 차이는 프로퍼티는 값. 추상 메서드는 동작 의미

from abc import abstractproperty
from vectors import add, scale

class CoordinateVector(Vector):
    @abstractproperty
    def dimension(self):
        pass

    def __init__(self, *coordinates):
        self.coordinates = tuple(x for x in coordinates)
    
    def add(self, other):
        return self.__class__(*add(self.coordinates, other.coordinates))
    
    def scale(self, scalar):
        return self.__class__(*scale(scalar, self.coordinates))
    
    def __repr__(self):
        return "{}{}".format(self.__class__.__qualname__, self.coordinates)

좌표벡터 상속받는 vec6만들기

이미 좌표 벡터에서 내용 대부분 구현되있어서 이건 간단

추상 프로퍼티로 설정한 dimension만 구현해주면 된다.

print vec6하면 repr대로 출력

print vec6.coordinates는 튜플

print vec6.dimension은 인트값 6

덧셈 연산도 vectors.py에 구현된 add 대로 잘 동작

class Vec6(CoordinateVector):
    def dimension(self):
        return 6

print(Vec6(1, 1, 1, 1, 4, 5))
print(Vec6(1, 1, 1, 1, 4, 5).coordinates)
print(Vec6(1, 1, 1, 1, 4, 5).dimension())
Vec6(1, 1, 1, 1, 4, 5) + Vec6(1, 2, 3, 4, 5, 6)

영벡터를 추상 프로퍼티로 갖는 메타 클래스 벡터 만들기

지난번에 만든 벡터에 추상 프로퍼티로 추가해주면 됨

추가로 scale을 이용한 neg 연산도 추가

class Vector(metaclass=ABCMeta):
    @abstractmethod
    def scale(self, scalar):
        pass

    @abstractmethod
    def add(self, other):
        pass

    def __add__(self, other):
        return self.add(other)
    
    def __mul__(self, scalar):
        return self.scale(scalar)
    
    def __rmul__(self, scalar):
        return self.scale(scalar)
    
    def substract(self, other):
        return self.add(-1 * other)
    
    def __sub__(self, other):
        return self.substract(other)

    @classmethod
    @abstractproperty
    def zero():
        pass

    def __neg__(self):
        return self.scale(-1)

추상 프로퍼티 zero를 구현한 Vec2 만들기

기존의 vec2에다가 추상 프로퍼티 zero(영벡터 반환)를 구현

메타클레스에서 __neg__ 구현한 덕분에 -Vec2 연산이 가능해짐

class Vec2(Vector):
    def __init__(self, x, y):
        self.x = x
        self.y = y
    
    def scale(self, scalar):
        return Vec2(self.x * scalar, self.y * scalar)
    
    def add(self, other):
        return Vec2(self.x + other.x, self.y + other.y)
    
    def zero():
        return Vec2(0, 0)
    
    def __repr__(self):
        return "Vec2({},{})".format(self.x, self.y)
    
    def __eq__(self, other):
        return self.x == other.x and self.y == other.y

print(Vec2.zero())
-Vec2(3, -4)

덧셈, 같은지 연산 전에 클래스 타입 확인하기

이전에 구현한 코드들은 vec2(1, 2) == vec3(1, 2, 3)은 같다로 반환하므로 

덧셈, 같음 연산전에 클래스 타입 확인하도록 수정이 필요

add와 eq 리턴 전에 assert 해주면 됨

같은 클래스면 문제없이 동작하지만 다르면 assertion error 발생

class Vec2(Vector):
    def __init__(self, x, y):
        self.x = x
        self.y = y
    
    def add(self, other):
        assert self.__class__ == other.__class__
        return Vec2(self.x + other.x, self.y + other.y)
    
    def scale(self, scalar):
        return Vec2(self.x * scalar, self.y *scalar)
    
    def __eq__(self, other):
        assert self.__class__ == other.__class__
        return self.x == other.x and self.y == other.y
    
    def __repr__(self):
        return "Vec2({},{})".format(self.x, self.y)
    
    def zero(self):
        return Vec2(0, 0)
 
print(Vec2(1, 2) == Vec2(1, 1))
print(Vec2(1, 2) + Vec2(3, 2))

Vec2(1, 2) == Vec3(1, 1, 1)
Vec2(1, 2) + Vec3(1, 1, 1)

메타클레스 벡터에 나누기 추가하기

추상 메서드로 덧셈, 스케일

추상 프로퍼티로 제로

그외는 덧셈, 그냥곱, 우측곱, 빼기, 부정 구현

추가된 나누기는 1.0/scalar를 곱하는 형태로 구현

class Vector(metaclass=ABCMeta):
    @abstractmethod
    def add(self, other):
        pass

    @abstractmethod
    def scale(self, scalar):
        pass

    def __add__(self, other):
        return self.add(other)

    def __mul__(self, scalar):
        return self.scale(scalar)
    
    def __rmul__(self, scalar):
        return self.scale(scalar)
    
    def substract(self, other):
        return self.add(-1 * other)
    
    def __sub__(self, other):
        return self.substract(other)

    @classmethod
    @abstractproperty
    def zero():
        pass

    def __neg__(self):
        return self.scale(-1)

    def __truediv__(self, scalar):
        return self.scale(1.0/scalar)

나눗셈 추가된 vec2 나누기 연산

나눗셈 연산은 vector에 추가되서 vec2에선 고칠게 없음

그냥 그대로 나눠보면 결과나옴

class Vec2(Vector):
    def __init__(self, x, y):
        self.x = x
        self.y = y
    
    def add(self, other):
        assert self.__class__ == other.__class__
        return Vec2(self.x + other.x and self.y + other.y)

    def __eq__(self, other):
        assert self.__class__ == other.__class__
        return self.x == other.x and self.y == other.y
    
    def scale(self, scalar):
        return Vec2(scalar * self.x, scalar * self.y)
    
    def zero():
        return Vec2(0, 0)
    
    def __repr__(self):
        return "Vec2({},{})".format(self.x, self.y)

print(Vec2(2, 3)/2)

1차원 벡터 정의하기

2차원 벡터때보다 더단순하다

class Vec1(Vector):
    def __init__(self, x):
        self.x = x
    
    def add(self, other):
        assert self.__class__ == other.__class__
        return Vec1(self.x + other.x)
    
    def __eq__(self, other):
        assert self.__class__ == other.__class__
        return self.x == other.x
    
    def scale(self, scalar):
        return Vec1(self.x * scalar)

    @classmethod
    def zero(cls):
        return Vec1(0)
    
    def __repr__(self):
        return "Vec1({})".format(self.x)

vec1_1 = Vec1(3)
vec1_2 = Vec1(5)
print(vec1_1)
print(vec1_1 * 3)
print(-vec1_2)
print(vec1_1/4)
print(Vec1.zero())

0벡터 다루기

담는 값이 없어서 더 단순

class Vec0(Vector):
    def __init__(self):
        pass

    def add(self, other):
        return Vec0()
    
    def scale(self, scalar):
        return Vec0()
    
    def __eq__(self, other):
        return self.__class__ == other.__class__ == Vec0
    
    def __repr__(self):
        return "Vec0()"
    
    @classmethod
    def zero(cls):
        return Vec0()

print(-3.14 * Vec0())
print(Vec0() + Vec0() + Vec0())

이거 원래 numpy로 간단하게 하던건데 직접 구현해서 쓰자니

행렬 그냥 적어서 풀때 비해서 잘 이해가는 편은 아님

행렬과 벡터의 곱

이 코드에서 zip(matrix) 하면 기존 matrix 전치행렬 나오고

선형 결합 함수에서는 스칼라곱해서 더하면 원래 벡터 x 행렬 결과가 잘 나온다.

대신 로직은 좀 이해안가는데 대충 넘어가야할듯

from vectors import *

def linaer_combination(scalars, *vectors):
    scaled = [scale(s, v) for s, v in zip(scalars, vectors)]
    return add(*scaled)

def multiply_matrix_vector(matrix, vector):
    return linear_combination(vector, *zip(*matrix))

B = (
    (0, 2, 1),
    (0, 1, 0),
    (1, 0, -1)
)
v = (3, -2, 5)

multiply_matrix_vector(B, v)

행렬간 곱셈 연산

그나마 행렬 곱 연산 로직은 벡터 x 행렬보단 이해는 간다

대충 위 a, b 곱 연산 흐름을 적어보면 이런듯

from vectors import *

def matrix_multiply(a, b):
    return tuple(
        tuple(dot(row, col) for col in zip(*b)) for row in a
    )

a = (
    (1, 1, 0),
    (1, 0, 1),
    (1, -1, 1)
)

b = (
    (0, 2, 1),
    (0, 1, 0),
    (1, 0, -1)
)

matrix_multiply(a, b)

주전자 회전시키기

animate_teapot.py 코드를 보면 안에 시간에 따라 z축 회전 행렬을 반환하는 함수와

pygame으로 주전자 모델 드로잉 하는 함수가 제공되는데

여기다 정리하자니 너무 복잡해질것같고 그냥 캡처

from teapot import load_triangles
from draw_model import draw_model
from math import sin,cos

def get_rotation_matrix(t): #1
    seconds = t/1000 #2
    return (
        (cos(seconds),0,-sin(seconds)),
        (0,1,0),
        (sin(seconds),0,cos(seconds))
    )

####################################################################
#### this code takes a snapshot to reproduce the exact figure 
#### shown in the book as an image saved in the "figs" directory
#### to run it, run this script with command line arg --snapshot
import sys
import camera
if '--snapshot' in sys.argv:
    camera.default_camera = camera.Camera('fig_5.4_draw_teapot',[0,1000,2000,3000,4000])
####################################################################

draw_model(load_triangles(), get_matrix=get_rotation_matrix)

3차원 공간에 공룡 그리기

이전에 2차원에 공룡 그린것과 큰 차이는 없고

기존의 공룡 벡터에다가 z축을 1로 추가해서 드로잉

from draw3d import draw3d, Points3D, Segment3D

dino_vectors = [(6,4), (3,1), (1,2), (-1,5), (-2,5), (-3,4), (-4,4),
    (-5,3), (-5,2), (-2,2), (-5,1), (-4,0), (-2,1), (-1,0), (0,-3),
    (-1,-4), (1,-4), (2,-3), (1,-2), (3,-1), (5,1)
]

def polygon_segments_3d(points, color="blue"):
    count = len(points)
    return [Segment3D(points[i], points[(i+1)%count], color=color) for i in range(0, count)]

dino_3d = [(x, y, 1) for x, y in dino_vectors]

draw3d(
    Points3D(*dino_3d, color="blue"),
    *polygon_segments_3d(dino_3d)
)

행렬로 공룡 평행이동 시키기

평행이동 행렬은 좌표 벡터와 곱하여 dx, dy값이 반영됨

translate_matrix = (
    (1, 0, 5),
    (0, 1, 2),
    (0, 0, 1)
)

translated_dino = [multiply_matrix_vector(translate_matrix, v) for v in dino_3d]

draw3d(
    Points3D(*dino_3d, color="C0"),
    *polygon_segments_3d(dino_3d, color="C0"),
    Points3D(*translated_dino, color="C3"),
    *polygon_segments_3d(translated_dino, color="C3")
)

행렬로 회전 편행이동 시키기

이전 코드와 거의 비슷한데 행렬에 xy 회전할수 있도록 값을 변경하면 됨

2차원 회전 행렬은 이런 식, 90도 회전하는경우 

cos 90 = 0, sin 90 = 1

R = [

       [ 0  -1],

       [1 0]

]

90도 회전하고, (5, 1)만큼 회전하는 경우 이런 식이 됨

rotate_and_translate_matrix = (
    (0, -1, 5),
    (1, 0, 1),
    (0, 0, 1)
)

rot_trans_dino = [multiply_matrix_vector(rotate_and_translate_matrix, v) for v in dino_3d]

draw3d(
    Points3D(*dino_3d, color="C0"),
    *polygon_segments_3d(dino_3d, color="C0"),
    Points3D(*rot_trans_dino, color="C3"),
    *polygon_segments_3d(rot_trans_dino, color="C3")
)

곤란한게 여기 3d 예제 대부분은 pygame과 주전자 모델 사용해서 직접 구현하기가 곤란하다.

아쉬운데로 8면체를 갖고 적용해보면

8면체로 x, y, z축 평행이동

이전에 3차원 8면체 드로잉한걸 조금 수정함

아깐 3차원이긴 하나 z=1로 고정된 공룡을 xy축으로 평행, z로 회전하다보니 변환행렬은 (3x3)로 충분했지만

xyz 축 전체 회전/평행이동을 위해선 각 점 길이가 4인 벡터 (x, y, z, 1)를 쓰며 변환행렬은 4 x 4 형태이다.

(0, 0, 0)을중심으로 한 8면체를 (5, 3, 2)로 이동한 결과

top = (0, 0, 3, 1)
bottom = (0, 0, -3, 1)
xy_plane = [(1, 0, 0, 1), (0, 1, 0, 1), (-1 , 0, 0, 1), (0, -1, 0, 1)]
edges = [Segment3D(top[:-1], p) for p in xy_plane] +\
    [Segment3D(bottom[:-1], p) for p in xy_plane] +\
    [Segment3D(xy_plane[i][:-1], xy_plane[(i + 1) % 4][:-1]) for i in range(0, 4)]
trans_matrix = (
    (1, 0, 0, 5),
    (0, 1, 0, 3),
    (0, 0, 1, 2),
    (0, 0, 0, 1)
)
def get_trans_octa(top, bottom, xy_plane, trans_matrix, color="red"):
    trans_top = multiply_matrix_vector(trans_matrix, top)[:-1]
    trans_bottom = multiply_matrix_vector(trans_matrix, bottom)[:-1]
    trans_xy_plane = [multiply_matrix_vector(trans_matrix, v)[:-1] for v in xy_plane]
    trans_edges = [Segment3D(trans_top, p, color=color) for p in trans_xy_plane] +\
        [Segment3D(trans_bottom, p, color=color) for p in trans_xy_plane] +\
        [Segment3D(trans_xy_plane[i], trans_xy_plane[(i + 1) % 4], color=color) for i in range(0, 4)]
    return trans_edges
octa_edges = get_trans_octa(top, bottom, xy_plane, trans_matrix)
draw3d(
    *edges,
    *octa_edges
)

8면체 회전, 평행이동

방금 평행이동 했으니 이번엔 추가로 회전

아깐 공룡으로 z축 90도 회전했는데 이번엔 y축으로 45도 회전시켜보자

cos(45) = 0.52

sin(45) = 0.85

y축으로 45도보다 더 기울어진것 같긴한데 대충 회전한건 맞으니 넘어가자

이전 글에서는 간단하게 점과 폴리곤을 그리고, 평행/회전이동, matplotlib을 이용한 드로잉 모듈을 만들었음

이번에는 3d 차원을 다뤄보자

3d 공간 띄우기

draw 3d 모듈 가져와서 실행하면 이런식으로 matplotlib에서 3d로 만든 공간이 나옴

from draw3d import *
draw3d()

점과 선분 그리기

여기서는 3차원 선분을 Segment3D로 제공

draw3d(
    Points3D((2, 2, 2), (1, -2, -2)),
    Segment3D((2, 2, 2), (1, -2, -2))
)

3차원 박스 그리기

box 3d를 쓰면 원점을 한 모서리로하는 3차원 박스 그릴수있음

draw3d(
    Points3D((2, 2, 2), (1, -2, -2)),
    Segment3D((2, 2, 2), (1, -2, -2)),
    Box3D(2, 2, 2),
    Box3D(1, -2, -2)
)

점과 선분으로 3d 사각형 그리기

이번엔 3차원 점과 3차원 선분으로 그림

pm1 = [1, -1]
vertices = [(x, y, z) for x in pm1 for y in pm1 for z in pm1]
print(vertices)

edges = [((-1,y,z),(1,y,z)) for y in pm1 for z in pm1] +\
                [((x,-1,z),(x,1,z)) for x in pm1 for z in pm1] +\
                [((x,y,-1),(x,y,1)) for x in pm1 for y in pm1]
print(edges)

draw3d(
    Points3D(*vertices),
    *[Segment3D(*edge) for edge in edges]
)

세 백터의 합 구하기, 3차원에서 길이 계산

zip은 리스트의 같은 번째 요소들끼리 묶음

(1, 1, 3), (2, 4, -4), (4, 2, -2) 0번 끼리, 1번끼리 2번끼리 묶어 (1, 2, 4), (1, 4,2). 식으로 됨.

길이야 2d 때처럼 유클리디안 거리로 구현

def add(*vecs):
    by_coord = zip(*vecs)
    coord_sum = [sum(coords) for coords in by_coord]
    return tuple(coord_sum)

print(list(zip(*[(1, 1, 3), (2, 4, -4), (4, 2, -2)])))
print([sum(coords) for coords in [(1, 2, 4), (1, 4, 2), (3, -4, -2)]])


def add(*vecs):
    return tuple(map(sum, zip(*vecs)))
add((1, 1, 3), (2, 4, -4), (4, 2, -2))



from math import sqrt
def length(v):
    return sqrt(sum([coord **2 for coord in v]))
length((3, 4, 5))

3차원 벡터 합 드로잉 하기

draw3d에 arrow3d로 드로잉

(4, 0, 3)과 (-1, 0, 1)을 합하여 (3, 0, 4) 가 나옴

draw3d(
    Arrow3D((4, 0, 3), color=red),
    Arrow3D((-1, 0, 1), color=blue),
    Arrow3D((3, 0, 4), color=purple)
)

삼각함수로 3차원 회전상승? 그리기

pi/6 만큼 상 이동하는 (sin, cos, 1/3) 벡터 리스트 vs 준비

벡터 합 연산을 누적 시켜 arrow 3d 준비

draw arrows로 회전상승하는 화살표 그리기

from math import sin, cos, pi
vs = [(sin(pi*t/6), cos(pi*t/6), 1.0/3) for t in range(0, 24)]

running_sum = (0, 0, 0)
arrows = []
for v in vs:
    next_sum = add(running_sum, v)
    arrows.append(Arrow3D(next_sum, running_sum))
    running_sum = next_sum

print(running_sum)
draw3d(*arrows)

3d 벡터 스칼라 곱 구현, 길이가 1로 조정하기

스칼라 곱으로 스케일 함수 구현

-1, -1, 2 유클리디안 거리 구하기

1/vec_len으로 몇을 곱하면 1이 되는지 찾기

(-1, -1, 2) 벡터에 0.408 스칼라곱하여 길이가 1이되도록 조정

def scale(scalar, v):
    return tuple(scalar * coord for coord in v)

vec_len = length((-1, -1, 2))
s = 1/vec_len
vec_scaled = scale(s, (-1, -1, 2))

print(vec_len)
print(s)
print(vec_scaled)

내적

- 벡터간 곱셈 연산 -> 스칼라 반환, 점곱/스칼라곱 이라고도함

- 두 벡터끼리 얼마나 같은 방향인지 나타내는 정도

- 내적 = u dot v = |u| x |v| x cos(theta)

 * 두 벡터 사이 각이 90인 경우 0 이 됨 (cos 90 = 0이므로)

def dot(u, v):
    return sum([coord1 * coord2 for coord1, coord2 in zip(u, v)])

print(dot((1, 0), (0, 2)))
print(dot((0, 3, 0), (0, 0, -5)))
print(dot((3, 4), (2, 3)))

두 벡터 사이간 각도

- 내적 공식 u dot v = |u| |v| cos theta 을 정리하면

 cos theta = u dot v / (|u| |v|)

theta = acos(u dot v / (|u| |v|) ) 로 정리 가능

벡터 (3, 4)와 (2, 3)의 내적 결과 18

두 벡터 사이 각도 0.0554 rad

|(3, 4)| x |(2, 3)| x cos(0.0554 rad) = 18

from math import cos, pi, acos

def angle_between(v1, v2):
    return acos(dot(v1, v2) / (length(v1) * length(v2)))

ang = angle_between((3,4), (2, 3))
print(ang)

print(dot((3, 4), (2, 3)))
dot_product = length((3, 4)) * length((2, 3)) * cos(ang)
print(dot_product)

외적

- 벡터 x 벡터 -> 벡터, 벡터곱이라고 부르며 벡터를 반환해서 크기와 방향을 가짐

- 방향은 두 벡터에 수직( 오른손 법칙)

- 크기는 두 벡터가 얼마나 포개지는지 다룬 내적과는 달리 얼마나 90도에 가까운지에 달림

=>  u x v = u cross v = |u| |v| sin theta

                 =  (uy * vz - uz * vy, uz*vx - ux*vz, ux*vy - uy*vx)

def cross(u, v):
    ux, uy, uz = u
    vx, vy, vz = v
    return (uy * vz - uz * vy, uz*vx - ux*vz, ux*vy - uy*vx)

print(cross((1, 0, 0), (0, 1, 0)))
print(cross((1, 0, 0), (0, 2, 0)))

외적 계산 예시보기

 외적은 두 벡터에 수직인 방향을 가지며, 크기는 90도일수록 큰 값을 가짐

(1, 0, 0)과 (0, 0, 1)의 외적 결과 (0, 0, 1)은 z축에 나란한 형태가 나옴

(1, 1, 1)과 (-1, 1, 1)의 외적결과도 두 벡터간 수직 방향 향함

두 벡터 다 크기가 같으나

두 벡터사이 각이 좁을떄는 작고

90도 되니 좁을때보다 커짐

3차원 벡터의 2차원 투영

3차원 벡터를 x, y 축 단위 벡터랑 점곱 연산해서 구함

z축 값이 0이 된다.

# (1, 0, 0)과 dot 연산
# (0, 1, 0)과 dot 연산
def component(v, direction):
    return (dot(v, direction) / length(direction))

# 3차원 백터를 xy 평면(2차원)에 투영하는 함수
def vector_to_2d(v):
    return (component(v, (1, 0, 0)), component(v, (0, 1, 0)))

u = (3, 2, 5)
projected_vec_2d = vector_to_2d(u)
projected_ved_3d = (projected_vec_2d[0], projected_vec_2d[1], 0)
print(u)
print(projected_ved_3d)
draw3d(
    Arrow3D(u),
    Arrow3D(projected_ved_3d, color="blue")
)

3차원 8면체 그리기

단순히 8면체 꼭지점 준비해서 연결만 하면 끝

top = (0, 0, 1)
bottom = (0, 0, -1)
xy_plane = [(1, 0, 0,), (0, 1, 0), (-1 , 0, 0), (0, -1, 0)]
edges = [Segment3D(top, p) for p in xy_plane] +\
    [Segment3D(bottom, p) for p in xy_plane] +\
    [Segment3D(xy_plane[i], xy_plane[(i + 1) % 4]) for i in range(0, 4)]
draw3d(*edges)

8면체 2차원 투영 준비하기

8면체의 면 face들은 세 점으로 구성

8면체 변수는 8개의 면 리스트

normal로 면에 수직 벡터 취득 후 unit으로 단위벡터화

unit_normal[2] > 0는 면의 수직 벡터가 위를 향함을 의미. 위에서 볼것이므로 윗 면만 나와야 함

face_to_2d로 3차원 면을 2차원 폴리곤으로 변환 후 전체 드로잉

from vectors import *
from draw2d import *

# 면에서 정점 가져옴
def vertices(faces):
    return list(set([vertex for face in faces for vertex in face]))

# 면을 2차원으로 투영
def face_to_2d(face):
    return [vector_to_2d(vertex) for vertex in face]

# 길이가 1인 단위벡터화
def unit(v):
    return scale(1./length(v), v)

# 두벡터간 차연산
def subtract(v1,v2):
    return tuple(v1-v2 for (v1,v2) in zip(v1,v2))

def normal(face):
    return cross(subtract(face[1], face[0]), subtract(face[2], face[0]))

blues = matplotlib.cm.get_cmap("Blues")
def render(faces, light = (1, 2, 3), color_map=blues, lines=None):
    polygons = []
    for face in faces:
        unit_normal = unit(normal(face))
        if unit_normal[2] > 0:
            c = color_map(1 - dot(unit(normal(face)), unit(light)))
            p = Polygon2D(*face_to_2d(face), fill=c, color=lines)
            polygons.append(p)
    draw2d(*polygons, axes=False, origin=False, grid=None)

octahedron = [
    [(1,0,0), (0,1,0), (0,0,1)],
    [(1,0,0), (0,0,-1), (0,1,0)],
    [(1,0,0), (0,0,1), (0,-1,0)],
    [(1,0,0), (0,-1,0), (0,0,-1)],
    [(-1,0,0), (0,0,1), (0,1,0)],
    [(-1,0,0), (0,1,0), (0,0,-1)],
    [(-1,0,0), (0,-1,0), (0,0,1)],
    [(-1,0,0), (0,0,-1), (0,-1,0)],
]

8면체 2차원에 투영하기

render(octahedron, color_map=matplotlib.cm.get_cmap('Blues'), lines=black)

파이썬수학 글에서 다루는 코드들은 개발자를 위한 수학 깃헙 베이스로 함.

https://github.com/orlandpm/Math-for-Programmers

공룡 그리기(점)

- 벡터리스트는 공룡 모양을 나타내는 점들이지만 선을 안그려서 공룡처럼 보이진 않음

from vector_drawing import *

dino_vecs = [
    (6,  4), (3, 1), (1, 2), (-1, 5), (-2, 5), (-3, 4), (-4, 4),
    (-5, 3), (-5, 2), (-2, 2), (-5, 1), (-4, 0), (-2, 1), (-1, 0), (0, -3),
    (-1, -4), (1, -4), (2, -3), (1, -2), (3, -1), (5, 1)
]

print(dino_vecs)
print(*dino_vecs)
draw(Points(*dino_vecs))

*의 의미

 그렇게 파이썬 쓰면서 어려운걸 안하니 변수앞에 *의미가 뭔지도 모르고, 포인턴가 하고 있었는데 아래 글에서 의미를 찾아봄. 가변인자 매개변수를 나타낸다고 한다.

https://hcnoh.github.io/2019-01-27-python-arguments-asterisk

이번엔 선하나 긋기

 코드 작성자가 draw 함수안에 matplotlib으로 다 구현해놔서, 골치아픈걸 신경안쓰고 이렇게 코드만 넣으면 동작되는데, 이번에는 points 뒤에 segment로 선 그을 두 좌표 넣으면 됨

draw(
    Points(*dino_vecs),
    Segment((6, 4), (3, 1))
)

폴리곤으로 공룡 그리기

 방금 segment 함수로 선하나 그렸지만 좌표 벡터를 그대로 폴리곤으로 그려주도록 기능이 이미 있다. 대충 polygon쓰면 공룡이 그려진다.

draw(
    Points(*dino_vecs),
    Polygon(*dino_vecs)
)

포물선 그리기

 x가 -5, 5인 구간에서 y = x^2인 포물선을 그려보자

 깃헙 자료에는 리스트 컴프리헨션으로 한줄에 점들 구하지만, 난 잘 못해서 그냥 풀어씀

거기다 선도 그려주려고 폴리곤 썻는데, 포물선이 닫히고 말았지만 그냥 넘어가자

vecs = []
for x in range(-5, 6):
    y = pow(x, 2)
    vecs.append((x, y))
print(vecs)
draw(
    Points(*vecs),
    Polygon(*vecs),
    grid=(1, 10),
    nice_aspect_ratio=False
)

벡터 합 연산으로 공룡을 평행 이동시키자

 두 벡터간 합 연산하는 함수를 정의하고

 기존 공룡 벡터리스트를 -2, -5 한뒤

 기존 공룡과 새 공룡 벡터로 드로잉하면 이렇게 나온다.

def add(v1, v2):
    return (v1[0] + v2[0], v1[1] + v2[1])
    
dino_vecs2 = []
for vec in dino_vecs:
    new_vec = add((-2, -5), vec)
    dino_vecs2.append(new_vec)

draw(
    Points(*dino_vecs),
    Polygon(*dino_vecs),
    Points(*dino_vecs2, color="red"),
    Polygon(*dino_vecs2, color="red")
)

두점 사이 거리 계산하기

 피타고라스 정리로 배운 두 점사이 거리 계산하는 코드

from math import sqrt

def length(v):
    return sqrt(v[0]**2 + v[1]**2)

여러개 좌표 벡터 더하기

아까 구현한 덧셈 연산은 좌표 벡터가 2개만 들어온경우를 생각하는데

이번엔 2개가 아니라 여러개 들어온 경우 덧셈하는 코드를 만들어보자

한줄 쓰기는 좋아하진 않지만 그냥 만든 사람 따라해야지

def add(*vecs):
    return (sum(v[0] for v in vecs), sum(v[1] for v in vecs))

여러개 평행이동 연산 하기

아까 구현한 공룡 이동은 -2, -5 한번 하고 끝냈는데

(-1, -1), (-2, -2), (-3, -3) 처럼 좌표 벡터 여러개 주면 여러개 연산 결과 나오게 하자

=> 3,3 에서  각각 순서대로 연산하면 이런 식으로 나온다

def translate(start, vecs):
    return [add(start, vec) for vec in vecs]
translate((3, 3), [(-1, -1), (-2, -2), (-3, -3)])

공룡 백마리 그리기

방금 구현한 평행이동 코드로 공룡 100마리 그려보자.

아깐 길이가 3인 평행이동 벡터리스트 만들었지만

길이가 100인 평행이동 벡터리스트 만들고

평행이동 함수로 공룡 폴리곤 100개 만들어 draw에 집어넣으면 됨

def hundred_dinos():
    translations = [
        (12 * x, 10 * y) for x in range(-5, 5) for y in range(-5, 5)
    ]

    dinos = [Polygon(*translate(t, dino_vecs)) for t in translations]
    draw(*dinos, grid=None, axes=None, origin=None)
hundred_dinos()

네모 그리고, 네모 폴리곤 길이 구하기

1. 두 점 벡터 차연산

2. 두 점 사이 거리 계산

3. 벡터들 들어오면 모든 벡터 사이 길이 계산 구현

def subtract(v1, v2):
    return (v1[0] - v2[0], v1[1] - v2[1])

def distance(v1, v2):
    #length는 피타고라스 정리로 계산했던 대각 길이
    return length(subtract(v1, v2))

def perimeter(vecs):
    dists = [distance(vecs[i], vecs[(i + 1) % len(vecs)]) for i in range(0, len(vecs))]
    return sum(dists)

square = [(0, 0), (0, 2), (1, 2), (1, 0)]
print(perimeter(square))
draw(Polygon(*square), grid=(0.5, 0.5), axes=None)

지금까지 직교좌표 다뤘고 이젠 극좌표 다루자

데카르트 좌표(카티지안) : (x, y)

극좌표 : (반지름 r, 각도 theta)

극좌표 -> 직교좌표

여기선 반지름 10, 45도(라디안 pi /4)로 직교 좌표로 변환

from math import tan, pi, sin, cos

def to_cartesian(polar_vec):
    len, ang = polar_vec[0], polar_vec[1]
    return (len * cos(ang), len * sin(ang))

angle = pi / 4
to_cartesian((10, angle))

tan, atan, 직교->극좌표 변환

tan(라디안) 시 기울기 반환

atan(y/x) 시 라디안 반환

atan2(y, x)시 라디안 반환

직교 -> 극좌표 바꾸는 to_polar() : 유클리디안거리(피타고라스정리), atan2로 극좌표 반환

극좌표 -> 직교 변환 및 드로잉

len : cos( 5 * 0 ~ 2pi)

ang : 0 ~ 2pi

사이 1000개 극좌표를 직교좌표로 변환하여 폴리곤 드로잉

polar_coords = [(cos(5 * x * pi / 500.0), 2 * pi * x / 1000.0) for x in range(0, 1000)]
#to_cartesian(p) p = (len, ang)
vecs = [to_cartesian(p) for p in polar_coords]
draw(Polygon(*vecs))

atan 드로잉 해보기

2, -3 극좌표 변환시 len 3.6, -0.98 rad(=-56 deg)

유클리디안 거리 2, -3 시 3.6

atan(-3/2)시 -0.98 rad

print(to_polar((2, -3))) #-0.98 rad = -56
print(length((2, -3)))
print(atan(-3/2)) # y/x
draw(Arrow((2, -3)), Points((2, -3)))

직교 -> 극좌표 -> 회전 덧셈 -> 직교 -> 드로잉 순으로 공룡 회전시키기

음.. 내가 알기론 직교좌표계에서 회전할때 회전행렬 사용했는데

여기선 회전행렬 얘기를 아직 안다루어서인지

직교 -> 극좌표 변환

극좌표(거리, 각도)의 각도 + 회전각 => 극좌표(거리, 원래각도 + 회전각)

바뀐극좌표 -> 직교 변환 => 폴리곤 드로잉

순으로 구현

rot_ang = pi / 4
dino_polar = [to_polar(v) for v in dino_vecs]
dino_rot_polar = [ (len, ang + rot_ang) for len, ang in dino_polar]
dino_rot_cart = [to_cartesian(p) for p in dino_rot_polar]

draw(
    Polygon(*dino_vecs, color=gray),
    Polygon(*dino_rot_cart, color=red)
)

공룡 회전 + 평행이동시키기

rot -> translate 하여 공룡 회전, 이동 구현

def rotate(rot_ang, vecs):
    polars = [to_polar(v) for v in vecs]
    return [to_cartesian((len, ang + rot_ang)) for len, ang in polars]

# 5 * pi / 3 = 5.2 rad = 297 deg
new_dino_vecs_notmove = translate((0, 0), rotate(5 * pi / 3, dino_vecs))
new_dino_vecs = translate((5, 5), rotate(5 * pi / 3, dino_vecs))
draw(
    Polygon(*dino_vecs),
    Polygon(*new_dino_vecs_notmove, color="red"),
    Polygon(*new_dino_vecs, color="green")
)

지금까지 사용한 2차원 드로잉 모듈 만들기

지금까지 구현한 코드는 원래 matplotlib으로 그릴려고하면 엄청 고생해야하는데

이 저장소 만드사람이 드로잉 하기 쉽게 다 만들어 두어서 쉽게 그릴수가 있었다.

하지만 다른데 쓰려면 직접 만들줄 알아야 하므로 이 드로잉 모듈을 따라 만들어봄

물론 이 분만큼 다 고려해서 만들기도 힘들고 금방 만들수 있게

자주 쓴것과 꼭 필요한 부분만 커스텀으로 따로 만들었음

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.pyplot import xlim, ylim

blue = 'C0'
black = 'k'
red = 'C3'
gray = 'gray'

class Polygon():
    def __init__(self, *vertices, color=blue, fill=None, alpha=0.4):
        self.vertices = vertices
        self.color = color
        self.fill = fill
        self.alpha = alpha

class Points():
    def __init__(self, *vectors, color=black):
        self.vectors = list(vectors)
        self.color = color

# helper function to extract all the vectors from list of objs
def extract_vectors(objs):
    for obj in objs:
        if type(obj) == Polygon:
            for v in obj.vertices:
                yield v
        elif type(obj) == Points:
            for v in obj.vectors:
                yield v
        else:
            raise TypeError("Unrecognized obj : {}".format(obj))

def draw(*objs, origin=True, axes=True,  nice_aspect_ratio=True, width=6):
    all_vecs = list(extract_vectors(objs))
    xs, ys = zip(*all_vecs)
    max_x, max_y, min_x, min_y = max(0, *xs), max(0, *ys), min(0, *xs), min(0, *ys)

    if origin:
        plt.scatter([0], [0], color='k', marker='x')
    if axes:
        plt.gca().axhline(linewidth=2, color='k')
        plt.gca().axvline(linewidth=2, color='k')
    
    for obj in objs:
        if type(obj) == Polygon:
            for i in range(0, len(obj.vertices)):
                x1, y1 = obj.vertices[i]
                x2, y2 = obj.vertices[(i+1)%len(obj.vertices)]
                plt.plot([x1, x2], [y1, y2], color=obj.color)
            if obj.fill:
                xs = [v[0] for v in obj.vertices]
                ys = [v[1] for v in obj.vertices]
                plt.gca().fill(xs, ys, obj.fill, alpha=obj.alpha)
        elif type(obj) == Points:
            xs = [v[0] for v in obj.vectors]
            ys = [v[1] for v in obj.vectors]
            plt.scatter(xs, ys, color=obj.color)
        else:
            raise TypeError("Unrecognized obj : {}".format(obj))
    
    fig = matplotlib.pyplot.gcf()

    if nice_aspect_ratio:
        coords_height = (ylim()[1] - ylim()[0])
        coords_width = (xlim()[1] - xlim()[0])
        fig.set_size_inches(width, width * coords_height / coords_width)
    plt.show()


dino_vecs = [
    (6,  4), (3, 1), (1, 2), (-1, 5), (-2, 5), (-3, 4), (-4, 4),
    (-5, 3), (-5, 2), (-2, 2), (-5, 1), (-4, 0), (-2, 1), (-1, 0), (0, -3),
    (-1, -4), (1, -4), (2, -3), (1, -2), (3, -1), (5, 1)
]

print(dino_vecs)
print(*dino_vecs)
draw(
    Polygon(*dino_vecs, color="red")
)

간소화 한다고 드로잉 요소로는 아까 많이쓴 포인트랑 폴리곤만 넣었고

드로우 함수에서 plt 그리드 파트가 너무 복잡하게 생겨 그냥 뺏음

그렇게 기존의 드로잉 모듈 전체 코드는 140줄 정도 되었는데

내가 수정한 부분은 대충 70줄 정도

맨 아래 공룡 드로잉 포함해서 80줄 정도로 된다.

이렇게 70줄 정도 만으로

2차원에서 폴리곤과 포인트를 간편하게 드로잉하는 모듈 만들었다.

그리드와 틱은 좀 아쉽긴한데

이 정도 구조면 나도 나중에 필요할때 만들어서 써먹을수 있을듯싶네