３Ｄプログラミング入門講座・その１９：物理演算その５・重力多体問題





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






　　




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;
    //惑星１周を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()で質点の相互重力計算をしています
こんな感じで重力多体問題を解きます





　■■■　３Ｄプログラミング入門講座　■■■　
ＮＡＳ６ＬＩＢ
ポリゴンテスト解説・x3dom使い方
その１：ベクトル
その２：行列
その３：クォータニオン
その４：基本総復習・実践
その５：応用その１・メッシュアニメーション、動的テクスチャ
その６：応用その２・ＧＬＳＬ、カスタムシェーダー、キーボード、ファイル
その７：応用その３・ゲームプログラミング、タグの動的追加
その８：応用その４・ＧＬＳＬ、シェーダー、その２
その９：物理演算その１・電卓で相対性理論を解く方法
その１０：物理演算その２・相対性理論的ニュートン力学
その１１：物理演算その３・ケプラー方程式で惑星軌道シミュレーターを作る

その１２：物理演算その４・ルンゲクッタ法で作った 相対性理論的ニュートン力学物理エンジンで惑星軌道シミュレーターを作る

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

その１４：プログラミングにおける配列テーブルテクニック
その１５：javascriptのクラス活用法
その１６：透視射影公式テスト

その１７：ケプラー方程式カプセルライブラリ使用法
その１８：ＣＳＶファイル処理
その１９：物理演算その５・重力多体問題
その２０：同次座標について（３Ｄ座標系の基本の基本）
その２１：おさらいコモンクラスの宣言
その２２：物理エンジンライブラリ解説(ケプラー方程式・ルンゲクッタ・相対論的万有引力)
その２３：サイト内のリンク先の更新確認がjavascriptで出来るようにする方法
その２４：グレゴリオ暦日付計算


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

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

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



＜＜prev　ＣＳＶファイル処理　：　同次座標について　next＞＞


戻る