Canvasの準備

前回までは左下原点の座標系でした。今回は中心原点の一般的なデカルト座標系を使います。また、円運動になるので長方形より正方形の方が良いでしょう。

<canvas id="canvas" width="640" height="640"></canvas>

const canvas = document.getElementById('canvas'),
      ctx = canvas.getContext('2d');

const dpr = window.devicePixelRatio || 1, // デバイスピクセルレシオ
      width = canvas.width,   // Canvas横幅
      height = canvas.height; // Canvas縦幅

/** 高解像度対応 */
// Canvasをピクセル比で拡大
canvas.width *= dpr;
canvas.height *= dpr;
// CSSで元のサイズに戻す
canvas.style.width = width + 'px';
canvas.style.height = height + 'px';
// Canvasの描画自体を拡大
ctx.scale(dpr, dpr);

/** 座標系の設定 */
// y座標を反転
ctx.scale(1, -1);
// x軸, y軸に沿って原点が中心にくるようにずらす
ctx.translate(width / 2, -height / 2);

scale と translate を合わせて一行で書くとこうです。

ctx.transform(dpr, 0, 0, -dpr, width * dpr / 2, height * dpr / 2);

考え方

前回までと同じく、時間を刻んで、位置の変化と速度変化を計算していくのですが、今回は加速度（重力）が位置によって変化します。

二つの物体間に働く重力（万有引力）は、次の式で計算できます。
（F は万有引力の大きさ、G 万有引力定数、M, m はそれぞれの物体の質量、r は物体間の距離）

$$F=G\frac{Mm}{r^2}$$

これを運動方程式に入れてみると、

$$ma=G\frac{Mm}{r^2}$$
$$a=G\frac{M}{r^2}$$

m が消せます。なので、重力のみを受けて運動する物体の加速度は、物体自身の質量は関係なく、重力源の質量のみに依存することがわかります（だから重いものも軽いものも落下速度は同じになります）。

つまり、ここでは人工衛星の質量は十分小さく、重力を及ぼさないと考えるので、質量も考えなくてよいことになります。必要なのは地球の質量と、地球中心からの距離だけです。また、地球中心は座標原点[0, 0]に固定とします。

単位と定数

単位として距離はkm、質量はkgを用いることとします。

万有引力定数の値として 6.674 × 10^-20km³s^-2kg^-1 を用います。他に地球の質量は 5.972 × 10²⁴kg、半径は 6.378 × 10³km とします。これらの定数は予め定義しておきましょう。

const G  = 6.674 * 1e-20, // 万有引力定数
      Me = 5.972 * 1e24,  // 地球質量
      Re = 6.378 * 1e3;   // 地球半径

今回は地球重力のみを考えるので、地球質量と万有引力定数をかけた地球重力係数も定数として用意します。それから1pxが何kmとなるかを決めるスケール定数も準備しておきましょう。

// 地球の重力係数
const K = G * Me;

// スケール
const scale = 30; // [km/px]

天体クラス

天体を扱うためのクラスを作ります。位置は内部的にはXY座標で扱いますが、初期化時は極座標（角度と距離）で指定できるようにします。

三角関数 Math.cos, Math.sin はラジアンを引数に取るので、直感的な度数法で計算できる補助関数を作っておくと便利です。

/** sin */
const sin = deg => {
  deg %= 360;

  if (deg === 180) {
    return 0; // Math.sin では誤差が出るので正確な値を返す
  }

  return Math.sin(deg / 180 * Math.PI);

};
/** cos */
const cos = deg => {
  deg %= 360;

  if (deg === 90 || deg === 270) {
    return 0; // Math.cos では誤差が出るので正確な値を返す
  }

  return Math.cos(deg / 180 * Math.PI);
};

中心からの距離を r 、x軸と成す角を deg とします。

let objectList = []; // 管理用配列

/** 天体クラス */
class AstroObject {
  constructor(attr) {
    // 位置
    this.x = attr.r * cos(attr.deg);
    this.y = attr.r * sin(attr.deg);
  }
}

速度も同様に、大きさと角度で指定できるようにします。角度の基準はどこでもいいのですが、ここではr軸（半径方向）を基準にします。極座標系の回転方向が常に90度となる向きです。この角度に位置の角度を足したものが、速度のx軸に対する角度となります。

速度の大きさをv、r軸に対する角度をvdegとします。

/** 天体クラス */
class AstroObject {
  constructor(attr) {
    ...
    // 速度
    this.vx = attr.v * cos(attr.vdeg + attr.deg);
    this.vy = attr.v * sin(attr.vdeg + attr.deg);

    // 配列に格納
    objectList.push(this);
  }
}

前述の通り質量は不要です。複数の衛星を扱えるように、コンストラクター内で配列に自分自身を格納しています。

中心からの距離rと角度rad, degを逆算するゲッターも設定しておきましょう。

...
// 中心からの距離
get r() {
  return Math.hypot(this.x, this.y);
}
// 角度
get rad() {
  // x = 0
  if (!this.x) {
    if (!this.y) {
      return 0; // [0, 0] の場合は0を返す
    }
    if (y < 0) {
      return Math.PI * 2.5; // 270°
    } else {
      return Math.PI * .5;  // 90°
    }
  }
  let rad = Math.atan(this.y / this.x); // arctan(y/x) : -π/2 〜 π/2
  if (this.x < 0) {
    // 第二・第三象限
    rad += Math.PI;
  } else if (this.x > 0 && this.y < 0) {
    // 第四象限 : 負数になってるので正数にする
    rad += Math.PI * 2;
  }
  return rad;
}
get deg() {
  return this.rad * 180 / Math.PI;
}
...

描画メソッドと更新メソッドを追加します。

