3Dプログラミング入門講座・その19:物理演算その5・重力多体問題






物理演算その5・重力多体問題






  




var TMan = new N6LTimerMan();  //タイマーマネージャー

var CNST_AU = 1.49597870700e+11;
var fFst = 1;
var dat;
var time;
var dt;
var Speed = 500000000000000.0;

var planet = new Array();
var mp = new Array();
var rk = new N6LRngKt();

var pnt = new Array(33);
var pnt2 = new Array(33);

var nzoom = 1.0;
var nz = 1.25;


function enter(){
  init();
  TMan.add();
  TMan.timer[0].setalerm(function() { GLoop(0); }, 50);  //メインループセット
}

//メインループ
function GLoop(id){
  onRunning();

  TMan.timer[id].setalerm(function() { GLoop(id); }, 50);  //メインループ再セット
}


function init() {
  nzoom = 1.0;
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;

  time = 0.0;
  dt = Speed * 60 * 60;

  dat = new Date();
  PlanetInit(dat);
  dt = Speed * 60 * 60;
  var pmp = new Array();
  var i;
  for(i = 0; i < 3; i++) pmp[i] = new N6LMassPoint(mp[i]);
  rk.Init(pmp, dt);
  calcline();

}

//惑星初期化
function PlanetInit(dat) {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var i;
  var j;

    //惑星初期化
    planet[0] = new N6LPlanet();
    planet[0].Create(0, 'P0', new Date(), new Date(), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 100, 1, 1);
    planet[1] = new N6LPlanet();
    planet[1].Create(1, 'P1', new Date(), new Date(), 1, 0, 0, 7.022e-15, 1, 1, 140361930687886.16, 0, 0, 0, 1, 1, 1);
    planet[2] = new N6LPlanet();
    planet[2].Create(2, 'P2', new Date(), new Date(), 0.6, 0.666, 0, 1.51e-14, 0.2, 1, 65234330399484.34, 0, 0, 0, 1, 1, 1);

    //質点初期化
    mp[0] = new N6LMassPoint(planet[0].x0, planet[0].v0, 100, 1, 1);
    mp[1] = new N6LMassPoint(planet[1].x0, planet[1].v0, 1, 1, planet[1].m_e);
    mp[2] = new N6LMassPoint(planet[2].x0, planet[2].v0, 1, 1, planet[2].m_e);

    planet[0].x0 = new N6LVector(3).ZeroVec();
    planet[0].v0 = new N6LVector(3).ZeroVec();
    mp[0] = new N6LMassPoint(planet[0].x0, planet[0].v0, planet[0].m_m, planet[0].m_r, planet[0].m_e);
    for(i = 1; i < 3; 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, planet[i].m_e);
  }
}

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

  //メインループ
  UpdateFrameRelative();
}


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 < 3; 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();
    }
    setmp();
  } 
  //if((1.8 * CNST_AU < mp[1].x.Abs())||(1.8 * CNST_AU < mp[2].x.Abs())){
    //init();
    //alert("Reset.");
  //}
}

function setmp() {
 //描画コンテキストの取得
 var canvas = document.getElementById('cnv');
 if (canvas.getContext) {
  var z = 1 / 200;
  var context = canvas.getContext('2d');
  context.fillStyle = 'rgb(0,0,0)';
  context.fillRect(0,0,600,600);
  context.lineWidth = 3;
  context.strokeStyle = 'rgb(192,192,0)'; 
  context.fillStyle = 'rgb(192,192,0)';
  context.beginPath();
  context.arc((mp[0].x.x[0] / CNST_AU / z) * nzoom + 300, (-mp[0].x.x[1] / CNST_AU / z) * nzoom + 300, 30, 0, Math.PI*2, false);
  context.fill();
  context.closePath();
  context.strokeStyle = 'rgb(0,192,192)'; 
  context.fillStyle = 'rgb(0,192,192)';
  context.beginPath();
  context.arc((mp[1].x.x[0] / CNST_AU / z) * nzoom + 300, (-mp[1].x.x[1] / CNST_AU / z) * nzoom + 300, 10, 0, Math.PI*2, false);
  context.fill();
  context.closePath();
  context.strokeStyle = 'rgb(0,0,192)'; 
  context.fillStyle = 'rgb(0,0,192)';
  context.beginPath();
  context.arc((mp[2].x.x[0] / CNST_AU / z) * nzoom + 300, (-mp[2].x.x[1] / CNST_AU / z) * nzoom + 300, 10, 0, Math.PI*2, false);
  context.fill();
  context.closePath();

  var i;
  context.lineWidth = 1;
  context.strokeStyle = 'rgb(128,255,192)'; 
  context.fillStyle = 'rgb(128,255,192)';
  context.beginPath();
  context.moveTo((pnt[0].x[0] / CNST_AU / z) * nzoom + 300, (pnt[0].x[1] / CNST_AU / z) * nzoom + 300);
   for(i = 0; i < 33; i++){
    context.lineTo((pnt[i].x[0] / CNST_AU / z) * nzoom + 300, (pnt[i].x[1] / CNST_AU / z) * nzoom + 300);
  }
  context.closePath();
  context.stroke();


  context.lineWidth = 1;
  context.strokeStyle = 'rgb(0,192,0)'; 
  context.fillStyle = 'rgb(0,192,0)';
  context.beginPath();
  context.beginPath();
  context.moveTo((pnt2[0].x[0] / CNST_AU / z) * nzoom + 300, (pnt2[0].x[1] / CNST_AU / z) * nzoom + 300);
   for(i = 0; i < 33; i++){
    context.lineTo((pnt2[i].x[0] / CNST_AU / z) * nzoom + 300, (pnt2[i].x[1] / CNST_AU / z) * nzoom + 300);
  }
  context.closePath();
  context.stroke();
 }
}




