2017-05-03프로그래밍 / 그래픽스

[JavaScript로 천체 구현하기] 행성의 움직임을 구현해보자

이번 포스팅에서는 저번 포스팅에 이어 실제 궤도의 모양과 크기, 위치, 방향을 정의하고 JavaScript코드로 작성을 해보려고 한다. 실제 어플리케이션을 작성할 때는 TypeScript를 사용하였으나, 편의상 JavsScript ES6로 포스팅을 진행한다.

궤도를 구하는 방법은 저번 포스팅의 용어 정리에서 언급했던 케플러 6요소를 이용하면 된다.

데이터 정의하기

지구의 궤도 데이터는 다음과 같다.

const AU = 149597870;
const EARTH_ORBIT = {
  base: {
    a: 1.00000261 * AU,
    e: 0.01671123,
    i: -0.00001531,
    o: 0.0,
    l: 100.46457166,
    lp: 102.93768193
  },
  cy: {
    a: 0.00000562 * AU,
    e: -0.00004392,
    i: -0.01294668,
    o: 0.0,
    l: 35999.37244981,
    lp: 0.32327364
  }
};

먼저 base프로퍼티에 들어있는 값은 천문학에서의 역기점인 2000년 1월 1일 정오, 즉 J2000때 측정된 궤도의 케플러 6요소들을 의미한다.

그리고 cy프로퍼티에 있는 값들은 1세기당 궤도 요소들의 변화량이다. a는 장반경, e는 이심률, i는 기울기, o는 승교점 적경, l은 평균 경도, lp는 근일점 경도를 의미한다.

이 중 장반경의 경우 태양과 지구 사이의 거리를 의미하는 AU 단위로 선언되어있기 때문에 1AU를 km로 환산한 값인 149597870를 곱해줘야 한다.

문제는 이 데이터에는 케플러 6요소 중 근일점 편각과 근일점 통과시각이 없다는 점이다.

먼저 근일점 편각은 이미 주어진 근일점 경도와 승교점 적경을 사용하여 계산할 수 있으며, 근일점 통과 시각은 어차피 궤도 계산에 필요한 평균근점이각(Mean anomaly)을 구할 때 필요한 것인데 평균근점이각을 구할 때는 근일점 통과 시각보다 근일점 경도와 평균 경도를 사용하는 방법이 훨씬 간단하기 때문에 오히려 편해진 상황이다.

시간 설정과 궤도 요소 계산하기

참고로 천문학에서는 우리가 평소 사용하는 그레고리력이 아닌 율리우스력을 사용하는데, 이건 기원전 4713년 1월 1일 월요일 정오를 기점으로 계산하는 것이기 때문에 현재 시점인 서기 2017년까지의 날짜를 다루기에는 자릿수가 너무 커진다. 그래서 이번 작업에는 1976년 IAU(국제천문연맹)에서 결정한 J2000을 사용하는 것이다.

이제 날짜 요소들을 선언해보자.

const DAY = 60 * 60 *24;
const YEAR = 365.25;
const CENTURY = 100 * YEAR;
const J2000 = new Date('2000-01-01T12:00:00-00:00');

let today = new Date();
let epochTime = (today - J2000) / 1000;
// Date객체끼리 연산하면 값이 ms로 나오기 때문에 1000으로 나눠서 초 단위로 나오게끔 환산해준다
let tDays = epochTime / DAY; // 일수로 환산
let T = tDAys / CENTURY; // 몇 세기나 지났는지 환산

이 글을 작성하고 있는 지금은 2017년 5월 3일 22시 27분이다. 계산 결과 J2000으로부터 약 6332.06028991일이 지났으며 0.17336236세기가 지났다는 것을 알 수 있었다.

그럼 이제 아까 선언한 데이터 중 지구의 한 세기당 변화량에 방금 구한 0.17336236을 곱한 후 궤도의 기본 요소 값에 더해주면 J2000이후 오늘 이 시점의 궤도 요소 값이 나올 것이다.

