物理演算その５・重力多体問題をマネージドクラス化する

この重力三体問題は黄色の重力源を小さな楕円を周回する青色の周回に合わせて
水色が青色で重力カタパルトにより円軌道がかき乱されてそのうち立場が入れ替わり
水色の周回に合わせて青色が重力カタパルトにより高度を得ていくという
万有引力シナリオです
さらにその先も立場が入れ替わるかもしれません

　　




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

const DataInstanceTemplate = Object({
  variablename: "PhysicsEngineData",
  thisproj: {
    CNST_AU: 1.49597870700e+11,
    fFst: 1,
    dat: 0,
    time: 0,
    dt: 0,
    Speed: 500000000000000.0,
    nzoom: 1.0,
    nz: 1.25,
    rk : null,
    planets : [],
    massPoints : [],
    Points : new Array(33),
    Points2 : new Array(33)
  }
});

var DataInstance = new N6LManagedClass(DataInstanceTemplate);

var proj;
var rk;
var planets;
var massPoints;
var Points;
var Points2;

/*
var rk = new N6LRngKt();
var planets = [];
var massPoints = [];
var Points = new Array(33);
var Points2 = new Array(33);
*/

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() {
  // 参照を一度取得する（これだけで短縮可能）
  proj = DataInstance.property.thisproj;
  rk = proj.rk;
  rk = new N6LRngKt();
  planets = proj.planets;
  massPoints = proj.massPoints;
  Points = proj.Points;
  Points2 = proj.Points2;

  proj.nzoom = 1.0;
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;

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

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

}

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

    //惑星初期化
    planets[0] = new N6LPlanet();
    planets[0].Create(0, 'P0', new Date(), new Date(), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 100, 1, 1);
    planets[1] = new N6LPlanet();
    planets[1].Create(1, 'P1', new Date(), new Date(), 1, 0, 0, 7.022e-15, 1, 1, 140361930687886.16, 0, 0, 0, 1, 1, 1);
    planets[2] = new N6LPlanet();
    planets[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);

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

    planets[0].x0 = new N6LVector(3).ZeroVec();
    planets[0].v0 = new N6LVector(3).ZeroVec();
    massPoints[0] = new N6LMassPoint(planets[0].x0, planets[0].v0, planets[0].m_m, planets[0].m_r, planets[0].m_e);
    for(i = 1; i < 3; i++){
      var dat0 = planets[i].m_dat0;
      var datt = proj.dat.getTime();
      var dat0t = dat0.getTime();
      var ddat = (datt - dat0t) / msecPerDay;
      var nday = ddat;

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

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

      massPoints[i] = new N6LMassPoint(planets[i].x0, planets[i].v0, planets[i].m_m, planets[i].m_r, planets[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(proj.Speed) * msecPerDay / 1000;
  var adt = Math.abs(proj.dt);
  var t;
  var i;

  if(proj.dt != 0.0) {
    for(t = adt; t <= tm; t += adt) {
      proj.time = proj.time + proj.dt * 1000;
      //質点アップデート
      rk.UpdateFrame();

      //太陽原点補正
      for(i = 1; i < 3; i++) {
        rk.mp[i].x = rk.mp[i].x.Sub(rk.mp[0].x);
        massPoints[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((massPoints[0].x.x[0] / proj.CNST_AU / z) * proj.nzoom + 300, (-massPoints[0].x.x[1] / proj.CNST_AU / z) * proj.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((massPoints[1].x.x[0] / proj.CNST_AU / z) * proj.nzoom + 300, (-massPoints[1].x.x[1] / proj.CNST_AU / z) * proj.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((massPoints[2].x.x[0] / proj.CNST_AU / z) * proj.nzoom + 300, (-massPoints[2].x.x[1] / proj.CNST_AU / z) * proj.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((Points[0].x[0] / proj.CNST_AU / z) * proj.nzoom + 300, (Points[0].x[1] / proj.CNST_AU / z) * proj.nzoom + 300);
   for(i = 0; i < 33; i++){
    context.lineTo((Points[i].x[0] / proj.CNST_AU / z) * proj.nzoom + 300, (Points[i].x[1] / proj.CNST_AU / z) * proj.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((Points2[0].x[0] / proj.CNST_AU / z) * proj.nzoom + 300, (Points2[0].x[1] / proj.CNST_AU / z) * proj.nzoom + 300);
   for(i = 0; i < 33; i++){
    context.lineTo((Points2[i].x[0] / proj.CNST_AU / z) * proj.nzoom + 300, (Points2[i].x[1] / proj.CNST_AU / z) * proj.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 / planets[i].m_nperday) * (j / n);
      var days = (proj.dat.getTime() - planets[i].m_dat0.getTime()) / msecPerDay;
      var nday = days + ad;
      var xx = new Array(new N6LVector(3));
      var vv = new Array(new N6LVector(3));
      var f = planets[i].kepler(nday, xx, vv);
      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){
        Points[j] = new N6LVector(3);
        Points[j].x[0] = x1.x[0];
        Points[j].x[1] = -x1.x[1];
        Points[j].x[2] = x1.x[2];
      }
      else{
        Points2[j] = new N6LVector(3);
        Points2[j].x[0] = x1.x[0];
        Points2[j].x[1] = -x1.x[1];
        Points2[j].x[2] = x1.x[2];
      }
    }
    if(i == 1){
      Points[j] = new N6LVector(3);
      Points[j].x[0] = x0.x[0];
      Points[j].x[1] = -x0.x[1];
      Points[j].x[2] = x0.x[2];
    }
    else{
      Points2[j] = new N6LVector(3);
      Points2[j].x[0] = x0.x[0];
      Points2[j].x[1] = -x0.x[1];
      Points2[j].x[2] = x0.x[2];
    }
  }
}


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

初期設定


function init() {
  // 参照を一度取得する（これだけで短縮可能）
  proj = DataInstance.property.thisproj;
  rk = proj.rk;
  rk = new N6LRngKt();
  planets = proj.planets;
  massPoints = proj.massPoints;
  Points = proj.Points;
  Points2 = proj.Points2;

  proj.nzoom = 1.0;
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;

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

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

}

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

    //惑星初期化
    planets[0] = new N6LPlanet();
    planets[0].Create(0, 'P0', new Date(), new Date(), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 100, 1, 1);
    planets[1] = new N6LPlanet();
    planets[1].Create(1, 'P1', new Date(), new Date(), 1, 0, 0, 7.022e-15, 1, 1, 140361930687886.16, 0, 0, 0, 1, 1, 1);
    planets[2] = new N6LPlanet();
    planets[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);

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

    planets[0].x0 = new N6LVector(3).ZeroVec();
    planets[0].v0 = new N6LVector(3).ZeroVec();
    massPoints[0] = new N6LMassPoint(planets[0].x0, planets[0].v0, planets[0].m_m, planets[0].m_r, planets[0].m_e);
    for(i = 1; i < 3; i++){
      var dat0 = planets[i].m_dat0;
      var datt = proj.dat.getTime();
      var dat0t = dat0.getTime();
      var ddat = (datt - dat0t) / msecPerDay;
      var nday = ddat;

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

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

      massPoints[i] = new N6LMassPoint(planets[i].x0, planets[i].v0, planets[i].m_m, planets[i].m_r, planets[i].m_e);
  }
}

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


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

  var dat1;
  var tm = Math.abs(proj.Speed) * msecPerDay / 1000;
  var adt = Math.abs(proj.dt);
  var t;
  var i;

  if(proj.dt != 0.0) {
    for(t = adt; t <= tm; t += adt) {
      proj.time = proj.time + proj.dt * 1000;
      //質点アップデート
      rk.UpdateFrame();

      //太陽原点補正
      for(i = 1; i < 3; i++) {
        rk.mp[i].x = rk.mp[i].x.Sub(rk.mp[0].x);
        massPoints[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()で質点の相互重力計算をしています
こんな感じで重力多体問題を解きます

さて今回の目玉のマネージドクラス化です

マネージドクラス導入のメリット
名前空間の保護: DataInstance という一つの窓口に全データを集約することで
グローバル変数の散乱を防ぎます。

状態の一元管理: シミュレーションのリセットや保存を行う際
DataInstance の中身を操作するだけで済むため
バグの混入リスクを劇的に減らせます。

コードのクリーン化: init() で参照を一度キャッシュ（代入）することで
処理コード自体は従来の書き方のまま、管理の安全性だけを向上させることが可能です。

【注意点】参照のライフサイクルについて
proj = DataInstance.property.thisproj; とすることで変数 proj に参照を渡していますが
DataInstance 側でオブジェクトの再生成（プロパティの入れ替え等）を行った場合は
再度 init() を呼ぶなどして参照を更新する必要があります。

リンク不良を避けるため、「データ構造の再構築時は、必ず参照変数も更新する」
という原則を心がけてください。

そして実践です
グローバルスコープにデータインスタンスのひな型と
実体とアクセス指定子を宣言します
var DataInstance = new N6LManagedClass(DataInstanceTemplate);
でN6LManagedClassでデータ実体を宣言しています


const DataInstanceTemplate = Object({
  variablename: "PhysicsEngineData",
  thisproj: {
    CNST_AU: 1.49597870700e+11,
    fFst: 1,
    dat: 0,
    time: 0,
    dt: 0,
    Speed: 500000000000000.0,
    nzoom: 1.0,
    nz: 1.25,
    rk : null,
    planets : [],
    massPoints : [],
    Points : new Array(33),
    Points2 : new Array(33)
  }
});

var DataInstance = new N6LManagedClass(DataInstanceTemplate);

var proj;
var rk;
var planets;
var massPoints;
var Points;
var Points2;

そしてinit()においてこのようにアクセス指定子の参照を設定して
以後このアクセス指定子でアクセスすればOKです
ただし参照切れ/リンク不良には注意が必要です


function init() {
  // 参照を一度取得する（これだけで短縮可能）
  proj = DataInstance.property.thisproj;
  rk = proj.rk;
  rk = new N6LRngKt();
  planets = proj.planets;
  massPoints = proj.massPoints;
  Points = proj.Points;
  Points2 = proj.Points2;

//・・・以下略
}

＜＜prev　分散自己安定化リレー制御[仮想アカウント版]　：　シュワルツシルト解の利用法　next＞＞

３Ｄプログラミング入門講座・その３０：物理演算その５・重力多体問題をマネージドクラス化する

物理演算その５・重力多体問題をマネージドクラス化する

さて今回の目玉のマネージドクラス化です