function calcline() {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var i;
  var j;
  var k;
  var n = 32;
  for(i = 1; i < 3; i++){
    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);
      if(i == 1){
        pnt[j] = new N6LVector(3);
        pnt[j].x[0] = x1.x[0];
        pnt[j].x[1] = -x1.x[1];
        pnt[j].x[2] = x1.x[2];
      }
      else{
        pnt2[j] = new N6LVector(3);
        pnt2[j].x[0] = x1.x[0];
        pnt2[j].x[1] = -x1.x[1];
        pnt2[j].x[2] = x1.x[2];
      }
    }
    if(i == 1){
      pnt[j] = new N6LVector(3);
      pnt[j].x[0] = x0.x[0];
      pnt[j].x[1] = -x0.x[1];
      pnt[j].x[2] = x0.x[2];
    }
    else{
      pnt2[j] = new N6LVector(3);
      pnt2[j].x[0] = x0.x[0];
      pnt2[j].x[1] = -x0.x[1];
      pnt2[j].x[2] = x0.x[2];
    }
  }
}


function zoom(sw) {
  if(sw < 0) {
    nzoom /= nz;
  }
  else {
    nzoom *= nz;
  }
}




初期設定



function init() {
  nzoom = 1.0;
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;

  time = 0.0;
  dt = Speed * 60 * 60;

  dat = new Date();
  PlanetInit(dat);
  dt = Speed * 60 * 60;
  var pmp = new Array();
  var i;
  for(i = 0; i < 3; i++) pmp[i] = new N6LMassPoint(mp[i]);
  rk.Init(pmp, dt);
  calcline();

}

//惑星初期化
function PlanetInit(dat) {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var i;
  var j;

    //惑星初期化
    planet[0] = new N6LPlanet();
    planet[0].Create(0, 'P0', new Date(), new Date(), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 100, 1, 1);
    planet[1] = new N6LPlanet();
    planet[1].Create(1, 'P1', new Date(), new Date(), 1, 0, 0, 7.022e-15, 1, 1, 140361930687886.16, 0, 0, 0, 1, 1, 1);
    planet[2] = new N6LPlanet();
    planet[2].Create(2, 'P2', new Date(), new Date(), 0.6, 0.666, 0, 1.51e-14, 0.2, 1, 65234330399484.34, 0, 0, 0, 1, 1, 1);

    //質点初期化
    mp[0] = new N6LMassPoint(planet[0].x0, planet[0].v0, 100, 1, 1);
    mp[1] = new N6LMassPoint(planet[1].x0, planet[1].v0, 1, 1, planet[1].m_e);
    mp[2] = new N6LMassPoint(planet[2].x0, planet[2].v0, 1, 1, planet[2].m_e);

    planet[0].x0 = new N6LVector(3).ZeroVec();
    planet[0].v0 = new N6LVector(3).ZeroVec();
    mp[0] = new N6LMassPoint(planet[0].x0, planet[0].v0, planet[0].m_m, planet[0].m_r, planet[0].m_e);
    for(i = 1; i < 3; 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, planet[i].m_e);
  }
}



PlanetInit()ではplanet[]をN6LPlanet.Create()で構築しています
物理演算のRngKtで利用するためnew N6LMassPoint()も構築します
基準時間と微笑増加時間でplanet[i].keplerをして初期位置、速度を求めます
それを最終的にN6LMassPointに入れます


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 < 3; 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();
    }
    setmp();
  } 
  //if((1.8 * CNST_AU < mp[1].x.Abs())||(1.8 * CNST_AU < mp[2].x.Abs())){
    //init();
    //alert("Reset.");
  }
}


N6LRngKt.UpdateFrame()で質点の相互重力計算をしています
こんな感じで重力多体問題を解きます






 ■■■ 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:おさらいコモンクラスの宣言


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

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

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



<<prev CSVファイル処理 : 同次座標について next>>






戻る