3Dプログラミング入門講座・その11:物理演算その3・ケプラー方程式で惑星軌道シミュレーターを作る







太陽系惑星軌道計算シミュレーターフォーム by javascript.htm

太陽系惑星軌道計算シミュレーターフォーム by javascript.zip



x3dom、自作ライブラリを用いた太陽系シミュレーターのjavascriptのサンプル

初期化
solarsystem.js

...省略...

function init(b) {
  bBBB = b;
  if(0 <= b) Speed = Number(F1.SPD.value);
  else Speed = Number(F1.SPD.value) * -1;
  Zoom = Number(F1.ZOM.value);
  if(Zoom < 0.0) Zoom *= -1.0;

  var radioList = document.getElementsByName("CALC");
  for(i = 0; i < radioList.length; i++){
      if(radioList[i].checked){
          CalcWay = Number(radioList[i].value);
          break;
      }
  }

  bWaiting = false;
  if(!CalcWay) InitKepler();
  else InitRelative();
  onWaiting();
}

function onWaiting() {
  if(!bWaiting) {
    fFst = -1;
    setmp();
    setline();
    if(bBBB) bRunning = true;
  }
}

function onRunning() {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;

  if(fFst != 0) {
    fFst = 0;
    time = 0.0;
    if(CalcWay) dt = Speed * 60 * 60;
    else dt = Speed * msecPerDay;
  }
 
  //メインループ
  if(CalcWay) UpdateFrameRelative();
  else UpdateFrameKepler();
}

function getnow() {
  var nt = new Date(Number(F1.T1.value), Number(F1.T2.value) - 1, Number(F1.T3.value), Number(F1.T4.value), Number(F1.T5.value), Number(F1.T6.value));
  return nt;
}

function InitKepler() {
  if(!bRead) fFst = 1;
  dat = getnow();
  PlanetInit(dat);
}


初期化実装
solarsystem.js

...省略...

//惑星初期化
function PlanetInit(dat) {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var i;
  var j;
  if(0 < fFst) {
    //データファイル読み込み
    if(!bWaiting) {
      bWaiting = true;
      readCSV('./javascripts/nas6lib/PData000.txt', 'analyzeCSV', 'readedCSV');
    }
    return;
  }
  else {
    for(i = 0; i < planetnum; i++) {
      var dat0 = planet[i].m_dat0;
      var datt = dat.getTime();
      var dat0t = dat0.getTime();
      var ddat = (datt - dat0t) / msecPerDay;
      var nday = ddat;

      var xx = new Array(new N6LVector(3));
      var f = planet[i].kepler(nday, xx);
      planet[i].x0 = new N6LVector(3);
      planet[i].x0.x[0] = xx[0].x[0];
      planet[i].x0.x[1] = xx[0].x[1];
      planet[i].x0.x[2] = 0.0;

      var xyz = new Array(new N6LVector(3));
      planet[i].ecliptic(planet[i].x0.x[0], planet[i].x0.x[1], planet[i].x0.x[2], xyz);
      if(isNaN(xyz[0].x[0]) || isNaN(xyz[0].x[1]) || isNaN(xyz[0].x[2])) {
        planet[i].x0.x[0] = 0.0;
        planet[i].x0.x[1] = 0.0;
        planet[i].x0.x[2] = 0.0;
      }
      else {
        planet[i].x0.x[0] = xyz[0].x[0];
        planet[i].x0.x[1] = xyz[0].x[1];
        planet[i].x0.x[2] = xyz[0].x[2];
      }

      planet[i].v0 = new N6LVector(3);
      
      //ケプラー方程式から軌道速度を求める
      var xyz2 = new Array(new N6LVector(3));

      var xxx = new Array(new N6LVector(3));
      planet[i].kepler(nday + (1.0 / (24.0 * 4.0) * planet[i].m_t), xxx);
      var vv = xxx[0].Sub(xx[0]);
      //速度微調整
      planet[i].v0.x[0] = (vv.x[0] / (60.0 * 60.0 * 24.0 / (24.0 * 4.0) * planet[i].m_t) / planet[i].CNST_C) * planet[i].m_mv;
      planet[i].v0.x[1] = (vv.x[1] / (60.0 * 60.0 * 24.0 / (24.0 * 4.0) * planet[i].m_t) / planet[i].CNST_C) * planet[i].m_mv;
      planet[i].v0.x[2] = 0.0;

      planet[i].ecliptic(planet[i].v0.x[0], planet[i].v0.x[1], planet[i].v0.x[2], xyz2);
      if(isNaN(xyz2[0].x[0]) || isNaN(xyz2[0].x[1]) || isNaN(xyz2[0].x[2])) {
        planet[i].v0.x[0] = 0.0;
        planet[i].v0.x[1] = 0.0;
        planet[i].v0.x[2] = 0.0;
      }
      else {
        planet[i].v0.x[0] = xyz2[0].x[0];
        planet[i].v0.x[1] = xyz2[0].x[1];
        planet[i].v0.x[2] = xyz2[0].x[2];
      }
      mp[i] = new N6LMassPoint(planet[i].x0, planet[i].v0, planet[i].m_m, planet[i].m_r);
    }
  }
}

