算額（その1116）

算額（その1116）

四 岩手県花巻市南笹間 東光寺 慶応2年(1866)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：球，3次元

盤の上に大球を載せ，その周りに小球を数個互いに接するように並べる。大球の直径と小球の個数が与えられたとき，小球の直径を求めよ。

注：単に大球と書いているが，図にあるように，半球である

大球の半径と中心座標を R, (0, 0, 0)
小球の半径と中心座標を r, (x, 0, r), (x*cos(2π/n), x*sin(2π/n)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms n::positiveinteger, R::positive, r::positive,
     x::positive
eq1 = x^2 + r^2 - (R + r)^2
eq2 = x*sind(Sym(180)/n) - r;

(r, x) = solve([eq1, eq2], (r, x))[2];
r |> println
x |> println

   R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))*sin(pi/n)
   R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))

小球の半径は大球の半径の (sqrt(sin(pi/n)^2 + 1) + sin(pi/n))*sin(pi/n) 倍である。

n = 12
R = 5
r = R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))*sin(pi/n)
x = R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))
@printf("大球の直径が %g で，%i 個の小球を並べるとき，小球の直径は %g である。\n小球の中心座標は以下の通りである。\n", 2R, n, 2r)
for i = 1:n
   println((x*cosd(360*i/n), x*sind(360*i/n)))
end

   大球の直径が 10 で，12 個の小球を並べるとき，小球の直径は 3.34335 である。
   小球の中心座標は以下の通りである。
   (5.593527388798881, 3.229424543642579)
   (3.229424543642579, 5.593527388798881)
   (0.0, 6.458849087285158)
   (-3.229424543642579, 5.593527388798881)
   (-5.593527388798881, 3.229424543642579)
   (-6.458849087285158, 0.0)
   (-5.593527388798881, -3.229424543642579)
   (-3.229424543642579, -5.593527388798881)
   (0.0, -6.458849087285158)
   (3.229424543642579, -5.593527388798881)
   (5.593527388798881, -3.229424543642579)
   (6.458849087285158, 0.0)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   n = 12
   R = 5
   r = R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))*sin(pi/n)
   x = R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))
   plot(showaxis=false)
   circle(0, 0, R, :blue, beginangle=0, endangle=180)
   circle(x, r, r)
   segment(x, R, x, -3R, :green)
   segment(R, R, R, -3R, :blue)
   segment(x + r, R, x + r, -3R, :red)
   segment(x - r, R, x - r, -3R, :red)
   segment(0, 0, x, r)
   segment(x, 0, x, r, :green)
   segment(-R, 0, R + 2r, 0)
   base = -1.7R
   circle(0, base, R, :blue)
   circle(0, base, x, :green)
   for i = 0:n - 1
       (x0, y0) = (x*cosd(i*360/n), base + x*sind(i*360/n))
       circle(x0, y0, r)
   end
   i = 0.5
   (x0, y0) = (x*cosd(i*360/n), base + x*sind(i*360/n))
   segment(0, base, x + r, base)
   segment(0, base, x, y0)
   segment(x, base, x0, y0)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R, 0, " R", :blue, :left, delta=-delta)
       point(x, r, " (x,r)", :red, :left, :vcenter)
       circle(0, base, 5delta, beginangle=0, endangle=360/2n)
       circle(0, base, 5.3delta, beginangle=0, endangle=360/2n)
       point(5delta, base + 1.5delta, "θ=π/n", mark=false)
       point(3delta, -20delta, "上から見た図", mark=false)
       point(3delta, 8delta, "横から見た図", mark=false)
       plot!(xlims=(-5delta, x + r + 5delta), ylims=(base - 8delta, R + 2delta))
   end
end;