const keys = Object.keys(EARTH.base);

let computed = { time: epochTime };

computed = keys.reduce((carry, el) => {
  const variation = EARTH.cy || 0;
  carry[el] = EARTH.base[el] + (variation * T);

  return carry;
}, computed);
/*
  a: 149598406.2031184
  e: 0.016703615925039474
  i: -0.0022597770311920135
  l: 6341.400827688619
  lp: 102.9937254119609
  o: 0.0
*/

부궤도 요소 계산하기

이제 현재 시점 기준의 궤도 요소 값을 구했으면 부궤도 요소를 구할 차례이다.

가장 먼저 케플러 6요소 중 하나지만 우리의 데이터에는 없는 근일점 편각을 구해보도록 하자. 근일점 편각 w는 근일점 경도에서 승교점 적경을 뺀 값이다. 하지만 지구의 승교점 적경 o는 0.0이므로 근일점 경도와 근일점 편각이 일치하게 된다.

computed.w = computed.lp - computed.o;
// w = 102.9937254119609 - 0.0 = 102.9937254119609

다음은 평균근점이각을 구해야한다. 평균근점이각은 어떤 물체가 공전 속도와 공전 주기를 유지한 채 정확한 원 궤도로 옮겨간다고 가정했을 때 물체와 궤도 근점간의 각거리를 의미한다.

즉 위 그림에서 실제 물체의 궤도는 회색 궤도이나, 저 궤도를 원이라고 가정한 빨간 궤도에서의 각거리를 구하는 것이다.

저 두 궤도의 공전 주기는 같기 때문에 같은 지역을 같은 시간 동안 지나간다. 근데 그림을 보면 원 궤도는 각 호들의 넓이가 일정하지만 타원 궤도인 회색 궤도는 일정하지 않다.

타원 궤도를 도는 물체는 초점과의 거리에 따라 각속도가 변하기 때문이다. 그래서 좀 더 구하기 쉬운 평균근점이각을 먼저 구한 후 타원 궤도에 이를 적용하게 된다. 이러한 접근 방법이 가능한 이유는 케플러 제2법칙인 면적속도 일정의 법칙 때문이다.

평균근점이각을 구하는 공식은 여러 개가 있으나 필자는 이미 주어진 근일점 경도와 평균 경도값을 가지고 평균근점이각을 구하는 공식을 사용하려고 한다.

평균근점이각 $M$ 은 평균 경도 $l$ 과 근일점 경도 $\omega$ 의 차로 나타내어 질 수 있다.

M = l-\omega

computed.M = computed.l - computed.lp
/* M = 6341.42938740066 - 102.99372566842652 =  6238.435661732234 */

평균근점이각까지 구했다면 이제는 편심이각(Eccentric anomaly)을 구할 차례이다. 이때 앞으로의 계산을 원활하게 하기 위해 Digree 단위로 표현되어있는 값들을 모두 Radian으로 변경해주겠다.

const DEG_TO_RAD = Math.PI / 180;

computed.a *= 1000; // 149598406204.91275
computed.i *= DEG_TO_RAD; // -0.000039441031811489665
computed.o *= DEG_TO_RAD; // 0.0
computed.w *= DEG_TO_RAD; // 1.7975796293754245
computed.M *= DEG_TO_RAD; // 108.8812424710587

편심이각은 타원 카플러 궤도를 따라 움직이는 물체의 위치를 결정하는 중요한 궤도 요소이다. 그림으로 보면 더 직관적으로 이해하기가 쉽다.

위 그림에서 P는 궤도를 돌고 있는 행성위 위치를 나타내며, 그에 따른 편심이각을 E로 나타내고 있다. 그리고 타원의 중심은 C이고 타원의 초점은 F이다.

이때 편심이각 E는 타원의 중심에 꼭지점 하나를 찍고 궤도 장반경과 같은 길이의 빗변 a를 그은 후, 장반경 e와 수직하면서도 P에 닿도록 선분을 그어 만들어진 직각삼각형에서 관찰되는 각이다.