function readedCSV(res) {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  bWaiting = false;
  bRead = true;

  for(i = 0; i < planetnum; i++) {

    var PlanetName = res[i][0];
    var PlanetNo = Number(res[i][1]);
    var EpochYY = Number(res[i][2]);
    var EpochMM = Number(res[i][3]);
    var EpochDD = Number(res[i][4]);
    var Epochh = Number(res[i][5]);
    var Epochm = Number(res[i][6]);
    var Epochs = Number(res[i][7]);
    var a = Number(res[i][8]);
    var e = Number(res[i][9]);
    var m0 = Number(res[i][10]);
    var npd = Number(res[i][11]);
    var ra = Number(res[i][12]);
    var rb = Number(res[i][13]);
    var p = Number(res[i][14]);
    var ss = Number(res[i][15]);
    var ii = Number(res[i][16]);
    var ww = Number(res[i][17]);
    var m = Number(res[i][18]);
    var r = Number(res[i][19]);
    var mv = Number(res[i][20]);

    var dat0 = new Date(EpochYY, EpochMM - 1, EpochDD, Epochh, Epochm, Epochs);
    var datt = dat.getTime();
    var dat0t = dat0.getTime();

    var ddat = (datt - dat0t) / msecPerDay;

    //惑星初期化
    planet[i] = new N6LPlanet();
    planet[i].Create(PlanetNo, PlanetName, ddat, dat0, a, e, m0, npd, ra, rb, p, ss, ii, ww, m, r, mv);

    //質点初期化
    mp[i] = new N6LMassPoint(planet[i].x0, planet[i].v0, m, r);
  }
  return true;
}

//惑星軌道線分設定
function setline() {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var i;
  var j;
  var k;
  var n = 32;
  var str;
  for(i = 1; i < planetnum; i++) {
    str = "";
    var x0;
    //惑星1周を32分割の線分設定
    for(j = 0; j < n; j++) {
      var ad = (360.0 * 360.0 / 365.2425 / planet[i].m_nperday) * (j / n);
      var days = (dat.getTime() - planet[i].m_dat0.getTime()) / msecPerDay;
      var nday = days + ad;
      var xx = new Array(new N6LVector(3));
      var f = planet[i].kepler(nday, xx);
      var x1 = new N6LVector(3);
      x1.x[0] = xx[0].x[0];
      x1.x[1] = xx[0].x[1];
      x1.x[2] = 0.0;
      if(j == 0) x0 = new N6LVector(x1);
      str += (x1.x[1] / CNST_AU / Zoom).toString() + " " + (-x1.x[0] / CNST_AU / Zoom).toString() + ", ";
    }
    str += (x0.x[1] / CNST_AU / Zoom).toString() + " " + (-x0.x[0] / CNST_AU / Zoom).toString();

    var ss = planet[i].m_s * planet[i].CNST_DR;
    var ii = planet[i].m_i * planet[i].CNST_DR;
    var ww = planet[i].m_w * planet[i].CNST_DR;

    var vec = new N6LVector(3);
    var mat = new N6LMatrix(3);
    mat = mat.UnitMat().RotAxis(vec.UnitVec(2), ss).RotAxis(vec.UnitVec(1).Mul(-1.0), ii).RotAxis(vec.UnitVec(2), ww);
    var VecWK = new N6LVector(4);
    var MatWK = new N6LMatrix(4);
    MatWK.x[0] = VecWK.UnitVec(0);
    MatWK.x[0].bHomo = false;
    for(k = 1; k < 4; k++) {
      MatWK.x[k] = mat.x[k - 1].NormalVec().ToHomo();
      MatWK.x[k].x[0] = 0.0;
      MatWK.x[k].bHomo = false;
    }
    VecWK = MatWK.NormalMat().Vector();


    var elm;
    var sp;

    elm = document.getElementById(IDL[i]);
    elm.setAttribute('lineSegments', new String(str));

    elm = document.getElementById(IDT[i]);
    sp = VecWK.ToX3DOM();
    elm.setAttribute('rotation', sp.toString());
  }
}
 

