3Dプログラミング入門講座・その12:物理演算その4・ルンゲクッタ法で作った 相対性理論的ニュートン力学物理エンジンで惑星軌道シミュレーターを作る





・理論解説
その9:物理演算その1・電卓で相対性理論を解く方法
その10:物理演算その2・相対性理論的ニュートン力学
その11:物理演算その3・ケプラー方程式で惑星軌道シミュレーターを作る

その17:ケプラー方程式カプセルライブラリ使用法
その19:物理演算その5・重力多体問題
その22:物理エンジンライブラリ解説(ケプラー方程式・ルンゲクッタ・相対論的万有引力)


相対性理論衛星軌道計算シミュレーターフォーム by javascript.htm

相対性理論衛星軌道計算シミュレーターフォーム by javascript.zip



x3dom、自作ライブラリを用いた相対性理論的ニュートン力学の物理エンジンのjavascriptのサンプル

初期化
satellite.js

...省略...

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

  var radioList = document.getElementsByName("PUTSEL");
  var i;
  for(i = 0; i

N6LMassPoint構築
惑星初期化
PlanetInit(dat);
軌道線分計算
setline();

初期化実装
satellite.js

...省略...

//惑星初期化
function PlanetInit(dat) {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var i;
  var j;
    for(i = 0; i < planetnum; i++) {
      if(mp[i].mass < 0.0) continue;
      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 setline() {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var a = new Date(1996,6,1,0,0,0);
  var ndayR = 0;
  var i;
  var j;
  var k;
  var n = 32;
  var str;
  for(i = 0; i < planetnum; i++) {
    str = "";
    if(mp[i].mass <= 0.0) {
      var rt = 0.001;
      var base = 4;
      var ofs = i * 0.1;
      var x;
      var y; 
      for(x = base, y = base; -base < x ; x--) 
        str += Number((x + ofs) * rt).toString() + " " + Number((y + ofs) * rt).toString() + ", ";
      for(; -base < y ; y--) 
        str += Number((x + ofs) * rt).toString() + " " + Number((y + ofs) * rt).toString() + ", ";
      for(; x < base ; x++) 
        str += Number((x + ofs) * rt).toString() + " " + Number((y + ofs) * rt).toString() + ", ";
      for(; y < base ; y++) 
        str += Number((x + ofs) * rt).toString() + " " + Number((y + ofs) * rt).toString() + ", ";
      str += Number((x + ofs) * rt).toString() + " " + Number((y + ofs) * rt).toString();
      var elm;
      var sp;

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

      elm = document.getElementById(IDT[i]);
      sp = new x3dom.fields.SFVec4f(0, 1, 0, 0);
      elm.setAttribute('rotation', sp.toString());
      continue;
    }
    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 eps = Math.PI / 32;
      if(Math.PI - eps < f && f < Math.PI + eps) ndayR = nday;
      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 a;
    var b = 0;
    var c;
    for(a = 0; a < planetnum; a++)
      if(0.0 < mp[a].mass){b++;c=a;}
    if(b == 2 && i) {
      var msecPerMinute = 1000 * 60;
      var msecPerHour = msecPerMinute * 60;
      var msecPerDay = msecPerHour * 24;
      document.F1.myFormTIME.value = (ndayR + planet[c].m_dat0.getTime() / msecPerDay) / 365.2425;
      dat = new Date(ndayR * msecPerDay + planet[c].m_dat0.getTime());
    }
*/
    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());
  }
}

...省略...


二つのケプラー方程式から初期位置と初期速度を求め
日心黄道直交座標系に変換

メインループ実装
satellite.js
	
...省略...

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

  var dat1;
  var tm = Math.abs(Speed) * msecPerDay / 1000;
  var adt = Math.abs(dt);
  var t;
  var i;
  if(dt != 0.0) {
    for(t = adt; t <= tm; t += adt) {
      time = time + dt * 1000;
      //質点アップデート
      rk.UpdateFrame()

      //太陽原点補正
      for(i = 1; i < planetnum; i++) {
        rk.mp[i].x = rk.mp[i].x.Sub(rk.mp[0].x);
        mp[i].x = new N6LVector(rk.mp[i].x);
      }
      rk.mp[0].x = rk.mp[0].x.ZeroVec();
    }
    var datt = dat.getTime();
    var dat1t = datt + time;
    var dat1 = new Date(dat1t);
    setmp();
    settime();
  } 
}

function settime() {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var days = time / msecPerDay;
  document.F1.myFormTIME.value = days / 365.2425;
}

function setmp() {
  var i;
  for(i = 0; i < planetnum; i++) {
    if(mp[i].mass < 0.0) {
      var elm = document.getElementById(IDTransA[i]);
      var sp = new x3dom.fields.SFVec3f(0, 0, 0);
      elm.setAttribute('translation', sp.toString());
      continue;
    }
    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());
  }
}

