３Ｄプログラミング入門講座・その１１：物理演算その３・ケプラー方程式で惑星軌道シミュレーターを作る

．．．省略．．．

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

  var radioList = document.getElementsByName("CALC");
  for(i = 0; i < radioList.length; i++){
      if(radioList[i].checked){
          CalcWay = Number(radioList[i].value);
          break;
      }
  }

  bWaiting = false;
  if(!CalcWay) InitKepler();
  else InitRelative();
  onWaiting();
}

function onWaiting() {
  if(!bWaiting) {
    fFst = -1;
    setmp();
    setline();
    if(bBBB) bRunning = true;
  }
}

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

  if(fFst != 0) {
    fFst = 0;
    time = 0.0;
    if(CalcWay) dt = Speed * 60 * 60;
    else dt = Speed * msecPerDay;
  }
 
  //メインループ
  if(CalcWay) UpdateFrameRelative();
  else UpdateFrameKepler();
}

function getnow() {
  var nt = new Date(Number(F1.T1.value), Number(F1.T2.value) - 1, Number(F1.T3.value), Number(F1.T4.value), Number(F1.T5.value), Number(F1.T6.value));
  return nt;
}

function InitKepler() {
  if(!bRead) fFst = 1;
  dat = getnow();
  PlanetInit(dat);
}

．．．省略．．．

//惑星初期化
function PlanetInit(dat) {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var i;
  var j;
  if(0 < fFst) {
    //データファイル読み込み
    if(!bWaiting) {
      bWaiting = true;
      readCSV('./javascripts/nas6lib/PData000.txt', 'analyzeCSV', 'readedCSV');
    }
    return;
  }
  else {
    for(i = 0; i < planetnum; 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 readedCSV(res) {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  bWaiting = false;
  bRead = true;

  for(i = 0; i < planetnum; i++) {

    var PlanetName = res[i][0];
    var PlanetNo = Number(res[i][1]);
    var EpochYY = Number(res[i][2]);
    var EpochMM = Number(res[i][3]);
    var EpochDD = Number(res[i][4]);
    var Epochh = Number(res[i][5]);
    var Epochm = Number(res[i][6]);
    var Epochs = Number(res[i][7]);
    var a = Number(res[i][8]);
    var e = Number(res[i][9]);
    var m0 = Number(res[i][10]);
    var npd = Number(res[i][11]);
    var ra = Number(res[i][12]);
    var rb = Number(res[i][13]);
    var p = Number(res[i][14]);
    var ss = Number(res[i][15]);
    var ii = Number(res[i][16]);
    var ww = Number(res[i][17]);
    var m = Number(res[i][18]);
    var r = Number(res[i][19]);
    var mv = Number(res[i][20]);

    var dat0 = new Date(EpochYY, EpochMM - 1, EpochDD, Epochh, Epochm, Epochs);
    var datt = dat.getTime();
    var dat0t = dat0.getTime();

    var ddat = (datt - dat0t) / msecPerDay;

    //惑星初期化
    planet[i] = new N6LPlanet();
    planet[i].Create(PlanetNo, PlanetName, ddat, dat0, a, e, m0, npd, ra, rb, p, ss, ii, ww, m, r, mv);

    //質点初期化
    mp[i] = new N6LMassPoint(planet[i].x0, planet[i].v0, m, r, e);
  }
  return true;
}

//惑星軌道線分設定
function setline() {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;
  var i;
  var j;
  var k;
  var n = 32;
  var str;
  for(i = 1; i < planetnum; i++) {
    str = "";
    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);
      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 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 UpdateFrameKepler() {
  var msecPerMinute = 1000 * 60;
  var msecPerHour = msecPerMinute * 60;
  var msecPerDay = msecPerHour * 24;

  var dat1;
  var day = time / msecPerDay;
  var tm = dt;
  if(dt != 0.0) {
    time = time + tm;
    var datt = dat.getTime();
    var dat1t = datt + time;
    var dat1 = new Date(dat1t);

    PlanetInit(dat1);    //新しい日時で惑星初期化

    setmp();
    setday(dat1);
  }
}

．．．省略．．．

function setmp() {
  var i;
  for(i = 0; i ＜ planetnum; i++) {
    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());
  }
}