구면삼각법과 적도좌표계

이번에는 구면삼각법을 공부하여, 실제 적도좌표계와 지평좌표계에서 적용해 보고 이해를 확장해 보려고 한다.

생각은 간결하지만 이것이 간단한 주제도 아니고 연관된 범위도 아주 넓다.

그래서 엄밀한 수식 증명은 생략하되, 언제나 그렇듯이 그래픽 구현 검증으로 직관을 높일 수 있는 것이 목적이다.

 

그럼, 먼저 기본적인 구면삼각법 정리부터 하자!

sphere 위에 구면삼각형

대원 : 구를 그 중심(o)을 지나는 평면으로 자를 때에 생기는 원 또는 그 둘레 A, B, C : 점 A에서 교차하는 두 대원이 이루는 이면각 a, b, c : 점 A, B, C의 대응호(BC, AC, AB)가 이루는 중심각     ex) b = ∠AoC 구면초과 E [rad] = A + B + C - π ⇒ 입체각 steradian[S/R2] 구면삼각형 면적 S = E * R2   ① sin(A)/sin(a) = sin(B)/sin(b) = sin(C)/sin(c) ② cos(a) = cos(b)*cos(c) + sin(b)*sin(c)*cos(A)

위 수식에 대한 자세한 증명은 https://blog.naver.com/sluggeryck/220577238348 참고 바랍니다!

 

 

이어서 각 좌표계의 정의를 살펴보고 이해를 돕기위해 내 주석을 추가해보겠다!

구면좌표계(spherical coordinate system) - (r, θ, ϕ) r은 0부터 ∞까지, 양의 방향의 z축과 이루는 각도 θ는 0부터 π까지, z축을 축으로 양의 방향의 x축과 이루는 각 ϕ는 0부터 2π까지의 값을 갖는다 ⇒ 전자기학 등에서 일반적으로 많이 쓰는데 기본이되는 좌표계로 써보면 편하다!

적도좌표계(equatorial coordinate system) 천구의 적도와 황도면과의 교점이 되는 춘분점을 경도의 원점으로 잡는다. 서쪽에서 동쪽으로 측정한 경도를 적경(RA, right ascension) 이라 하고 천구의 적도면에서 북쪽을 +로, 남쪽을 -로 하여 0~90°까지의 각으로 나타낸 것을 적위(DEC, declination)라 한다. ⇒ 위 구에서 원점o는 지구, 노란껍질은 천구면으로 상상하자. 파란색 z축은 지구 자전축(DEC=90°)이고 빨강색 x축은 춘분선(RA=0h)이다. z축에 있는 북극성(B)에서 지구를 바라봤을 때, 지구는 반시계방향으로 자전과 공전을하지만 별들은 천구에 붙박이다. 춘분일때 태양은 빨강색 x축 춘분선에 있다. 예시로 북극성(B, DEC=89°16', RA=2h31m), 직녀성(C, DEC=38°47', RA=18h36m), 아크투르스(A, DEC=19°11', RA=14h15m)를 천구면에 표시하였다. 시간과 장소에 따라 좌표값이 변하지 않기 때문에 천체의 위치를 나타내는 데 가장 널리 쓰인다.

지평좌표계(horizontal coordinate system) 지평좌표계는 관측자의 지평면(horizon)을 기준면으로 사용하는 좌표계이다. 고도(altitude 또는 elevation, 0° ~ 90°)와 방위각(azimuth, 0° ~ 360°)을 사용하여 천체의 위치를 정의한다. 관측자 중심으로 천체의 위치를 나타내므로 특정한 시각과 장소에서 천체의 위치를 쉽게 찾을 수 있다 ⇒ 마찬가지로 위 구에서 원점o는 지구, 노란껍질은 천구면으로 상상하자. 지구 위 관찰자가 천정에 A별이 있을 때, C별을 바라본 다면 고도와 방위각이 어떻게 될까? 아마 고도는 90°-(b=59°)=31°, 방위각은 진북 A=55.7° 일 것이다.

 

구면삼각법으로 별의 좌표계 계산까지 가능한 것을 보았으니, 이제 아래 시간개념을 더 하여 응용한다면 간단한 별자리판이나 허접하게나마 스텔라리움을 흉내낼 수 있을 것이다. 물론 프로그래밍이 쉬운 것은 아니지만...

광막한 우주와 별들이 단순히 신비롭기만 한 것에서 이제 경이로움을 느끼게 되고 더 친밀해 지는 순간이다. 한편 이제 우주에서 내가 밟고 있는 이 세상이 얼마나 초라하게 변방인지 첫 눈을 떴다고 할까?!

 

협정세계시(Coordinated Universal time, UTC)는 초는 원자시계에 따른 초를 그대로 쓰고, 시각은 평균태양시(태양이 남중하는 때가 12시 정오)로부터 멀어지지 않도록 조정하는 표준시이다. 한국 표준시(KST, Korea Standard Time)는 그리니치천문대를 지나는 본초자오선(prime meridian)으로부터 일본 표준시와 같은 동경 135도를 기준으로 하여 UTC보다 9시간 빠른 표준시(UTC+09:00)이다.

 

마지막으로 컴퓨터로 구면삼각법을 계산하고 좌표계에서 표시하여 검증하는 것으로 이번 주제는 끝내겠다.

 

from numpy import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import matplotlib.animation as animation

def ang2rad(x) : return x*pi/180
def rad2ang(x) : return x*180/pi
def h2rad(x) : return x*pi/12
def rad2h(x) : return x*12/pi
def s2r(r, th, ph) : return r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)
def r2s(x, y, z) : return sqrt(x**2+ y**2+ z**2), arctan2(sqrt(x**2+ y**2), z), arctan2(y, x)