function viewp() {
  if(!x3domRuntime) return;
  var selid = F1.VP.selectedIndex;
  var elm = document.getElementById('viewp000');

  var SWM = x3domRuntime.viewMatrix().inverse(); //ワールド回転行列取得
  var WM = new N6LMatrix().FromX3DOM(SWM);
  var Seye = SWM.multMatrixPnt(new x3dom.fields.SFVec3f(0, 0, 0)); //視点位置取得
  var sp = new x3dom.fields.SFVec3f(mp[selid].x.x[1] / CNST_AU / Zoom, -mp[selid].x.x[0] / CNST_AU / Zoom, mp[selid].x.x[2] / CNST_AU / Zoom);
  var Sat = x3dom.fields.SFVec3f.copy(sp);
  var lookat = new N6LVector([1.0, Sat.x, Sat.y, Sat.z], true);
  var LAM = WM.LookAtMat2(lookat);
  var Vec = LAM.Vector();
  var ori = Vec.ToX3DOM();

  elm.setAttribute('position', Seye.toString());
  elm.setAttribute('orientation', ori.toString());
  elm.setAttribute('centerOfRotation', sp.toString());
}

...省略...


rk.Init(pmp, dt);

N6LRngKtクラスをN6LMassPointクラスのArrayと
微時間dtで初期化したら

//質点アップデート
rk.UpdateFrame()

この関数は相対性理論的ニュートン力学の物理エンジンが
カプセル化されているのでただ呼び出すだけで
N体重力計算を行ってくれます

これをちょっと改変してロケットの軌道シミュレーターとか 作れないわけではないけれど
弾道ミサイルとかの計算に使われたら嫌だからそこまでは作らないし また作ったとしても公開はしません


rk.UpdateFrame2()

は簡略高速化でmp[0]つまり太陽とほかの質点とだけの重力関係の質点アップデータです






 ■■■ 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・ルンゲクッタ法で作った 相対性理論的ニュートン力学物理エンジンで惑星軌道シミュレーターを作る

その13:経路探索(A*:A-STAR)&巡回セールスマン問題 : 巨大サイズ : くろにゃんこ

その14:プログラミングにおける配列テーブルテクニック
その15:javascriptのクラス活用法
その16:透視射影公式テスト

その17:ケプラー方程式カプセルライブラリ使用法
その18:CSVファイル処理
その19:物理演算その5・重力多体問題
その20:同次座標について(3D座標系の基本の基本)
その21:おさらいコモンクラスの宣言
その22:物理エンジンライブラリ解説(ケプラー方程式・ルンゲクッタ・相対論的万有引力)
その23:サイト内のリンク先の更新確認がjavascriptで出来るようにする方法
その24:グレゴリオ暦日付計算


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

THREE.JS (半透明シェーダー)

THREE.JS 3D演算で必要な計算(具体例)★とても重要★
THREE.JS THREE-VRM をいじってみた



<<prev ケプラー方程式で惑星軌道シミュレーターを作る : 経路探索(A*:A-STAR) next>>




戻る