평균근점이각과 마찬가지로 편심이각도 여러 개의 정의로 나타내어질 수 있는데, 필자는 이미 평균근점이각을 구했기 때문에 평균근점이각으로부터 유도되는 공식을 사용하겠다.

E = \frac{M+e{sinM}}{1-e{cosM}}

여기서 e는 이심률, M은 평균근점이각을 의미한다. 그러나 한가지 슬픈 사실은 여기서 나오는 E값이 근사값이라는 것이다. 그렇기 때문에 보통 이 공식을 통해 나온 해는 Newton-Raphson method를 사용하여 오차를 줄여줘야한다.

Newton-Raphson method는 미분 가능한 함수의 해를 찾는 방법인데, 계속적인 급수 형태로 계산되어지며 계산을 반복할수록 더 정확한 근사치를 뱉어준다고 보면 된다.

값의 오차는 이심률이 높을 수록 더 커지게 되는데, 데이터 상 지구의 이심률은 0.0167703정도니까 오차율이 그렇게 크진 않을 것이라고 예상된다. 그럼 코드를 작성해보자.

function getEccentricity(callback, x0, maxCount) {
  let x = 0;
  let x2 = x0;
  for(let i = 0; i < maxCount; i++) {
    x = x2;
    x2 = callback(x);
  }
}

function kepler(e, M) {
  return x => {
    return x + (M + e * Math.sin(x) - x) / (1 - e * Math.cos(x));
  };
}

computed.E = getEccentricity(kepler(computed.e, computed.M), computed.M, 6);
/* E =108.8962365500302 */

행성의 위치 도출하기

이렇게 특정 날짜를 기준으로 행성의 위치를 알아내기 위해 필요한 값들을 모두 구했다.

이 값들이 있다면 근일점 쪽을 X축으로 하는 황도좌표평면에 대한 직각 좌표 값을 계산할 수 있으며, 이 좌표 값의 유클리드거리 r과 항성과 궤도의 근일점 기준으로 어느 각도에 행성이 위치하고 있는지를 나타내는 진근점이각 v를 통해 행성의 위치를 표현할 수 있다.

즉, 행성은 항성으로부터 r만큼의 거리만큼 떨어져 있고 궤도의 근일점으로부터 v만큼 돌아간 위치에 존재하고 있다는 것이다. 이렇게 거리랑 각만 알면 3D Scene내에서의 로컬좌표를 정의할 수 있게된다.

const RAD_TO_DEG = 180 / Math.PI;

computed.pos = {x: null, y: null};
computed.pos.x = computed.a * (Math.cos(computed.E) - computed.e);
computed.pos.y = computed.a * (Math.sqrt(1 - (computed.e * computed.e))) * Math.sin(computed.E);

computed.r = Math.sqrt(Math.pow(computed.pos.x, 2) + Math.pow(computed.pos.y, 2));
computed.v = Math.atan2(computed.pos.y, computed.pos.x);

computed.r /= (1000 * AU);
computed.v *= RAD_TO_DEG;

/*
  r = 1.0081888306835929 AU
  v = 120.17208256525967 도
*/

이렇게 이 포스팅 작성 시점을 기준으로 지구는 태양으로부터 1.0081...AU 떨어진 위치에 지구 궤도의 근일점으로부터 약 120도 돌아간 곳에 위치하고 있다는 결과를 얻을 수 있다.

마치며

사실 제일 마지막 pos를 구하는 공식과 진근점이각을 구하는 공식은 아직 이해가 잘 안된다. 하지만 중요한 것은 이런 공식을 사용해서 내가 원하는 뭔가를 만들어 볼 수 있다는 것이 아닐까 생각한다.

이 공식들을 사용하여 필자가 만든 시뮬레이터는 여기에서 확인해볼 수 있고, 전체 소스는 Solarsystem 프로젝트 깃허브 레파지토리에서 확인해볼 수 있다.

이상으로 행성 위치 게산하기 포스팅을 마친다.