参考: パソコンで見る天体の動き, pp.166-175
原点に固定してある超重量物との距離が 0.65 よりも近くなったときだけ假想時間軸に移る。
TwoBodies = class() function TwoBodies:init(funcs, t0, inits, h, numOfDiv, inWhichDomain) self.unpack = unpack or table.unpack self.Sqrt = math.sqrt self.inWhichDomain = inWhichDomain self.numOfDiv = numOfDiv or 1 self.funcs = funcs self.step = h self.h = self.step / self.numOfDiv self.t0 = t0 self.inits = inits self.convert2Reals = function(r, inits) local u, uPrm, v, vPrm = inits[1], inits[2], inits[3], inits[4] local x = u^2 - v^2 local y = 2 * u * v local xDot = (2/r) * (u * uPrm - v * vPrm) local yDot = (2/r) * (v * uPrm + u * vPrm) return x, xDot, y, yDot end self.setInits = function() if self.inWhichDomain == "Real" then self.x, self.xDot, self.y, self.yDot = unpack(self.inits) self.tVirtual = 0 self.tReal = self.t0 self.r = self.Sqrt((self.x - 0)^2 + (self.y - 0)^2) self.initsReal = self.inits elseif self.inWhichDomain == "Virtual" then self.u, self.uPrm, self.v, self.vPrm, self.tReal = unpack(self.inits) self.tVirtual = self.t0 self.r = self.u^2 + self.v^2 self.x, self.xDot, self.y, self.yDot = self.convert2Reals(self.r, {self.u, self.uPrm, self.v, self.vPrm}) end end self.list2oneColMat = function(list) local oneColMat = {} for i = 1, #list do oneColMat[i] = {} oneColMat[i][1] = list[i] end return oneColMat end self.newMat = function(numRows, numCols, val) local mat = {} for i = 1, numRows do mat[i] = {} for j = 1, numCols do mat[i][j] = val or nil end end return mat end self.midPoints = function(funcs, t0, inits, h) local dim = #funcs local S = 7 local temp = {} local t = {} local x = self.newMat(dim, 1, _) local X = self.newMat(dim, 1, _) for s = 1, S do local N = 2^s local hs = h/N t[1] = t0 + hs for i = 1, dim do x[i][1] = inits[i] + hs * funcs[i](t0, self.unpack(inits)) end for n = 0, N - 1 do for i = 1, dim do temp[i] = x[i][n+1] end t[n+2] = (t[n] or t0) + 2 * hs for i = 1, dim do x[i][n+2] = (x[i][n] or inits[i]) + 2 * hs * funcs[i](t[n+1], self.unpack(temp)) end end for i = 1, dim do X[i][s] = (1/4) * (x[i][N-1] + 2 * x[i][N] + x[i][N+1]) end end return X end self.neville = function(listOfLists) local numOfLists = #listOfLists local numOfHs = #listOfLists[1] local X = {} local extrapo = {} for n = 1, numOfLists do X[n] = self.list2oneColMat(listOfLists[n]) for i = 2, numOfHs do for j = 2, i do X[n][i][j] = X[n][i][j-1] + (X[n][i][j-1] - X[n][i-1][j-1]) / (4^(j-1) - 1) end end extrapo[n] = X[n][numOfHs][numOfHs] end return extrapo end self.setInits() end function TwoBodies:updateEX() for i = 1, self.numOfDiv do self.inits = self.neville(self.midPoints(self.funcs, self.t0, self.inits, self.h)) self.t0 = self.t0 + self.h self.setInits() end return self end function TwoBodies:print() print(self.t0, table.concat(self.inits, "\t, ")) return self end --- 確かめる。 do local unpack = unpack or table.unpack local Sqrt = math.sqrt local Ceil = math.ceil local Ln = math.log local function sign(n) return (n >= 0) and 1 or -1 end local function convert2Virtuals(t, r, inits) local x, xDot, y, yDot = inits[1], inits[2], inits[3], inits[4] local u = Sqrt((r + x)/2) local v = Sqrt((r - x)/2) * sign(y) local uPrm = 0.5 * ( u * xDot + v * yDot) local vPrm = 0.5 * (-v * xDot + u * yDot) return u, uPrm, v, vPrm, t end local function varStore(t_real, t_virtual, x, y, r) var.store("t_real", t_real) var.store("t_virtual", t_virtual) var.store("x", x) var.store("y", y) var.store("r", r) end local function ceilB(num) if (num > 0) and (num < 1) then return 2^Ceil(Ln(num)/Ln(2)) else return Ceil(num) end end local e = 1 - (1e-14) -- 離心率 local m = {x = -(1 + e), y = 0, xDot = 0, yDot = -Sqrt(1 - e^2)/(1 + e)} -- 初期条件 local r = Sqrt((m.x - 0)^2 + (m.y - 0)^2) -- 初期条件での 2 体間の距離 local E = (m.xDot^2 + m.yDot^2)/2 - (1/r) -- 全エネルギー local R = 0.65 -- 2 体間の距離の閾値 local hReal = 0.0625 local hVirtual = ceilB(hReal/R) -- 現実時間軸で解く微分方程式 local function xDot(t, x, xDot, y, yDot) return xDot end local function xDotDot(t, x, xDot, y, yDot) return -x / (Sqrt(x^2 + y^2))^3 end local function yDot(t, x, xDot, y, yDot) return yDot end local function yDotDot(t, x, xDot, y, yDot) return -y / (Sqrt(x^2 + y^2))^3 end local funcsReal = {xDot, xDotDot, yDot, yDotDot} -- 假想時間軸で解く微分方程式 local function uPrm(t, u, uPrm, v, vPrm) return uPrm end local function uPrmPrm(t, u, uPrm, v, vPrm) return (E/2) * u end local function vPrm(t, u, uPrm, v, vPrm) return vPrm end local function vPrmPrm(t, u, uPrm, v, vPrm) return (E/2) * v end local function tPrm(t, u, uPrm, v, vPrm) return u^2 + v^2 end local funcsVirtual = {uPrm, uPrmPrm, vPrm, vPrmPrm, tPrm} local tb local function updateScreen() varStore(tb.tReal, tb.tVirtual, tb.x, tb.y, tb.r) end local function reset() local initsReal = {m.x, m.xDot, m.y, m.yDot} local initsVirtual = {convert2Virtuals(0, r, initsReal)} if r >= R then tb = TwoBodies(funcsReal , 0, initsReal , hReal , _, "Real" ) elseif r < R then tb = TwoBodies(funcsVirtual, 0, initsVirtual, hVirtual, _, "Virtual") end updateScreen() end local reset_off_on = nil function on.construction() reset() var.store("reset_off_on", 0) var.monitor("reset_off_on") end function on.varChange() reset_off_on = var.recall("reset_off_on") if reset_off_on > 0.7 then timer.start(tb.step) elseif reset_off_on < -0.7 then reset() else timer.stop() end end function on.timer() tb:updateEX() -- 更新するたびに if (tb.inWhichDomain == "Real") and (tb.r < R) then -- 今現実時間軸にいて、かつ 2 体間の距離が近づきすぎたら、 local initsVirtual = {convert2Virtuals(tb.tReal, tb.r, tb.initsReal)} -- 現実値を假想値に変換して、 tb = TwoBodies(funcsVirtual,0 ,initsVirtual , hVirtual, _, "Virtual") -- 假想時間軸で 2 体をインスタンス化し直すが、 elseif (tb.inWhichDomain == "Virtual") and (tb.r >= R) then -- 今假想時間軸にいて、かつ 2 体間の距離が十分に離れたら、 tb = TwoBodies(funcsReal, tb.tReal, {tb.x, tb.xDot, tb.y, tb.yDot}, hReal, _, "Real" ) -- 現実時間軸で 2 体をインスタンス化し直す。 end updateScreen() end end
