ピタゴラス3体問題4 of 4

function on.resize()
   W, H = platform.window:width(), platform.window:height()
   X0, Y0 = math.floor(W/2), math.floor(H/2)
   UNIT = math.floor(H/8)
   DT = {1/(2^4), 1/(2^18)}
   TOL = 1/10000
end
function on.paint(gc)
   local x31, vx31, y31, vy31, x41, vx41, y41, vy41, x51, y51, dt1, tol = 1, 0, 3, 0, -2, 0, -1, 0, 1, -1, DT[1], TOL   
   for i = 0, ★★★★★★★ do
      local dt2, x32, vx32, y32, vy32, x42, vx42, y42, vy42, x52, y52 = pyt(x31, vx31, y31, vy31, x41, vx41, y41, vy41, dt1, TOL)
      gc:setColorRGB(255, 0, 0)
      gc:drawLine(xy2screen(x31, y31, x32, y32, X0, Y0, UNIT))
      gc:setColorRGB(0, 255, 0)
      gc:drawLine(xy2screen(x41, y41, x42, y42, X0, Y0, UNIT))
      gc:setColorRGB(0, 0, 255)
      gc:drawLine(xy2screen(x51, y51, x52, y52, X0, Y0, UNIT))
      local r34 = math.sqrt(((x32 - x42)^2) + ((y32 - y42)^2))
      local r35 = (1/5) * math.sqrt(((8 * x32 + 4 * x42)^2) + ((8 * y32 + 4 * y42)^2))
      local r45 = (1/5) * math.sqrt(((3 * x32 + 9 * x42)^2) + ((3 * y32 + 9 * y42)^2))
      local dist = math.min(r34, r35, r45)
      if     0.5  <= dist then dt2 = DT[1]
      elseif dist <  0.5  then dt2 = DT[2]
      end
      dt1, x31, vx31, y31, vy31, x41, vx41, y41, vy41, x51, y51 = dt2, x32, vx32, y32, vy32, x42, vx42, y42, vy42, x52, y52
   end
end

function pyt(x31, vx31, y31, vy31, x41, vx41, y41, vy41, dt1, tol)
   var.store("x31", x31); var.store("vx31", vx31); var.store("y31", y31); var.store("vy31", vy31);
   var.store("x41", x41); var.store("vx41", vx41); var.store("y41", y41); var.store("vy41", vy41);
   var.store("dt1", dt1); var.store("tol", tol)
   local temp = math.eval("pyt(x31, vx31, y31, vy31, x41, vx41, y41, vy41, dt1, tol)")
   return  temp[1][2], -- dt2
           temp[2][2], -- x32
           temp[3][2], -- vx32
           temp[4][2], -- y32
           temp[5][2], -- vy32
           temp[6][2], -- x42
           temp[7][2], -- vx42
           temp[8][2], -- y42
           temp[9][2], -- vy42
           -(3 * temp[2][2] + 4 * temp[6][2])/5, -- x52
           -(3 * temp[4][2] + 4 * temp[8][2])/5  -- y52
end
----------------------------------------------
-- 2 組の xy 座標をスクリーン座標に変換する --
----------------------------------------------
function xy2screen(x1, y1, x2, y2, X0, Y0, UNIT)
   return (X0 + x1 * UNIT), (Y0 - y1 * UNIT), (X0 + x2 * UNIT), (Y0 - y2 * UNIT)
end