def ina(v0, v1) : # included angle
    return arccos(dot(v0, v1)/sqrt(dot(v0, v0)*dot(v1, v1)))
def rot_mat(axis, theta) : # rotation matrix
    a = cos(theta / 2.0)
    b, c, d = -array(axis) * sin(theta / 2.0)
    aa, bb, cc, dd = a * a, b * b, c * c, d * d
    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
    return array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                  [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                  [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
def arc(v0, v1) :
    ang = ina(v0, v1)
    n = cross(v0, v1)
    n /= sqrt(dot(n, n))
    ps = array([])
    for th in linspace(0, ang, 16) :
        p = dot(rot_mat(n, th), v0)
        ps = append(ps, p)
    return reshape(ps, (int(len(ps)/3),3))
def st(a, b, c) :    # spherical trigonometry
    A = arccos((cos(a)-cos(b)*cos(c))/(sin(b)*sin(c)))
    B = arccos((cos(b)-cos(a)*cos(c))/(sin(a)*sin(c)))
    C = arccos((cos(c)-cos(a)*cos(b))/(sin(a)*sin(b)))
    return A,B,C

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal')

u = linspace(0, 2 * pi, 64)
v = linspace(0, pi, 64)

x = outer(cos(u), sin(v))     # outer(ones((5,)), linspace(-2, 2, 5))
y = outer(sin(u), sin(v))
z = outer(ones(size(u)), cos(v))     # outer(arange(0,5), arange(0,5))

ax.plot_surface(x, y, z, rstride =4, cstride =4, color ='y', alpha=0.6)

right_ascension = ax.plot([0, 2],[0, 0],[0, 0], 'r', linewidth=1) # x - axis
rotation_axis = ax.plot([0, 0],[0, 0],[0, 2], 'b', linewidth=1) # z - axis

# (r, th, ph) <- (r, pi/2-ang2rad(Dec), h2rad(RA))
v0 = array(s2r(1, pi/2-ang2rad(89+16/60), h2rad(2+31/60)))     # Polaris
v1 = array(s2r(1, pi/2-ang2rad(38+47/60), h2rad(18+36/60)))     # Vega
v2 = array(s2r(1, pi/2-ang2rad(19+11/60), h2rad(14+15/60)))     # Arcturus

a = arc(v0, v1)
b = arc(v1, v2)
c = arc(v2, v0)
ax.plot(a[:,0], a[:,1], a[:,2], 'g')
ax.plot(b[:,0], b[:,1], b[:,2], 'g')
ax.plot(c[:,0], c[:,1], c[:,2], 'g')
ax.plot([v0[0]], [v0[1]], [v0[2]], 'k.')
ax.plot([v1[0]], [v1[1]], [v1[2]], 'k.')
ax.plot([v2[0]], [v2[1]], [v2[2]], 'k.')

print('length :', 'a:',ina(v0,v1), 'b:', ina(v1,v2), 'c:',ina(v2,v0))
A, B, C = st(ina(v0,v1), ina(v1,v2), ina(v2,v0))
print('area   :', A+B+C-pi/2)
ax.text(v0[0], v0[1], v0[2], 'B '+str(round(rad2ang(B),1))+"˚")
ax.text(v1[0], v1[1], v1[2], 'C '+str(round(rad2ang(C),1))+"˚")
ax.text(v2[0], v2[1], v2[2], 'A '+str(round(rad2ang(A),1))+"˚")
ax.text(0, 0, 0, 'o')

ax.view_init(elev=50, azim=260)              # angle

plt.show(block=False)

 

이것도 막상 정리하려니 시간도 많이 걸리고 고생스럽네ㅎㅎㅎ

그래서 나 자신에게 토닥토닥^^;