CSVファイルから惑星諸定数を読み込み
planet,masspoint各クラスの構造体に保存します

CSVファイルから読み込んだ構造体から
planetのインスタンス
//惑星初期化
planet[i] = new N6LPlanet();
planet[i].Create(PlanetNo, PlanetName, ddat, dat0, a, e, m0, npd, ra, rb, p, ss, ii, ww, m, r, mv);

ケプラー方程式を呼び出す

var xx = new Array(new N6LVector(3));
var f = planet[i].kepler(nday, xx);

planet.js
//kepler//ケプラー方程式
N6LPlanet.prototype.kepler = function(nday, xx, num)

nday:1996/1/1からの経過日数(実数)
xx:座標を返す(N6LVector(3))
num:近似式繰り返し回数(デフォルト100)

↑で得た座標を日心黄道直交座標系に変換

planet.js
//xyz to ecliptic//xyz系を日心黄道直交座標系に変換
N6LPlanet.prototype.ecliptic = function(x, y, z, xyz)

これらの中身が気になる方はplanet.jsに潜ってください
参考文献は
天体の位置計算増補版長沢工著地人書館
理科年表
です

メインループ
solarsystem.js

...省略...

function UpdateFrameKepler() {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;

  var dat1;
  var day = time / msecPerDay;
  var tm = dt;
  if(dt != 0.0) {
    time = time + tm;
    var datt = dat.getTime();
    var dat1t = datt + time;
    var dat1 = new Date(dat1t);

    PlanetInit(dat1);    //新しい日時で惑星初期化

    setmp();
    setday(dat1);
  }
}

...省略...

function setmp() {
  var i;
  for(i = 0; i < planetnum; i++) {
    var elm = document.getElementById(IDTransA[i]);
    var sp = new x3dom.fields.SFVec3f(mp[i].x.x[1] / CNST_AU / Zoom, -mp[i].x.x[0] / CNST_AU / Zoom, mp[i].x.x[2] / CNST_AU / Zoom);
    elm.setAttribute('translation', sp.toString());
  }
}


PlanetInit(dat1); //新しい日時で惑星初期化
これで惑星の位置が更新されます

N6LPlanet,N6LMassPointクラスを使えば
ややこしい計算がカプセル化されているので
比較的容易な三次元座標計算位で
ケプラー方程式の惑星軌道計算が出来ます
相対性理論的ニュートン力学はN6LRngKtクラスに
カプセル化されています





3Dプログラミング入門講座
NAS6LIB
ポリゴンテスト解説・x3dom使い方
その1:ベクトル
その2:行列
その3:クォータニオン
その4:基本総復習・実践
その5:応用その1・メッシュアニメーション、動的テクスチャ
その6:応用その2・GLSL、カスタムシェーダー、キーボード、ファイル
その7:応用その3・ゲームプログラミング、タグの動的追加
その8:応用その4・GLSL、シェーダー、その2
その9:物理演算その1・電卓で相対性理論を解く方法
その10:物理演算その2・相対性理論的ニュートン力学
その11:物理演算その3・ケプラー方程式で惑星軌道シミュレーターを作る
その12:物理演算その4・ルンゲクッタ法で作った 相対性理論的ニュートン力学物理エンジンで惑星軌道シミュレーターを作る


 ■■■ THREE.JSプログラミング講座 ■■■ 
THREE.JSテスト解説・THREE.JS使い方
THREE.JS examplesをいじってみた(フレネル反射透過シェーダー)


<<prev 相対性理論的ニュートン力学 : ルンゲクッタ法で作った 相対性理論的ニュートン力学物理エンジンで惑星軌道シミュレーターを作る next>>




戻る