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

．．．省略．．．

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

．．．省略．．．

//惑星初期化
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;
    //惑星１周を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());
  }
}

．．．省略．．．

	
．．．省略．．．

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());
}

．．．省略．．．