...
draw() {
  // 4px × 4px の正方形で描画する
  ctx.beginPath();
  ctx.rect(this.x / scale - 2, this.y / scale - 2, 4, 4);
  ctx.fill();
}
update(dt) {
  const gravity = -K / this.r ** 2; // 現在位置の重力加速度

  this.x  += this.vx * dt;
  this.y  += this.vy * dt;
  this.vx += gravity * cos(this.rad) * dt;
  this.vy += gravity * sin(this.rad) * dt;
}
...

組み立て

Canvasへの描画や更新は前回までと基本的に同じです。ちょっとだけ付け足して、地球を描画する処理を組み込みます。

const fps = 30,     // FPS
      dt = 1 / fps; // 1フレーム当たりの秒数
  
let startTime,      // 開始時間
    frameCount = 0; // 経過フレーム数

/** 地球描画 */
const drawEarth = () => {
  ctx.save();

  ctx.fillStyle = '#25A8FF'; // 塗り色を青に
  ctx.beginPath();
  ctx.arc(0, 0, Re / scale, 0, Math.PI * 2); // 中心に円を描く
  ctx.fill();

  ctx.restore();
}

/** tick */
const tick = (time) => {
  if (!startTime) startTime = time; // 初回呼び出し時に開始時間を設定
  const currentFrame = Math.floor((time - startTime) * fps / 1000); // 現在の経過フレーム数
  let df = currentFrame - frameCount // 前回からのフレーム増分

  // フレーム増分だけ更新処理を繰り返す
  while (df) {
    // オブジェクトの更新
    objectList.forEach((obj, i) => {
      if (obj) {
        obj.update(dt);

        // 衝突判定
        if (obj.r <= Re) {
          objectList[i] = null; // 中心星に衝突した物体は取り除く
        }
      }
    });
    df--;
  }

  objectList = objectList.filter(Boolean); // 衝突した物体を配列から除去

  // 画面の消去
  ctx.clearRect(-width / 2, -height / 2, width, height);

  // 再描画
  drawEarth();
  objectList.forEach(obj => {
    obj.draw();
  });

  frameCount = currentFrame; // フレーム数を更新
  requestAnimationFrame(tick);
};

これで準備ができました。

実践

実際に人工衛星を飛ばして見ましょう。ここではISSと同じ高度400kmの円軌道に乗せることとします。

円軌道の軌道速度は以下の式で計算できます。

$$v=\sqrt{\frac{mG}{r}}$$

計算すると 約7.67km/s です。この値を使って天体クラス AstroObject のインスタンスを作ります。

const satellite = new AstroObject({
     r: Re + 400, // 地球半径 + 高度
   deg: 90,
     v: 7.67,
  vdeg: 90
});

あとは tick() を起動すればシミュレート開始です。

See
the Pen
Orbital Motion by SHIN Inc. (@shin-inc)
on CodePen.

どうでしょうか？　円と小さな四角が表示されてるだけに見えるかもしれません。30kmを1pxで表示しているので、秒速7.67kmだと1px移動するのに約4秒かかることになります。スケールに対して速度が遅いので、ほとんど動いていないように見えるのです。高度400kmでは一周するのに約一時間半かかります。

このまま一時間半だらだら見続けるのは大変なので、倍速で動かせるように改良します。やり方は単純に1フレーム内での更新処理を倍速度分繰り返すだけです。

const speed = 200; // 再生速度
...
/** tick */
const tick = (time) => {
  ...
  let df = currentFrame - frameCount // 前回からのフレーム増分

  df *= speed; // 速度分処理回数を増やす

  // フレーム増分だけ更新処理を繰り返す
  while (df) {
    ...
  }
  ...
}

これで高速で再生されるようになります。

See
the Pen
Orbital Motion x200 by SHIN Inc. (@shin-inc)
on CodePen.

円軌道を描いて地球の周りを回るのがよくわかると思います。速度や高度を変えて色々と試してみてください。

複数個の衛星を同時に飛ばすこともできますが、当然数が増えるほど負荷は増えます。また、速度も上げるほど負荷が増えて処理落ちするので注意してください。

課題

今回は地球周回軌道のシミュレーションを作成しました。最後に今後の改良点をいくつか挙げておきます。

中心星の変更

地球固定ではなく、中心星も自由にパラメーターをいじれるようにオブジェクト化すると応用幅がぐっと広がります。火星や金星はもちろん、架空の惑星の周りの軌道を自由にシミュレートできるようになります。

複数重力源の対応

人工衛星のような小さい天体ではなく、地球と月のように複数の重力源を扱えるようにしたいです。地球周回軌道程度なら月の重力は無視できますが、月に近い軌道を扱いたい場合は両方の重力を考えなければなりません。

複数の重力源を扱う場合、重力による加速度計算と位置の更新は別々に扱う必要があり、処理が複雑になります。

数値計算の改良

現在の積分計算はオイラー法と呼ばれるもので、これは精度があまりよくありません。そこでより精度の高いリープ・フロッグ法などに置き換えることが考えられます。厳密な計算でなければオイラー法でも問題ないので、優先度は低いです。

おまけ

地上での直感に反するような、軌道上における物体の振る舞いを紹介します。

前を行く赤い物体は高度400kmを速度7.67km/sで飛行しています。後ろから青い物体が同じ高度で速度8km/sで追いかける様子です。

See
the Pen
Orbital Chaser Motion by SHIN Inc. (@shin-inc)
on CodePen.

ご覧の通り、最初こそ距離が縮まるように思えますが、次第に距離が離れていき、結局青い物体は赤い物体に追いつけないのがわかります。軌道上で後ろから追いかけようとしたら、単純に速度を上げるだけではダメなのです。

ではどうしたらいいのか……考えてみてください。